| 1 |
************************************************************************* |
| 2 |
* This version of momanhough performs tracker CALIBRATIONS: |
| 3 |
* - eta functions |
| 4 |
* - charge correlation |
| 5 |
* - ADC-to-MIP conversion |
| 6 |
* ------------------------------------------------------------------- |
| 7 |
* It works as it follow: |
| 8 |
* - it reads cluster information (read a list of LEVEL1 ntuples, |
| 9 |
* output of reduction.f) |
| 10 |
* - it performs track identification, select good tracks (6 hit planes) |
| 11 |
* and perform first-order fit |
| 12 |
* - it calculate calibration variables |
| 13 |
* - it creates a LEVEL2 ntuple with selected events (in order to check) |
| 14 |
* |
| 15 |
************************************************************************* |
| 16 |
program momanhough |
| 17 |
|
| 18 |
include '../common/commontracker.f' |
| 19 |
include '../common/common_momanhough.f' |
| 20 |
include '../common/momanhough_init.f' |
| 21 |
include '../common/common_level2debug.f' |
| 22 |
include '../common/common_mech.f' |
| 23 |
include '../common/common_xyzPAM.f' |
| 24 |
include '../common/common_mini_2.f' |
| 25 |
include '../common/calib.f' |
| 26 |
include '../common/level1.f' |
| 27 |
include '../common/level2.f' |
| 28 |
include '../common/common_calibration.f' |
| 29 |
|
| 30 |
* flag to set debug mode |
| 31 |
logical DEBUG |
| 32 |
common/dbg/DEBUG |
| 33 |
logical DONTDO |
| 34 |
data DONTDO/.true./ |
| 35 |
c$$$* flag to chose PFA |
| 36 |
c$$$ character*10 PFA |
| 37 |
c$$$ common/FINALPFA/PFA |
| 38 |
|
| 39 |
|
| 40 |
c parameter (inf=1.e8) !just a huge number... |
| 41 |
|
| 42 |
* external functions |
| 43 |
external npl |
| 44 |
external acoordsi,coordsi,nld,coord,dcoord |
| 45 |
c------------------------------------------------------------------------ |
| 46 |
c |
| 47 |
c local variables |
| 48 |
c |
| 49 |
c------------------------------------------------------------------------ |
| 50 |
|
| 51 |
character*24 processing_date |
| 52 |
|
| 53 |
parameter (nfile_max=500) |
| 54 |
parameter (lun_data_level1=71) !data file id number |
| 55 |
parameter (lun_data_level2=72) !data file id number |
| 56 |
parameter (lun_data_calib=74) !data file id number |
| 57 |
c parameter (lun_file_out=75) !data file id number |
| 58 |
|
| 59 |
character*74 data_file !data file name |
| 60 |
character*74 data_dir ! |
| 61 |
character*74 out_dir ! |
| 62 |
character*74 data_file_calib |
| 63 |
character*74 data_file_level1 |
| 64 |
character*74 data_file_level2 |
| 65 |
character*10 input_string(nfile_max) |
| 66 |
character*74 file_out |
| 67 |
character*74 fname_param |
| 68 |
character*7 rzdir |
| 69 |
|
| 70 |
|
| 71 |
|
| 72 |
# ifndef TEST2003 |
| 73 |
parameter(ncalibmax=50) |
| 74 |
character*40 file_calib(ncalibmax) |
| 75 |
parameter(lun_calib_list=66) !calibration list file id |
| 76 |
integer which_calib_last |
| 77 |
# endif |
| 78 |
|
| 79 |
integer tottr |
| 80 |
integer minevent !first event to be analysed |
| 81 |
logical FIMAGE ! |
| 82 |
|
| 83 |
|
| 84 |
COMMON/QUEST/IQUEST(100) |
| 85 |
c !permette di ottenere ntuple funzionanti nonostante |
| 86 |
c ! il messaggio dei 64K di RZOUT...!??? |
| 87 |
|
| 88 |
# ifdef TEST2003 |
| 89 |
parameter (nang=3) !number of angular bins |
| 90 |
parameter (angstep=4.) |
| 91 |
parameter (angmin=-6.) |
| 92 |
print*,' ' |
| 93 |
print*,'Set angular binning TEST 2003' |
| 94 |
print*,' ' |
| 95 |
|
| 96 |
nangbin=nang |
| 97 |
do ia=1,nangbin |
| 98 |
angL(ia)=angmin+angstep*(ia-1) |
| 99 |
angR(ia)=angmin+angstep*ia |
| 100 |
print*,ia,' - ',angL(ia),angR(ia) |
| 101 |
enddo |
| 102 |
print*,' ' |
| 103 |
|
| 104 |
|
| 105 |
# else |
| 106 |
c$$$ parameter (nang=21) !number of angular bins |
| 107 |
c$$$ parameter (angstep=2.) |
| 108 |
c$$$ parameter (angmin=-21.) |
| 109 |
c$$$ print*,' ' |
| 110 |
c$$$ print*,'Set angular binning Rome-2005' |
| 111 |
c$$$ print*,' ' |
| 112 |
c$$$ |
| 113 |
c$$$ nangbin=nang |
| 114 |
c$$$ do ia=1,nangbin |
| 115 |
c$$$ angL(ia)=angmin+angstep*(ia-1) |
| 116 |
c$$$ angR(ia)=angmin+angstep*ia |
| 117 |
c$$$ print*,ia,' - ',angL(ia),angR(ia) |
| 118 |
c$$$ enddo |
| 119 |
c$$$ print*,' ' |
| 120 |
|
| 121 |
parameter (nang=17) !number of angular bins |
| 122 |
* NB angL and angR are created with length 21 |
| 123 |
data angL |
| 124 |
$ /-21.,-17.,-13.,-11.,-9.,-7.,-5.,-3.,-1. |
| 125 |
$ ,1.,3.,5.,7.,9.,11.,13.,17.,0.,0.,0.,0./ |
| 126 |
data angR |
| 127 |
$ /-17.,-13.,-11.,-9.,-7.,-5.,-3.,-1. |
| 128 |
$ ,1.,3.,5.,7.,9.,11.,13.,17.,21.,0.,0.,0.,0./ |
| 129 |
print*,' ' |
| 130 |
print*,'Set angular binning Rome-2005' |
| 131 |
print*,' ' |
| 132 |
nangbin=nang |
| 133 |
do ia=1,nangbin |
| 134 |
print*,ia,' - ',angL(ia),angR(ia) |
| 135 |
enddo |
| 136 |
print*,' ' |
| 137 |
|
| 138 |
# endif |
| 139 |
|
| 140 |
c------------------------------------------------------------------------ |
| 141 |
c |
| 142 |
c HBOOK initialization |
| 143 |
c |
| 144 |
c------------------------------------------------------------------------ |
| 145 |
|
| 146 |
call HLIMIT(NWPAWC) |
| 147 |
|
| 148 |
|
| 149 |
c------------------------------------------------------------------------ |
| 150 |
c |
| 151 |
c reads input informations |
| 152 |
c |
| 153 |
c------------------------------------------------------------------------ |
| 154 |
call fdate(processing_date) |
| 155 |
write(*,101) |
| 156 |
$ processing_date |
| 157 |
101 format(/ |
| 158 |
$ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/ |
| 159 |
$ ,'* *',/ |
| 160 |
$ ,'* CALIBRATION *',/ |
| 161 |
$ ,'* *',/ |
| 162 |
$ ,'*** *** *** *** *** *** *** *** *** *** *** *** ***',/ |
| 163 |
$ ,a24,/ |
| 164 |
$ ) |
| 165 |
|
| 166 |
print*,'Number of files to be analysed:' |
| 167 |
read(*,*)nfile |
| 168 |
print*,nfile |
| 169 |
111 format(a) |
| 170 |
print*,'Input data directory:' |
| 171 |
read(*,111)data_dir |
| 172 |
print*,data_dir |
| 173 |
print*,'Output data directory:' |
| 174 |
read(*,111)out_dir |
| 175 |
print*,out_dir |
| 176 |
print*,'List of file identifiers: (DATE_NUM)' |
| 177 |
300 format(A10) |
| 178 |
do i=1,nfile |
| 179 |
read(*,300)input_string(i) |
| 180 |
print*,input_string(i) |
| 181 |
enddo |
| 182 |
minevent=1 |
| 183 |
print*,'Maximum number of events to be analized:' !20000 |
| 184 |
read(*,*) ntotev |
| 185 |
print*,ntotev |
| 186 |
print*,'DEBUG mode: (T, F)' |
| 187 |
11 format(l1) |
| 188 |
read(*,11)DEBUG |
| 189 |
print*,DEBUG |
| 190 |
print*,'---------------------------------------------------' |
| 191 |
|
| 192 |
|
| 193 |
****** INITIALIZATIONS ************************************* |
| 194 |
|
| 195 |
|
| 196 |
* 2) read z coordinates of the planes |
| 197 |
print*,' ' |
| 198 |
print*,'- read z coordinates of the planes' |
| 199 |
print*,' ' |
| 200 |
call mech_sensor !reads sensors centres coordinates |
| 201 |
do ip=1,nplanes |
| 202 |
fitz(ip)=z_mech_sensor(ip,1,1)*0.1 !cm |
| 203 |
* gets planes mechanical z positions |
| 204 |
* (in mm) and sets them in micrometers |
| 205 |
enddo |
| 206 |
|
| 207 |
* 4) read magnetic field map |
| 208 |
print*,' ' |
| 209 |
print*,'- read magnetic field map' |
| 210 |
print*,' ' |
| 211 |
call read_B |
| 212 |
|
| 213 |
* 5) read allignment parameters |
| 214 |
print*,' ' |
| 215 |
print*,'- read aligment parameters' |
| 216 |
print*,' ' |
| 217 |
call readalignparam |
| 218 |
|
| 219 |
************************************************************ |
| 220 |
************************************************************ |
| 221 |
************************************************************ |
| 222 |
************************************************************ |
| 223 |
|
| 224 |
print*,' ' |
| 225 |
print*,'---------------------------------------------------' |
| 226 |
print*,' ' |
| 227 |
print*,'OPENING OUTPUT-HISTO FILES:' |
| 228 |
|
| 229 |
1010 format('LADDER',i1) |
| 230 |
do ilad=1,3 |
| 231 |
|
| 232 |
write(rzdir,1010)ilad |
| 233 |
file_out= |
| 234 |
$ out_dir(1:LNBLNK(out_dir))//'trk-' |
| 235 |
$ //rzdir//'-histos.rz' |
| 236 |
print*,' -> '//file_out |
| 237 |
IQUEST(10)=65000 |
| 238 |
call HROPEN(lun_file_out+ilad,rzdir,file_out,'PQN',4096,istat) |
| 239 |
if(istat.ne.0) goto 16 |
| 240 |
call HCDIR('//'//rzdir,' ') |
| 241 |
|
| 242 |
call HCDIR('//PAWC','') |
| 243 |
call HMDIR(rzdir,'S') |
| 244 |
|
| 245 |
call book_eta_histos |
| 246 |
call book_charge_histos |
| 247 |
call book_check_histos |
| 248 |
call HCDIR('//PAWC','') |
| 249 |
|
| 250 |
enddo |
| 251 |
|
| 252 |
************************************************************ |
| 253 |
************************************************************ |
| 254 |
************************************************************ |
| 255 |
************************************************************ |
| 256 |
************************************************************ |
| 257 |
************************************************************ |
| 258 |
************************************************************ |
| 259 |
************************************************************ |
| 260 |
************************************************************ |
| 261 |
************************************************************ |
| 262 |
************************************************************ |
| 263 |
************************************************************ |
| 264 |
|
| 265 |
c------------------------------------------------------------------------ |
| 266 |
c |
| 267 |
c LEVEL 2 ntuple booking |
| 268 |
c |
| 269 |
c------------------------------------------------------------------------ |
| 270 |
|
| 271 |
|
| 272 |
503 format(a,'DW__level2-calib.rz') |
| 273 |
write(data_file_level2,503) |
| 274 |
$ out_dir(1:LNBLNK(out_dir)) |
| 275 |
c $ ,data_file(1:LNBLNK(data_file)) |
| 276 |
|
| 277 |
print*,'' |
| 278 |
print*,'------------------------------------' |
| 279 |
print*,' Creating LEVEL2 rz file' |
| 280 |
print*,' ',data_file_level2 |
| 281 |
print*,'------------------------------------' |
| 282 |
|
| 283 |
|
| 284 |
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
| 285 |
C largest RZ file: IQUEST(10) records x LREC words x 4 byte |
| 286 |
C with the following settings: 65000 x 4096 x 4 = 1G |
| 287 |
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - |
| 288 |
IQUEST(10)=65000 |
| 289 |
c !permette di ottenere ntuple funzionanti nonostante |
| 290 |
c ! il messaggio dei 64K di RZOUT...!??? |
| 291 |
call HROPEN(lun_data_level2, |
| 292 |
$ 'LEVEL2',data_file_level2,'QNP',4096,istat) !opens rz |
| 293 |
if(istat.ne.0) goto 19 |
| 294 |
call book_level2 |
| 295 |
|
| 296 |
print*,' ' |
| 297 |
print*,'---------------------------------------------------' |
| 298 |
print*,' * * * START ANALYSIS * * *' |
| 299 |
print*,'---------------------------------------------------' |
| 300 |
|
| 301 |
c------------------------------------------------------------------------ |
| 302 |
c |
| 303 |
c loop on files |
| 304 |
c |
| 305 |
c------------------------------------------------------------------------ |
| 306 |
|
| 307 |
goodtr = 0 |
| 308 |
tottr = 0 |
| 309 |
do ifile = 1,nfile !files loop |
| 310 |
|
| 311 |
print*,' ' |
| 312 |
data_file = input_string(ifile) |
| 313 |
print*,'--> FILE ',ifile,': '//data_file(1:LNBLNK(data_file)) |
| 314 |
|
| 315 |
c------------------------------------------------------------------------ |
| 316 |
c |
| 317 |
c fills bad variables with online DSP bad strips from calib histograms |
| 318 |
c for bad strip exclusion |
| 319 |
c |
| 320 |
c------------------------------------------------------------------------ |
| 321 |
|
| 322 |
# ifdef TEST2003 |
| 323 |
c------------------------------------------------------------------------ |
| 324 |
c |
| 325 |
c opens calib file |
| 326 |
c |
| 327 |
c------------------------------------------------------------------------ |
| 328 |
|
| 329 |
505 format(a,'output_',a,'_calib.rz') |
| 330 |
write(data_file_calib,505) |
| 331 |
$ data_dir(1:LNBLNK(data_dir)) |
| 332 |
$ ,data_file(1:LNBLNK(data_file)) |
| 333 |
|
| 334 |
print*,' ' |
| 335 |
print*,'OPENING CALIB FILE...', data_file_calib |
| 336 |
|
| 337 |
IQUEST(10)=65000 |
| 338 |
call HROPEN |
| 339 |
$ (lun_data_calib,'IN',data_file_calib,'QP',4096,istat) |
| 340 |
if(istat.ne.0) goto 17 |
| 341 |
call HCDIR('//IN',' ') |
| 342 |
|
| 343 |
call HRIN(0,9999,0) !puts histograms in memory |
| 344 |
|
| 345 |
print*,' ' |
| 346 |
|
| 347 |
print*,' ' |
| 348 |
print*,'READING PEDESTAL, SIGMA AND BADSTRIP HISTOGRAMS...' |
| 349 |
print*,' ' |
| 350 |
|
| 351 |
call fillpedsig |
| 352 |
|
| 353 |
call HREND('IN') |
| 354 |
close (lun_data_calib) |
| 355 |
# else |
| 356 |
|
| 357 |
print*,' ' |
| 358 |
print*,'OPENING CALIBRATION-LIST FILE:' |
| 359 |
501 format(a,'DW_',a,'_calib.txt') |
| 360 |
write(data_file_calib,501) |
| 361 |
$ data_dir(1:LNBLNK(data_dir)) |
| 362 |
$ ,data_file(1:LNBLNK(data_file)) |
| 363 |
print*,data_file_calib |
| 364 |
open(lun_calib_list, |
| 365 |
$ FILE=data_file_calib(1:LNBLNK(data_file_calib)) |
| 366 |
$ ,STATUS='UNKNOWN' |
| 367 |
$ ,IOSTAT=iostat |
| 368 |
$ ) |
| 369 |
|
| 370 |
113 format(i5,' ',a25) |
| 371 |
|
| 372 |
do i=1,ncalibmax |
| 373 |
|
| 374 |
read(lun_calib_list,113,IOSTAT=iostat) n_cal_list |
| 375 |
$ ,file_calib(i) |
| 376 |
c file_calib(i)=file_calib(i)(1:LNBLNK(file_calib(i))) |
| 377 |
if(iostat.ne.0)then |
| 378 |
ncal=i-1 |
| 379 |
goto 2000 |
| 380 |
endif |
| 381 |
print*,n_cal_list,' - ' |
| 382 |
$ ,file_calib(i)(1:LNBLNK(file_calib(i))) |
| 383 |
|
| 384 |
enddo |
| 385 |
|
| 386 |
2000 close(lun_calib_list) |
| 387 |
which_calib_last=0 |
| 388 |
|
| 389 |
# endif |
| 390 |
|
| 391 |
c------------------------------------------------------------------------ |
| 392 |
c |
| 393 |
c opens level1 file |
| 394 |
c |
| 395 |
c------------------------------------------------------------------------ |
| 396 |
|
| 397 |
# ifdef TEST2003 |
| 398 |
504 format(a,'output_',a,'_level1.rz') |
| 399 |
# else |
| 400 |
504 format(a,'DW_',a,'_level1.rz') |
| 401 |
# endif |
| 402 |
write(data_file_level1,504) |
| 403 |
$ data_dir(1:LNBLNK(data_dir)) |
| 404 |
$ ,data_file(1:LNBLNK(data_file)) |
| 405 |
print*,'' |
| 406 |
print*,'OPENING LEVEL1 FILE:' |
| 407 |
print*,data_file_level1 |
| 408 |
IQUEST(10)=65000 |
| 409 |
c !permette di ottenere ntuple funzionanti nonostante |
| 410 |
c ! il messaggio dei 64K di RZOUT...!??? |
| 411 |
call HROPEN(lun_data_level1, |
| 412 |
$ 'LEVEL1',data_file_level1,'QP',4096,istat) !opens rz |
| 413 |
if(istat.ne.0) goto 19 |
| 414 |
call HRIN(ntp_level1,9999,20) |
| 415 |
call access_level1 |
| 416 |
c call HPRNTU(ntp_level1+20) |
| 417 |
call HNOENT(ntp_level1+20,iemax0) |
| 418 |
print*,'( ',iemax0,' events)' |
| 419 |
|
| 420 |
|
| 421 |
************************************************************ |
| 422 |
************************************************************ |
| 423 |
************************************************************ |
| 424 |
* |
| 425 |
* start track analysis |
| 426 |
* |
| 427 |
************************************************************ |
| 428 |
************************************************************ |
| 429 |
************************************************************ |
| 430 |
nselect = 0 |
| 431 |
maxevent=minevent+ntotev-1 |
| 432 |
do iev = minevent,MIN(iemax0,maxevent) !loop on events |
| 433 |
|
| 434 |
call HCDIR('//LEVEL1',' ') |
| 435 |
call HGNT(ntp_level1+20,iev,ierr) !reads an event |
| 436 |
if(ierr.ne.0) goto 21 |
| 437 |
|
| 438 |
tottr = tottr +1 |
| 439 |
|
| 440 |
*------------------------------------------------------ |
| 441 |
* LEVEL2 N-TUPLE INITIALIZATIONS |
| 442 |
call init_level2 |
| 443 |
if(.not.good1)then |
| 444 |
* goto 8800 !fill nt-uple and go to next event |
| 445 |
goto 100 !go to next event |
| 446 |
endif |
| 447 |
*------------------------------------------------------ |
| 448 |
|
| 449 |
|
| 450 |
# ifndef TEST2003 |
| 451 |
|
| 452 |
if(which_calib.ne.which_calib_last.and. |
| 453 |
$which_calib.ne.0)then |
| 454 |
|
| 455 |
data_file_calib= |
| 456 |
$ data_dir(1:LNBLNK(data_dir))// |
| 457 |
$ file_calib(which_calib) |
| 458 |
$ (1:LNBLNK(file_calib(which_calib))) |
| 459 |
c print*,data_file_calib |
| 460 |
print*,'' |
| 461 |
print*, |
| 462 |
$ '@ event ',nev2 |
| 463 |
$ ,' (CPU pkt N.',pkt_num1,')' |
| 464 |
print*,'--> ',data_file_calib |
| 465 |
IQUEST(10)=65000 |
| 466 |
call HROPEN(lun_data_calib, |
| 467 |
$ 'CALIB',data_file_calib,'QP',4096,istat) !opens |
| 468 |
if(istat.ne.0) goto 19 |
| 469 |
call HRIN(0,9999,0) |
| 470 |
call fillpedsig |
| 471 |
do iview=1,nviews |
| 472 |
call HDELET(id_hi_bad+iview) |
| 473 |
call HDELET(id_hi_ped+iview) |
| 474 |
call HDELET(id_hi_sig+iview) |
| 475 |
enddo |
| 476 |
call HREND('CALIB') |
| 477 |
close(lun_data_calib) |
| 478 |
which_calib_last=which_calib |
| 479 |
|
| 480 |
elseif(which_calib.eq.0)then |
| 481 |
|
| 482 |
nocalib=nocalib+1 |
| 483 |
good2=.false. |
| 484 |
goto 8800 !fill nt-uple and go to next event |
| 485 |
|
| 486 |
endif |
| 487 |
|
| 488 |
# endif |
| 489 |
|
| 490 |
|
| 491 |
*------------------------------------------------------ |
| 492 |
* cut on maximum number of clusters |
| 493 |
*------------------------------------------------------ |
| 494 |
* if(nclstr1.gt.nclstrmax_level2)then |
| 495 |
if(nclstr1.gt.20)then |
| 496 |
* goto 8800 !fill nt-uple and go to next event |
| 497 |
goto 100 !go to next event |
| 498 |
endif |
| 499 |
|
| 500 |
do i=1,nclstr1 |
| 501 |
cl_used(i)=0 !init mask of clusters associated to a track |
| 502 |
enddo |
| 503 |
|
| 504 |
if(DEBUG)then |
| 505 |
print*,'----------------------------------' |
| 506 |
print*,iev,' ** ',nev2 |
| 507 |
endif |
| 508 |
|
| 509 |
*------------------------------------------------------------------------------- |
| 510 |
* STEP 1 |
| 511 |
*------------------------------------------------------------------------------- |
| 512 |
* X-Y cluster association |
| 513 |
* |
| 514 |
* Clusters are associated to form COUPLES |
| 515 |
* Clusters not associated in any couple are called SINGLETS |
| 516 |
* |
| 517 |
* Track identification (Hough transform) and fitting is first done on couples. |
| 518 |
* Hence singlets are possibly added to the track. |
| 519 |
* |
| 520 |
* Variables assigned by the routine "cl_to_couples" are those in the |
| 521 |
* common blocks: |
| 522 |
* - common/clusters/cl_good |
| 523 |
* - common/couples/clx,cly,ncp_plane,ncp_tot,cp_useds1,cp_useds2 |
| 524 |
* - common/singlets/ncls,cls,cl_single |
| 525 |
* |
| 526 |
* NB - The routine version used for calibration is different from that used |
| 527 |
* for event processing. Different cluster selection is applied: |
| 528 |
* - "ngoodstr" good strip around MAXS required (ngoodstr=2 => 5 good strips) |
| 529 |
* - no charge correlation required |
| 530 |
*------------------------------------------------------------------------------- |
| 531 |
*------------------------------------------------------------------------------- |
| 532 |
|
| 533 |
iflag=0 |
| 534 |
call cl_to_couples_nocharge(iflag) |
| 535 |
if(iflag.eq.1)then !bad event |
| 536 |
* goto 880 !fill ntp and go to next event |
| 537 |
goto 100 !go to next event |
| 538 |
endif |
| 539 |
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 540 |
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 541 |
* PRE-SELEZIONE DI TRACCE PULITE |
| 542 |
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 543 |
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
| 544 |
* at least a couple per plane |
| 545 |
do ip=1,nplanes |
| 546 |
if(ncp_plane(ip).eq.0)goto 100 !next event |
| 547 |
enddo |
| 548 |
c if(good2.eq.0)goto 880!fill ntp and go to next event |
| 549 |
|
| 550 |
*----------------------------------------------------- |
| 551 |
*----------------------------------------------------- |
| 552 |
* HOUGH TRASFORM |
| 553 |
*----------------------------------------------------- |
| 554 |
*----------------------------------------------------- |
| 555 |
|
| 556 |
|
| 557 |
*------------------------------------------------------------------------------- |
| 558 |
* STEP 2 |
| 559 |
*------------------------------------------------------------------------------- |
| 560 |
* |
| 561 |
* Association of couples to form |
| 562 |
* - DOUBLETS in YZ view |
| 563 |
* - TRIPLETS in XZ view |
| 564 |
* |
| 565 |
* Variables assigned by the routine "cp_to_doubtrip" are those in the |
| 566 |
* common blocks: |
| 567 |
* - common/hough_param/ |
| 568 |
* $ alfayz1, !Y0 |
| 569 |
* $ alfayz2, !tg theta-yz |
| 570 |
* $ alfaxz1, !X0 |
| 571 |
* $ alfaxz2, !tg theta-xz |
| 572 |
* $ alfaxz3 !1/r |
| 573 |
* - common/doublets/ndblt,cpyz1,cpyz2 |
| 574 |
* - common/triplets/ntrpt,cpxz1,cpxz2,cpxz3 |
| 575 |
*------------------------------------------------------------------------------- |
| 576 |
*------------------------------------------------------------------------------- |
| 577 |
|
| 578 |
iflag=0 |
| 579 |
call cp_to_doubtrip(iflag) |
| 580 |
if(iflag.eq.1)then !bad event |
| 581 |
* goto 880 !fill ntp and go to next event |
| 582 |
goto 100 !go to next event |
| 583 |
endif |
| 584 |
|
| 585 |
|
| 586 |
*------------------------------------------------------------------------------- |
| 587 |
* STEP 3 |
| 588 |
*------------------------------------------------------------------------------- |
| 589 |
* |
| 590 |
* Classification of doublets and triplets to form CLOUDS, |
| 591 |
* according to distance in parameter space. |
| 592 |
* |
| 593 |
* cloud = cluster of points (doublets/triplets) in parameter space |
| 594 |
* |
| 595 |
* |
| 596 |
* |
| 597 |
* Variables assigned by the routine "doub_to_YZcloud" are those in the |
| 598 |
* common blocks: |
| 599 |
* - common/clouds_yz/ |
| 600 |
* $ nclouds_yz |
| 601 |
* $ ,alfayz1_av,alfayz2_av |
| 602 |
* $ ,ptcloud_yz,db_cloud,cpcloud_yz |
| 603 |
* |
| 604 |
* Variables assigned by the routine "trip_to_XZcloud" are those in the |
| 605 |
* common blocks: |
| 606 |
* common/clouds_xz/ |
| 607 |
* $ nclouds_xz |
| 608 |
* $ ,alfaxz1_av,alfaxz2_av,alfaxz3_av |
| 609 |
* $ ,ptcloud_xz,tr_cloud,cpcloud_xz |
| 610 |
*------------------------------------------------------------------------------- |
| 611 |
*------------------------------------------------------------------------------- |
| 612 |
|
| 613 |
iflag=0 |
| 614 |
call doub_to_YZcloud(iflag) |
| 615 |
if(iflag.eq.1)then !bad event |
| 616 |
* goto 880 !fill ntp and go to next event |
| 617 |
goto 100 !go to next event |
| 618 |
endif |
| 619 |
iflag=0 |
| 620 |
call trip_to_XZcloud(iflag) |
| 621 |
if(iflag.eq.1)then !bad event |
| 622 |
* goto 880 !fill ntp and go to next event |
| 623 |
goto 100 !go to next event |
| 624 |
endif |
| 625 |
|
| 626 |
|
| 627 |
*------------------------------------------------------------------------------- |
| 628 |
* STEP 4 (ITERATED until any other physical track isn't found) |
| 629 |
*------------------------------------------------------------------------------- |
| 630 |
* |
| 631 |
* YZ and XZ clouds are combined in order to obtain the initial guess |
| 632 |
* of the candidate-track parameters. |
| 633 |
* A minimum number of matching couples between YZ and XZ clouds is required. |
| 634 |
* |
| 635 |
* A TRACK CANDIDATE is defined by |
| 636 |
* - the couples resulting from the INTERSECTION of the two clouds, and |
| 637 |
* - the associated track parameters (evaluated by performing a zero-order |
| 638 |
* track fitting) |
| 639 |
* |
| 640 |
* The NTRACKS candidate-track parameters are stored in common block: |
| 641 |
* |
| 642 |
* - common/track_candidates/NTRACKS,AL_STORE |
| 643 |
* $ ,XV_STORE,YV_STORE,ZV_STORE |
| 644 |
* $ ,XM_STORE,YM_STORE,ZM_STORE |
| 645 |
* $ ,RESX_STORE,RESY_STORE |
| 646 |
* $ ,AXV_STORE,AYV_STORE |
| 647 |
* $ ,XGOOD_STORE,YGOOD_STORE |
| 648 |
* $ ,CP_STORE,RCHI2_STORE |
| 649 |
* |
| 650 |
*------------------------------------------------------------------------------- |
| 651 |
*------------------------------------------------------------------------------- |
| 652 |
ntrk=0 !counter of identified physical tracks |
| 653 |
|
| 654 |
11111 continue !<<<<<<< come here when performing a new search |
| 655 |
|
| 656 |
iflag=0 |
| 657 |
call clouds_to_ctrack(iflag) |
| 658 |
if(iflag.eq.1)then !no candidate tracks found |
| 659 |
* goto 880 !fill ntp and go to next event |
| 660 |
goto 100 !go to next event |
| 661 |
endif |
| 662 |
|
| 663 |
|
| 664 |
*------------------------------------------------------ |
| 665 |
* cut on maximum number of track candidates |
| 666 |
* (to have clean events) |
| 667 |
* NB: must be at least 2 to do not discard tracks |
| 668 |
* contained in a sensor column on Y view |
| 669 |
* |
| 670 |
* cut replaced with cut on the # cluster out of the track |
| 671 |
*------------------------------------------------------ |
| 672 |
* if(ntracks.gt.2)then |
| 673 |
* goto 8800 !fill nt-uple and go to next event |
| 674 |
* goto 100 !go to next event |
| 675 |
* endif |
| 676 |
|
| 677 |
|
| 678 |
* if(DEBUG)print*,'*GOOD!*' |
| 679 |
|
| 680 |
FIMAGE=.false. !processing best track (not track image) |
| 681 |
ibest=0 !best track among candidates |
| 682 |
iimage=0 !track image |
| 683 |
* ----------------------------------------------------- |
| 684 |
* --------------- select the best track --------------- |
| 685 |
* selection cuts: |
| 686 |
* - NX = NY = 6 |
| 687 |
* ----------------------------------------------------- |
| 688 |
* |
| 689 |
rchi2best=100000. |
| 690 |
do i=1,ntracks |
| 691 |
|
| 692 |
npx=0 |
| 693 |
npy=0 |
| 694 |
do ip=1,nplanes |
| 695 |
npx=npx+XGOOD_STORE(ip,i) |
| 696 |
npy=npy+YGOOD_STORE(ip,i) |
| 697 |
enddo |
| 698 |
* print*, |
| 699 |
* $ 'Candidate ',i,' nx,ny ',npx,npy,' chi2 ',RCHI2_STORE(i) |
| 700 |
if(npx.eq.nplanes.and.npy.eq.nplanes)then |
| 701 |
if(RCHI2_STORE(i).lt.rchi2best.and. |
| 702 |
$ RCHI2_STORE(i).gt.0)then |
| 703 |
ibest=i |
| 704 |
rchi2best=RCHI2_STORE(i) |
| 705 |
endif |
| 706 |
endif |
| 707 |
enddo |
| 708 |
* if(ibest.eq.0)goto 880 !>> no good candidates |
| 709 |
if(ibest.eq.0)goto 100 !>> no good candidates |
| 710 |
|
| 711 |
* print*,'BEST TRACK ',ibest,' nx,ny',npx,npy |
| 712 |
|
| 713 |
* print*,'>>>> ',iev |
| 714 |
|
| 715 |
* ----------------------------------------------------- |
| 716 |
* ------------- store the best track only ------------- |
| 717 |
* ----------------------------------------------------- |
| 718 |
ntr=1 |
| 719 |
|
| 720 |
ntrk=ntr |
| 721 |
image(ntr)=0 |
| 722 |
good2=.true. |
| 723 |
chi2_nt(ntr) = RCHI2_STORE(ibest) |
| 724 |
c print*,chi2_nt(ntr) |
| 725 |
|
| 726 |
* **************************************************** |
| 727 |
* normalize the angles |
| 728 |
* **************************************************** |
| 729 |
pig=acos(-1.) |
| 730 |
phi = AL_STORE(4,ibest) |
| 731 |
sinth = AL_STORE(3,ibest) |
| 732 |
if(sinth.lt.0)then |
| 733 |
sinth = -sinth |
| 734 |
phi = phi + pig |
| 735 |
endif |
| 736 |
npig = aint(phi/(2*pig)) |
| 737 |
phi = phi - npig*2*pig |
| 738 |
if(phi.lt.0) |
| 739 |
$ phi = phi + 2*pig |
| 740 |
AL_STORE(4,ibest) = phi |
| 741 |
AL_STORE(3,ibest) = sinth |
| 742 |
* **************************************************** |
| 743 |
* ... then assign the other variables |
| 744 |
* **************************************************** |
| 745 |
do i=1,5 |
| 746 |
al_nt(i,ntr) = AL_STORE(i,ibest) |
| 747 |
c print*,i,'al_nt(i,ntr)',al_nt(i,ntr) |
| 748 |
do j=1,5 |
| 749 |
coval(i,j,ntr) = 0 ! TEMPORANEO!!!!!!!!!!!! |
| 750 |
enddo |
| 751 |
enddo |
| 752 |
do ip=1,nplanes ! loop on planes |
| 753 |
xgood_nt(ip,ntr) = XGOOD_STORE(ip,ibest) |
| 754 |
ygood_nt(ip,ntr) = YGOOD_STORE(ip,ibest) |
| 755 |
xm_nt(ip,ntr) = XM_STORE(ip,ibest) |
| 756 |
ym_nt(ip,ntr) = YM_STORE(ip,ibest) |
| 757 |
zm_nt(ip,ntr) = ZM_STORE(ip,ibest) |
| 758 |
RESX_nt(IP,ntr) = RESX_STORE(ip,ibest) |
| 759 |
RESY_nt(IP,ntr) = RESY_STORE(ip,ibest) |
| 760 |
xv_nt(ip,ntr) = XV_STORE(ip,ibest) |
| 761 |
yv_nt(ip,ntr) = YV_STORE(ip,ibest) |
| 762 |
zv_nt(ip,ntr) = ZV_STORE(ip,ibest) |
| 763 |
axv_nt(ip,ntr) = AXV_STORE(ip,ibest) |
| 764 |
ayv_nt(ip,ntr) = AYV_STORE(ip,ibest) |
| 765 |
|
| 766 |
c print*,ip,' xm_nt(ip,ntr)', xm_nt(ip,ntr) |
| 767 |
c print*,ip,' ym_nt(ip,ntr)', ym_nt(ip,ntr) |
| 768 |
c print*,ip,' zm_nt(ip,ntr)', zm_nt(ip,ntr) |
| 769 |
c print*,ip,' RESX_nt(ip,ntr)', RESX_nt(ip,ntr) |
| 770 |
c print*,ip,' RESY_nt(ip,ntr)', RESY_nt(ip,ntr) |
| 771 |
c print*,ip,' xv_nt(ip,ntr)',xv_nt(ip,ntr) |
| 772 |
c print*,ip,' yv_nt(ip,ntr)',yv_nt(ip,ntr) |
| 773 |
c print*,ip,' zv_nt(ip,ntr)',zv_nt(ip,ntr) |
| 774 |
c print*,ip,' axv_nt(ip,ntr)',axv_nt(ip,ntr) |
| 775 |
c print*,ip,' ayv_nt(ip,ntr)',ayv_nt(ip,ntr) |
| 776 |
|
| 777 |
id = CP_STORE(ip,ibest) |
| 778 |
icp = icp_cp(id) |
| 779 |
iclx = clx(nplanes-ip+1,icp) |
| 780 |
icly = cly(nplanes-ip+1,icp) |
| 781 |
which_cl_x(ip) = iclx |
| 782 |
which_cl_y(ip) = icly |
| 783 |
dedx_x(ip,ntr) = dedx(iclx) |
| 784 |
dedx_y(ip,ntr) = dedx(icly) |
| 785 |
enddo |
| 786 |
|
| 787 |
c call CalcBdL(100,xxxx,IFAIL) |
| 788 |
c if(ifps(xxxx).eq.1)BdL(ntr) = xxxx |
| 789 |
BdL(ntr) = 0 |
| 790 |
|
| 791 |
* --- then remove selected clusters (ibest+iimage) from clouds ---- |
| 792 |
call clean_XYclouds(ibest) |
| 793 |
if(iflag.eq.1)then !bad event |
| 794 |
* goto 880 !fill ntp and go to next event |
| 795 |
goto 100 !go to next event |
| 796 |
endif |
| 797 |
|
| 798 |
* >>>>>>>>>>>>>>>>>>> NT-UPLE filling <<<<<<<<<<<<<<<<<<<< |
| 799 |
* + + + + + + + + + + + + + + + + + |
| 800 |
* / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / \ / |
| 801 |
* + + + + + + + + + + + + + + + + + |
| 802 |
|
| 803 |
880 continue |
| 804 |
* ********************************************************** |
| 805 |
* stores info about clusters not associated with any track |
| 806 |
* ********************************************************** |
| 807 |
c call fill_level2_siglets |
| 808 |
|
| 809 |
if(DEBUG)then |
| 810 |
|
| 811 |
print*,'' |
| 812 |
print*,'DONE!' |
| 813 |
print*,'' |
| 814 |
print*,'* summary *' |
| 815 |
print*,'tracks ',ntrk |
| 816 |
print*,'cl used ',(cl_used(i),i=1,nclstr1) |
| 817 |
c$$$ print*,'cl unused (x-y)' |
| 818 |
c$$$ do ip=1,nplanes |
| 819 |
c$$$ print*,ip,' << ',nclsx(ip),nclsy(ip) |
| 820 |
c$$$ enddo |
| 821 |
print*,'' |
| 822 |
print*,'' |
| 823 |
endif |
| 824 |
|
| 825 |
|
| 826 |
8800 continue |
| 827 |
|
| 828 |
|
| 829 |
* ----------------------------------------------------- |
| 830 |
* +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 831 |
* calibrations: |
| 832 |
* +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 833 |
* ----------------------------------------------------- |
| 834 |
|
| 835 |
|
| 836 |
* ================================== |
| 837 |
* ////////////////////////////////// |
| 838 |
* Track selection |
| 839 |
* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ |
| 840 |
* ================================== |
| 841 |
* - cut on number of planes |
| 842 |
* (fatto a monte... non mi ricordo perche') |
| 843 |
c$$$ npx=0 |
| 844 |
c$$$ npy=0 |
| 845 |
c$$$ do ip=1,nplanes |
| 846 |
c$$$ npx=npx+XGOOD_STORE(ip,i) |
| 847 |
c$$$ npy=npy+YGOOD_STORE(ip,i) |
| 848 |
c$$$ enddo |
| 849 |
c$$$ if(.not.(npx.ge.5..and.npy.ge.5))goto 100 |
| 850 |
* - chi^2 selection |
| 851 |
rig = abs(1./al_nt(5,1)) |
| 852 |
c$$$ chicut = 5. !100. |
| 853 |
c$$$ chicut = chicut * ( 1.48 + 4.64/rig + 18.35/rig**2 ) |
| 854 |
c$$$ if(chi2_nt(1).gt.chicut)goto 100 !next event |
| 855 |
chicut = 20 |
| 856 |
chicut = chicut *( |
| 857 |
c $ 3.9633*exp(-0.54485E-01*rig) |
| 858 |
c $ +0.99231 |
| 859 |
c $ +10.835*exp(-0.67155*rig) ) |
| 860 |
$ 3.7032*exp(-0.41911E-01*rig) |
| 861 |
$ +0.85355 |
| 862 |
$ +10.905*exp(-0.61893*rig) ) |
| 863 |
if(chi2_nt(1).gt.chicut.or.chi2_nt(1).lt.0)goto 100 !next event |
| 864 |
* - cut on number of cluster out of the track (noise) |
| 865 |
* (no constraints on the last plane becouse of the |
| 866 |
* back-scattering from calorimeter) |
| 867 |
c$$$ do ip=1,5 !6 |
| 868 |
c$$$ nclstot=nclsx(ip)+nclsy(ip) |
| 869 |
c$$$ if(nclstot.ge.1)goto 100 !next event |
| 870 |
c$$$ enddo |
| 871 |
nclstot = 0 |
| 872 |
do isx=1,nclsx |
| 873 |
c if(planex(isx).ne.6)nclstot = nclstot +1 |
| 874 |
nclstot = nclstot +1 |
| 875 |
enddo |
| 876 |
do isy=1,nclsy |
| 877 |
c if(planey(isy).ne.6)nclstot = nclstot +1 |
| 878 |
nclstot = nclstot +1 |
| 879 |
enddo |
| 880 |
if(nclstot.ne.0)print*,'iev ',iev,' ncls ',nclstot |
| 881 |
c if(nclstot.ge.1)goto 100 |
| 882 |
* ================================== |
| 883 |
* ////////////////////////////////// |
| 884 |
* Track selection (end) |
| 885 |
* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ |
| 886 |
* ================================== |
| 887 |
|
| 888 |
call HCDIR('//LEVEL2',' ') |
| 889 |
|
| 890 |
goodtr = goodtr +1 |
| 891 |
nselect = nselect +1 |
| 892 |
c print*,goodtr,' ****',iev |
| 893 |
call HFNT(ntp_level2) !fill LEVEL2 nt-uple |
| 894 |
|
| 895 |
|
| 896 |
call HCDIR('//PAWC','') |
| 897 |
|
| 898 |
* ----------------------------------------------------- |
| 899 |
* +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 900 |
* calibrations (sum over events) |
| 901 |
* +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 902 |
* ----------------------------------------------------- |
| 903 |
|
| 904 |
if(goodtr.lt.ncorrmax)then |
| 905 |
do ip=1,nplanes !1-->6 |
| 906 |
|
| 907 |
iviewx = (nplanes - ip)*2+2 !12-->2 |
| 908 |
ilx=nld(maxs(which_cl_x(ip)),iviewx) |
| 909 |
iviewy = (nplanes - ip)*2+1 !11-->1 |
| 910 |
ily=nld(maxs(which_cl_y(ip)),iviewy) |
| 911 |
if(ilx.ne.ily)print*,'LADDERS DO NOT MATCH!!!' |
| 912 |
|
| 913 |
chargexy(ip,ilx,goodtr,1)=dedx(which_cl_x(ip)) |
| 914 |
chargexy(ip,ily,goodtr,2)=dedx(which_cl_y(ip)) |
| 915 |
|
| 916 |
enddo |
| 917 |
endif |
| 918 |
* compute average amgle in each bin |
| 919 |
do i=1,nviews |
| 920 |
do ia=1,nangbin |
| 921 |
if(navangbin(ia,i).ne.0) |
| 922 |
$ avangbin(ia,i)= |
| 923 |
$ avangbin(ia,i)/navangbin(ia,i) |
| 924 |
enddo |
| 925 |
enddo |
| 926 |
|
| 927 |
* apply a cut on the measured rigidity: |
| 928 |
if(rig.gt.rigmax)then |
| 929 |
call fill_eta_histos |
| 930 |
call fill_check_histos |
| 931 |
endif |
| 932 |
|
| 933 |
* ----------------------------------------------------- |
| 934 |
* +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 935 |
* calibrations (end) |
| 936 |
* +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 937 |
* ----------------------------------------------------- |
| 938 |
100 continue |
| 939 |
|
| 940 |
enddo !end loop on events |
| 941 |
|
| 942 |
|
| 943 |
print*,'' |
| 944 |
print *,'- Stored ',nselect,' entries in LEVEL2 nt-uple' |
| 945 |
call HCDIR('//LEVEL1',' ') |
| 946 |
call HREND('level1') |
| 947 |
close(lun_data_level1) |
| 948 |
|
| 949 |
call HCDIR('//PAWC',' ') |
| 950 |
call HDELET(ntp_level1+20) |
| 951 |
|
| 952 |
enddo !end loop on files |
| 953 |
|
| 954 |
print*,'' |
| 955 |
call HCDIR('//LEVEL2',' ') |
| 956 |
call RZPURG(-1) |
| 957 |
call HROUT(ntp_level2,ICYCLE,'T') |
| 958 |
call HREND('level2') |
| 959 |
close(lun_data_level2) |
| 960 |
|
| 961 |
call HCDIR('//PAWC',' ') |
| 962 |
call HDELET(ntp_level2) |
| 963 |
|
| 964 |
|
| 965 |
* =================================================== |
| 966 |
* in the following the charge correlation parameters |
| 967 |
* and the eta correction functions are evaluated |
| 968 |
* and stored in txt files. |
| 969 |
* However I chosed to do it "off-line", by means of |
| 970 |
* a kumac, so that I can group more files together. |
| 971 |
* =================================================== |
| 972 |
|
| 973 |
* ----------------------------------------------------- |
| 974 |
* +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 975 |
* calibrations |
| 976 |
* +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 977 |
* ----------------------------------------------------- |
| 978 |
|
| 979 |
do il=1,3 |
| 980 |
write(rzdir,1010)il |
| 981 |
************************************************** |
| 982 |
************************************************** |
| 983 |
* CHARGE CORRELATION CALIBRATION |
| 984 |
************************************************** |
| 985 |
************************************************** |
| 986 |
if(DONTDO)goto 8881 |
| 987 |
call charge_calibration(il) |
| 988 |
107 format(/ |
| 989 |
$ ,/ |
| 990 |
$ ,/ |
| 991 |
$ ' >>> CHARGE-CORRELATION CALIBRATION <<<',/ |
| 992 |
$ ,/ |
| 993 |
$ 'Save charge correlation parameters --->',/ |
| 994 |
$ ,a,/ |
| 995 |
$ ,/ |
| 996 |
$ ,'Parameterization ------> Sy = K * Sx + C',/ |
| 997 |
$ ,'Plane K EK C EC',/ |
| 998 |
$ ,'-----------------------------------------') |
| 999 |
c write(lun_file_log,107)file_out(1:LNBLNK(file_out)) |
| 1000 |
|
| 1001 |
file_out= |
| 1002 |
$ out_dir(1:LNBLNK(out_dir))//'trk-' |
| 1003 |
$ //rzdir//'-charge.dat' |
| 1004 |
write(*,107)file_out(1:LNBLNK(file_out)) |
| 1005 |
open(55, |
| 1006 |
$ FILE=file_out(1:LNBLNK(file_out)), |
| 1007 |
c $ STATUS='NEW', |
| 1008 |
$ FORM='FORMATTED') |
| 1009 |
|
| 1010 |
108 format(i5,4f9.5) |
| 1011 |
do ip=1,nplanes |
| 1012 |
c write(lun_file_log,108) !log file |
| 1013 |
c $ ip,syksx(ip),esyksx(ip),sycsx(ip),esycsx(ip) |
| 1014 |
write(*,108) !screen |
| 1015 |
$ ip,syksx(ip),esyksx(ip),sycsx(ip),esycsx(ip) |
| 1016 |
write(55,*) !data |
| 1017 |
$ ip,syksx(ip),esyksx(ip),sycsx(ip),esycsx(ip) |
| 1018 |
enddo |
| 1019 |
close(55) |
| 1020 |
|
| 1021 |
8881 call fill_charge_histos(il) |
| 1022 |
************************************************** |
| 1023 |
************************************************** |
| 1024 |
* ETA CALIBRATION |
| 1025 |
************************************************** |
| 1026 |
************************************************** |
| 1027 |
file_out= |
| 1028 |
$ out_dir(1:LNBLNK(out_dir))//'trk-' |
| 1029 |
$ //rzdir |
| 1030 |
|
| 1031 |
109 format('Creating files:',/ |
| 1032 |
$ ,' ',a,'-eta2-binning.dat',/ |
| 1033 |
$ ) |
| 1034 |
write(*,109)file_out(1:LNBLNK(file_out)) |
| 1035 |
open(lun_file_eta, |
| 1036 |
$ FILE=file_out(1:LNBLNK(file_out))//'-eta2-binning.dat' |
| 1037 |
$ ,FORM='FORMATTED' |
| 1038 |
$ ) |
| 1039 |
do ia=1,nangbin |
| 1040 |
write(lun_file_eta,*)ia,angL(ia),angR(ia) |
| 1041 |
enddo |
| 1042 |
close(lun_file_eta) |
| 1043 |
|
| 1044 |
|
| 1045 |
if(DONTDO)goto 8882 |
| 1046 |
do ia=1,nang |
| 1047 |
|
| 1048 |
201 format('-eta2-bin',i1,'.dat') |
| 1049 |
202 format('-eta2-bin',i2,'.dat') |
| 1050 |
if(ia.lt.10)write(fname_param,201)ia |
| 1051 |
if(ia.ge.10)write(fname_param,202)ia |
| 1052 |
|
| 1053 |
open(lun_file_eta, |
| 1054 |
$ FILE=file_out(1:LNBLNK(file_out)) |
| 1055 |
$ //fname_param(1:LNBLNK(fname_param)) |
| 1056 |
$ ,FORM='FORMATTED' |
| 1057 |
$ ) |
| 1058 |
|
| 1059 |
call eta_calibration(il,ia) |
| 1060 |
|
| 1061 |
* ********************** write ************************** |
| 1062 |
203 format(13f10.4) |
| 1063 |
do ie=1,nhbin2+1 |
| 1064 |
write(lun_file_eta,203) |
| 1065 |
$ eta2abs(ie),(eta2corr(ie,i),i=1,nviews) |
| 1066 |
enddo |
| 1067 |
close(lun_file_eta) |
| 1068 |
* ********************** write ************************** |
| 1069 |
|
| 1070 |
enddo !end loop on angular bins |
| 1071 |
8882 continue |
| 1072 |
|
| 1073 |
enddo !end loop on LADDERs |
| 1074 |
|
| 1075 |
* ----------------------------------------------------- |
| 1076 |
* +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 1077 |
* calibrations |
| 1078 |
* +++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 1079 |
* ----------------------------------------------------- |
| 1080 |
8888 continue |
| 1081 |
|
| 1082 |
|
| 1083 |
c------------------------------------------------------------------------ |
| 1084 |
c |
| 1085 |
c no error exit |
| 1086 |
c |
| 1087 |
c------------------------------------------------------------------------ |
| 1088 |
|
| 1089 |
goto 9000 !happy ending |
| 1090 |
|
| 1091 |
|
| 1092 |
c------------------------------------------------------------------------ |
| 1093 |
c |
| 1094 |
c data file opening error |
| 1095 |
c |
| 1096 |
c------------------------------------------------------------------------ |
| 1097 |
19 continue |
| 1098 |
|
| 1099 |
print*,' ' |
| 1100 |
print*,'ERROR OPENING DATA FILE: ',data_file |
| 1101 |
print*,' ' |
| 1102 |
print*,' ' |
| 1103 |
|
| 1104 |
goto 9000 !the end |
| 1105 |
|
| 1106 |
c------------------------------------------------------------------------ |
| 1107 |
c |
| 1108 |
c level1 ntuple event reading error |
| 1109 |
c |
| 1110 |
c------------------------------------------------------------------------ |
| 1111 |
|
| 1112 |
21 continue |
| 1113 |
|
| 1114 |
print*,' ' |
| 1115 |
print*,'ERROR WHILE READING LEVEL1 NTUPLE, AT EVENT |
| 1116 |
$ : ',iev |
| 1117 |
print*,' ' |
| 1118 |
print*,' ' |
| 1119 |
|
| 1120 |
goto 9000 !the end |
| 1121 |
|
| 1122 |
c------------------------------------------------------------------------ |
| 1123 |
c |
| 1124 |
c calib file opening error |
| 1125 |
c |
| 1126 |
c------------------------------------------------------------------------ |
| 1127 |
|
| 1128 |
17 continue |
| 1129 |
|
| 1130 |
print*,' ' |
| 1131 |
print*,'preanalysis: ERROR OPENING INPUT FILE: ',data_file_calib |
| 1132 |
print*,' ' |
| 1133 |
print*,' ' |
| 1134 |
|
| 1135 |
goto 9000 !the end |
| 1136 |
|
| 1137 |
16 continue |
| 1138 |
|
| 1139 |
print*,' ' |
| 1140 |
print*,'preanalysis: ERROR OPENING OUTPUT FILE: ',file_out |
| 1141 |
print*,' ' |
| 1142 |
print*,' ' |
| 1143 |
|
| 1144 |
goto 9000 !the end |
| 1145 |
|
| 1146 |
c------------------------------------------------------------------------ |
| 1147 |
c |
| 1148 |
c closes files and exits |
| 1149 |
c |
| 1150 |
c------------------------------------------------------------------------ |
| 1151 |
|
| 1152 |
9000 continue |
| 1153 |
|
| 1154 |
c$$$ if(DEBUG)call HPRNTU(ntp_level2+1) |
| 1155 |
c$$$ call HPRNTU(ntp_level2) |
| 1156 |
c$$$ print*,'' |
| 1157 |
c$$$ call HCDIR('//LEVEL2',' ') |
| 1158 |
c$$$ call HROUT(ntp_level2,ICYCLE,'T') |
| 1159 |
c$$$ print *,'- Stored LEVEL2 nt-uple (',nev2,' entries )' |
| 1160 |
c$$$ if(DEBUG)call HROUT(ntp_level2+1,ICYCLE,'T') |
| 1161 |
c$$$ if(DEBUG)print *,'- Stored DEBUG nt-uple ' |
| 1162 |
c$$$ call HREND('level2') |
| 1163 |
c$$$ close(lun_data_level2) |
| 1164 |
|
| 1165 |
c$$$ call HCDIR('//LEVEL1',' ') |
| 1166 |
c$$$ call HREND('level1') |
| 1167 |
c$$$ close(lun_data_level1) |
| 1168 |
|
| 1169 |
write(*,105) |
| 1170 |
$ tottr |
| 1171 |
$ ,goodtr |
| 1172 |
105 format(/ |
| 1173 |
$ ,/ |
| 1174 |
$ ,'|------------------------------------------------|',/ |
| 1175 |
$ ,'| Number of processed events |',/ |
| 1176 |
$ ,'| Total ',i8,' |'/ |
| 1177 |
$ ,'| Good tracks ',i8,' |'/ |
| 1178 |
$ ,'|------------------------------------------------|',/ |
| 1179 |
$ ) |
| 1180 |
|
| 1181 |
|
| 1182 |
|
| 1183 |
do ilad=1,3 |
| 1184 |
write(rzdir,1010)ilad |
| 1185 |
|
| 1186 |
call HCDIR('//PAWC/'//rzdir,' ') |
| 1187 |
call HCDIR('//'//RZDIR,' ') |
| 1188 |
|
| 1189 |
call HROUT(id_hi_chi,ICYCLE,'T') |
| 1190 |
call HROUT(id_hi_chi+100,ICYCLE,'T') |
| 1191 |
call HROUT(id_hi_chi+200,ICYCLE,'T') |
| 1192 |
call HROUT(id_hi_rig,ICYCLE,'T') |
| 1193 |
call HROUT(id_hi_charge,ICYCLE,'T') |
| 1194 |
do i=1,nviews |
| 1195 |
call HROUT(id_hi_charge+i*100,ICYCLE,'T') |
| 1196 |
call HROUT(id_hi_charge+i*100+1,ICYCLE,'T') |
| 1197 |
call HROUT(id_hi_charge+i*100+2,ICYCLE,'T') |
| 1198 |
call HROUT(80000+id_hi_charge+i*100+2,ICYCLE,'T') |
| 1199 |
call HROUT(id_hi_mult+i*100,ICYCLE,'T') |
| 1200 |
call HROUT(80000+id_hi_mult+i*100,ICYCLE,'T') |
| 1201 |
do ia=1,nang |
| 1202 |
call HROUT(id_hi_eta2+i*100+ia,ICYCLE,'T') |
| 1203 |
call HROUT(id_hi_eta3+i*100+ia,ICYCLE,'T') |
| 1204 |
call HROUT(id_hi_eta4+i*100+ia,ICYCLE,'T') |
| 1205 |
c call HROUT(id_hi_charge+i*100+ia,ICYCLE,'T') |
| 1206 |
call HROUT(id_hi_mult+i*100+ia,ICYCLE,'T') |
| 1207 |
call HROUT(id_hi_mult+i*100+ia+50000,ICYCLE,'T') |
| 1208 |
enddo |
| 1209 |
c call HROUT(id_hi_chi+100+ia,ICYCLE,'T') |
| 1210 |
c call HROUT(id_hi_chi+200+ia,ICYCLE,'T') |
| 1211 |
enddo |
| 1212 |
do ip=1,nplanes |
| 1213 |
call HROUT(id_hi_angle+100+ip,ICYCLE,'T') |
| 1214 |
call HROUT(id_hi_angle+200+ip,ICYCLE,'T') |
| 1215 |
call HROUT(id_hi_angle_stat+100+ip,ICYCLE,'T') |
| 1216 |
call HROUT(id_hi_angle_stat+200+ip,ICYCLE,'T') |
| 1217 |
call HROUT(id_hi_chargeco+ip*100+4,ICYCLE,'T') |
| 1218 |
call HROUT(id_hi_chargeco+ip*100,ICYCLE,'T') |
| 1219 |
do is=1,2 |
| 1220 |
call HROUT(id_hi_beam+ip*100+is+200000,ICYCLE,'T') |
| 1221 |
call HROUT(id_hi_beam+ip*100+is+100000,ICYCLE,'T') |
| 1222 |
enddo |
| 1223 |
call HROUT(id_hi_bx+ip,ICYCLE,'T') |
| 1224 |
call HROUT(id_hi_by+ip,ICYCLE,'T') |
| 1225 |
call HROUT(id_hi_bz+ip,ICYCLE,'T') |
| 1226 |
enddo |
| 1227 |
|
| 1228 |
call HREND(RZDIR) |
| 1229 |
close(lun_file_out+ilad) |
| 1230 |
|
| 1231 |
enddo |
| 1232 |
|
| 1233 |
|
| 1234 |
|
| 1235 |
|
| 1236 |
|
| 1237 |
|
| 1238 |
stop |
| 1239 |
end |
| 1240 |
|
| 1241 |
|
| 1242 |
# include "momanhough-subroutines.F" |
| 1243 |
|
| 1244 |
|
| 1245 |
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** |
| 1246 |
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* |
| 1247 |
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** |
| 1248 |
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* |
| 1249 |
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** |
| 1250 |
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* |
| 1251 |
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** |
| 1252 |
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* |
| 1253 |
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** |
| 1254 |
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* |
| 1255 |
****** * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ****** |
| 1256 |
******* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ******* |
| 1257 |
|
| 1258 |
*************************************************************** |
| 1259 |
* |
| 1260 |
* |
| 1261 |
* |
| 1262 |
* |
| 1263 |
* |
| 1264 |
*************************************************************** |
| 1265 |
subroutine book_eta_histos |
| 1266 |
|
| 1267 |
include '../common/commontracker.f' |
| 1268 |
include '../common/calib.f' |
| 1269 |
include '../common/common_calibration.f' |
| 1270 |
|
| 1271 |
character*64 title !histos title |
| 1272 |
|
| 1273 |
502 format('Eta-2 (view ',i2,', angle bin ',i2,' )') |
| 1274 |
503 format('Eta-3 (view ',i2,', angle bin ',i2,' )') |
| 1275 |
504 format('Eta-4 (view ',i2,', angle bin ',i2,' )') |
| 1276 |
|
| 1277 |
do i=1,nviews |
| 1278 |
do ia=1,nangbin |
| 1279 |
write(title,502)i,ia |
| 1280 |
call HBOOK1(id_hi_eta2+i*100+ia |
| 1281 |
$ ,title |
| 1282 |
c $ ,nhbin2,0.,1.,0.) |
| 1283 |
$ ,nhbin2,-0.5,0.5,0.) |
| 1284 |
write(title,503)i,ia |
| 1285 |
call HBOOK1(id_hi_eta3+i*100+ia |
| 1286 |
$ ,title |
| 1287 |
c $ ,nhbin3,0.,2.,0.) |
| 1288 |
$ ,nhbin2,-0.5,0.5,0.) |
| 1289 |
write(title,504)i,ia |
| 1290 |
call HBOOK1(id_hi_eta4+i*100+ia |
| 1291 |
$ ,title |
| 1292 |
c $ ,nhbin4,0.,3.,0.) |
| 1293 |
c $ ,nhbin3,-1.,1.,0.) |
| 1294 |
$ ,nhbin3,-0.5,.5,0.) |
| 1295 |
enddo |
| 1296 |
enddo |
| 1297 |
|
| 1298 |
|
| 1299 |
return |
| 1300 |
end |
| 1301 |
C............................................................ |
| 1302 |
subroutine fill_eta_histos |
| 1303 |
|
| 1304 |
include '../common/commontracker.f' |
| 1305 |
include '../common/calib.f' |
| 1306 |
include '../common/common_calibration.f' |
| 1307 |
include '../common/level1.f' |
| 1308 |
include '../common/level2.f' |
| 1309 |
|
| 1310 |
character*10 rzdir |
| 1311 |
|
| 1312 |
* flag to tag low-multilpicity cluster |
| 1313 |
* (according to Samu definition) |
| 1314 |
logical LMSX,LMSY,DRAY |
| 1315 |
|
| 1316 |
|
| 1317 |
10 format('LADDER',i1) |
| 1318 |
* start loop on planes to fill histos |
| 1319 |
do ip=1,nplanes !1-->6 |
| 1320 |
|
| 1321 |
*--------------------------------------------------------- |
| 1322 |
* main x-y cluster parameters |
| 1323 |
*--------------------------------------------------------- |
| 1324 |
xxx = xm_nt(ip,1) |
| 1325 |
yyy = ym_nt(ip,1) |
| 1326 |
call whichsensor(ip,xxx,yyy,iladder,isensor) |
| 1327 |
if(iladder.eq.0.or.isensor.eq.0)then |
| 1328 |
goto 10102 |
| 1329 |
endif |
| 1330 |
|
| 1331 |
multx = mult(which_cl_x(ip)) |
| 1332 |
multy = mult(which_cl_y(ip)) |
| 1333 |
thetax = axv_nt(ip,1) |
| 1334 |
thetay = ayv_nt(ip,1) |
| 1335 |
LMSX=.false. |
| 1336 |
LMSY=.false. |
| 1337 |
* =============================================== |
| 1338 |
* multiplicity selection of SAMUELE |
| 1339 |
c$$$ if( |
| 1340 |
c$$$ $ (thetax.le.15..and.multx.le.2).or. |
| 1341 |
c$$$ $ (thetax.gt.15..and.multx.le.3).or. |
| 1342 |
c$$$ $ .false.)LMSX=.true. |
| 1343 |
c$$$ if( |
| 1344 |
c$$$ $ (multx.le.2).or. |
| 1345 |
c$$$ $ .false.)LMSY=.true. |
| 1346 |
* =============================================== |
| 1347 |
* looser cut ... after LANDI warning |
| 1348 |
* (I also require that both clx anc cly are low-multiplicity) |
| 1349 |
if( |
| 1350 |
$ (abs(thetax).le.10..and.multx.le.3).or. |
| 1351 |
$ ((abs(thetax).gt.10..and.abs(thetax).le.15.).and.multx.le.4).or. |
| 1352 |
$ (abs(thetax).gt.15..and.multx.le.5).or. |
| 1353 |
$ .false.)LMSX=.true. |
| 1354 |
if( |
| 1355 |
$ (multy.le.3).or. |
| 1356 |
$ .false.)LMSY=.true. |
| 1357 |
|
| 1358 |
|
| 1359 |
* charge scaling factor (to vertical trajectory) |
| 1360 |
pi=acos(-1.) |
| 1361 |
fact= |
| 1362 |
$ sqrt( |
| 1363 |
$ tan(pi*thetax/180.)**2 + |
| 1364 |
$ tan(pi*thetay/180.)**2 + |
| 1365 |
$ 1. |
| 1366 |
$ ) |
| 1367 |
dedxscalx = dedx(which_cl_x(ip)) / fact |
| 1368 |
dedxscaly = dedx(which_cl_y(ip)) / fact |
| 1369 |
dedxscal = (dedxscalx+dedxscaly)/2 |
| 1370 |
|
| 1371 |
DRAY=.false. |
| 1372 |
if(dedxscal.gt.dedxmax)DRAY=.true. |
| 1373 |
|
| 1374 |
*--------------------------------------------------------- |
| 1375 |
|
| 1376 |
|
| 1377 |
write(rzdir,10)iladder |
| 1378 |
call HCDIR('//'//rzdir,'') |
| 1379 |
call HCDIR('//PAWC/'//rzdir,'') |
| 1380 |
|
| 1381 |
*--------------------------------------------------------- |
| 1382 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 1383 |
* X VIEWS |
| 1384 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 1385 |
*--------------------------------------------------------- |
| 1386 |
|
| 1387 |
if( thetax.lt.angL(1).or. |
| 1388 |
$ thetax.ge.angR(nangbin).or. |
| 1389 |
$ .not.LMSX.or. |
| 1390 |
$ .not.LMSY.or. |
| 1391 |
$ DRAY.or. |
| 1392 |
$ .false.)goto 10101 ! skip |
| 1393 |
|
| 1394 |
iview = (nplanes - ip)*2+2 !12-->2 |
| 1395 |
e2=cog(2,which_cl_x(ip)) |
| 1396 |
e3=cog(3,which_cl_x(ip)) |
| 1397 |
e4=cog(4,which_cl_x(ip)) |
| 1398 |
do ia=1,nangbin |
| 1399 |
if (thetax.lt.angR(ia))goto 111 |
| 1400 |
enddo |
| 1401 |
111 ibin_ang = ia !angular bin |
| 1402 |
|
| 1403 |
avangbin(ibin_ang,iview)=avangbin(ibin_ang,iview)+thetax |
| 1404 |
navangbin(ibin_ang,iview)=navangbin(ibin_ang,iview)+1 |
| 1405 |
|
| 1406 |
*!!!!!!!!NB non sempre MAXS coincide con la strip massima !!!! |
| 1407 |
if(e2.ge.0.5)e2=e2-1. |
| 1408 |
if(e2.lt.-0.5)e2=e2+1. |
| 1409 |
if(e3.ge.0.5)e3=e3-1. |
| 1410 |
if(e3.lt.-0.5)e3=e3+1. |
| 1411 |
c if(e4.ge.1.)e4=e4-2. |
| 1412 |
c if(e4.lt.-1.)e4=e4+2. |
| 1413 |
if(e4.ge.0.5)e4=e4-1. |
| 1414 |
if(e4.lt.-0.5)e4=e4+1. |
| 1415 |
*!!!!!!!!NB |
| 1416 |
|
| 1417 |
if(e2.ne.0.) |
| 1418 |
$ call hfill(id_hi_eta2+iview*100+ibin_ang, |
| 1419 |
$ e2,0.,1.) |
| 1420 |
if(e3.ne.0.) |
| 1421 |
$ call hfill(id_hi_eta3+iview*100+ibin_ang, |
| 1422 |
$ e3,0.,1.) |
| 1423 |
if(e4.ne.0.) |
| 1424 |
$ call hfill(id_hi_eta4+iview*100+ibin_ang, |
| 1425 |
$ e4,0.,1.) |
| 1426 |
|
| 1427 |
10101 continue |
| 1428 |
|
| 1429 |
*--------------------------------------------------------- |
| 1430 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 1431 |
* Y VIEWS |
| 1432 |
*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ |
| 1433 |
*--------------------------------------------------------- |
| 1434 |
if( thetay.lt.angL(1).or. |
| 1435 |
$ thetay.ge.angR(nangbin).or. |
| 1436 |
$ .not.LMSY.or. |
| 1437 |
$ .not.LMSX.or. |
| 1438 |
$ DRAY.or. |
| 1439 |
$ .false.)goto 10102 !skip |
| 1440 |
|
| 1441 |
iview = (nplanes - ip)*2+1 !11-->1 |
| 1442 |
e2=cog(2,which_cl_y(ip)) |
| 1443 |
e3=cog(3,which_cl_y(ip)) |
| 1444 |
e4=cog(4,which_cl_y(ip)) |
| 1445 |
do ia=1,nangbin |
| 1446 |
c if (thetay.lt.(ang(ia)+angstep/2.))goto 112 |
| 1447 |
c if (thetay.lt.angbin(ia))goto 112 |
| 1448 |
if (thetay.lt.angR(ia))goto 112 |
| 1449 |
enddo |
| 1450 |
112 ibin_ang = ia |
| 1451 |
|
| 1452 |
avangbin(ibin_ang,iview)=avangbin(ibin_ang,iview)+thetay |
| 1453 |
navangbin(ibin_ang,iview)=navangbin(ibin_ang,iview)+1 |
| 1454 |
c print*,'****',avangbin(ibin_ang,iview),thetay |
| 1455 |
|
| 1456 |
|
| 1457 |
*!!!!!!!!NB non sempre MAXS coincide con la strip massima !!!! |
| 1458 |
if(e2.gt.0.5)e2=e2-1. |
| 1459 |
if(e2.lt.-0.5)e2=e2+1. |
| 1460 |
if(e3.ge.0.5)e3=e3-1. |
| 1461 |
if(e3.lt.-0.5)e3=e3+1. |
| 1462 |
c if(e4.ge.1.)e4=e4-2. |
| 1463 |
c if(e4.lt.-1.)e4=e4+2. |
| 1464 |
if(e4.ge.0.5)e4=e4-1. |
| 1465 |
if(e4.lt.-0.5)e4=e4+1. |
| 1466 |
*!!!!!!!!NB |
| 1467 |
|
| 1468 |
if(e2.ne.0.) |
| 1469 |
$ call hfill(id_hi_eta2+iview*100+ibin_ang, |
| 1470 |
$ e2,0.,1.) |
| 1471 |
if(e3.ne.0.) |
| 1472 |
$ call hfill(id_hi_eta3+iview*100+ibin_ang, |
| 1473 |
$ e3,0.,1.) |
| 1474 |
if(e4.ne.0.) |
| 1475 |
$ call hfill(id_hi_eta4+iview*100+ibin_ang, |
| 1476 |
$ e4,0.,1.) |
| 1477 |
|
| 1478 |
10102 continue |
| 1479 |
enddo |
| 1480 |
call HCDIR('//PAWC','') |
| 1481 |
return |
| 1482 |
|
| 1483 |
end |
| 1484 |
|
| 1485 |
C............................................................ |
| 1486 |
subroutine book_charge_histos |
| 1487 |
|
| 1488 |
include '../common/commontracker.f' |
| 1489 |
include '../common/calib.f' |
| 1490 |
include '../common/common_calibration.f' |
| 1491 |
|
| 1492 |
character*64 title !histos titleinfo -f g77 M LEX |
| 1493 |
|
| 1494 |
c 500 format('Charge correlation (ladder ',i1,' plane ',i2,')') |
| 1495 |
500 format('Charge correlation ( plane ',i2,')') |
| 1496 |
501 format('Charge residuals ( plane',i2,')') |
| 1497 |
|
| 1498 |
do ip=1,nplanes |
| 1499 |
write(title,500)ip |
| 1500 |
call HBOOK2(id_hi_chargeco+ip*100+4 |
| 1501 |
$ ,title |
| 1502 |
$ ,200,0.,1000. |
| 1503 |
$ ,200,0.,1000.,0.) |
| 1504 |
write(title,501)ip |
| 1505 |
call HBOOK1(id_hi_chargeco+ip*100 |
| 1506 |
$ ,title |
| 1507 |
$ ,200,-100.,100.,0.) |
| 1508 |
enddo |
| 1509 |
|
| 1510 |
|
| 1511 |
|
| 1512 |
return |
| 1513 |
end |
| 1514 |
C............................................................ |
| 1515 |
subroutine fill_charge_histos(iladder) |
| 1516 |
|
| 1517 |
include '../common/commontracker.f' |
| 1518 |
include '../common/calib.f' |
| 1519 |
include '../common/common_calibration.f' |
| 1520 |
include '../common/level1.f' |
| 1521 |
include '../common/level2.f' |
| 1522 |
* start loop on planes |
| 1523 |
|
| 1524 |
c parameter (ncorrmax=20000) |
| 1525 |
c real chargexy(nplanes,ncorrmax,2) |
| 1526 |
c integer whichplane,goodtr |
| 1527 |
c common/chargecorr/chargexy,whichplane,goodtr |
| 1528 |
|
| 1529 |
character*10 rzdir |
| 1530 |
|
| 1531 |
1010 format('LADDER',i1) |
| 1532 |
write(rzdir,1010)iladder |
| 1533 |
call HCDIR('//'//rzdir,'') |
| 1534 |
call HCDIR('//PAWC/'//rzdir,'') |
| 1535 |
|
| 1536 |
do ip=1,nplanes !1-->6 |
| 1537 |
|
| 1538 |
do ie=1,min(goodtr,ncorrmax) |
| 1539 |
|
| 1540 |
sx = chargexy(ip,iladder,ie,1) |
| 1541 |
sy = chargexy(ip,iladder,ie,2) |
| 1542 |
|
| 1543 |
if(sx.ne.0..and.sy.ne.0)then |
| 1544 |
|
| 1545 |
call HFILL(id_hi_chargeco+ip*100+4,sx,sy,1.) |
| 1546 |
d= sy - syksx(ip)*sx - sycsx(ip) |
| 1547 |
d=d/sqrt(syksx(ip)**2+1) |
| 1548 |
call HFILL(id_hi_chargeco+ip*100, |
| 1549 |
$ d,0.,1.) |
| 1550 |
|
| 1551 |
endif |
| 1552 |
|
| 1553 |
enddo |
| 1554 |
enddo |
| 1555 |
|
| 1556 |
return |
| 1557 |
end |
| 1558 |
C............................................................ |
| 1559 |
|
| 1560 |
subroutine book_check_histos |
| 1561 |
* |
| 1562 |
* some histos to check the sample tracks |
| 1563 |
* and define quality cuts |
| 1564 |
* |
| 1565 |
include '../common/commontracker.f' |
| 1566 |
include '../common/calib.f' |
| 1567 |
include '../common/common_calibration.f' |
| 1568 |
|
| 1569 |
character*64 title !histos title |
| 1570 |
real xbins(nangbin+1) |
| 1571 |
|
| 1572 |
parameter (nrigbin=500) |
| 1573 |
real xrigbin(nrigbin+1) |
| 1574 |
|
| 1575 |
* * * * * * * * * * * * * * * * * * * * * * * * |
| 1576 |
* CHI**2 + RIGIDITY |
| 1577 |
|
| 1578 |
|
| 1579 |
rmin=0.2 |
| 1580 |
rmax=1000. |
| 1581 |
step=abs(log10(rmin)-log10(rmax))/nrigbin |
| 1582 |
xrigbin(1)=rmin |
| 1583 |
do i=2,nrigbin+1 |
| 1584 |
xrigbin(i)=10**(log10(xrigbin(i-1))+step) |
| 1585 |
* print*,'### ',i,xrigbin(i-1),xrigbin(i) |
| 1586 |
enddo |
| 1587 |
call HBOOKB(id_hi_rig |
| 1588 |
$ ,'Rigidity' |
| 1589 |
$ ,nrigbin,xrigbin,0.) |
| 1590 |
|
| 1591 |
|
| 1592 |
call HBOOK1(id_hi_chi |
| 1593 |
$ ,'Chi**2' |
| 1594 |
$ ,1000,0.,100.,0.) |
| 1595 |
call HBPROF(id_hi_chi+100 |
| 1596 |
$ ,'Chi**2 vs Rigidity' |
| 1597 |
$ ,150,rmin,0.1*rmax,0.,100.,0.) |
| 1598 |
call HBPROF(id_hi_chi+200 |
| 1599 |
$ ,'Chi**2 vs Angle x' |
| 1600 |
$ ,50,0.,25.,0.,100.,0.) |
| 1601 |
|
| 1602 |
|
| 1603 |
* * * * * * * * * * * * * * * * * * * * * * * * |
| 1604 |
* ANGLE |
| 1605 |
|
| 1606 |
nbin=200. |
| 1607 |
from=angL(1)!angbin(1) |
| 1608 |
to=angR(nangbin)!angbin(nang+1) |
| 1609 |
do ip=1,nplanes |
| 1610 |
call HBOOK1(id_hi_angle+100+ip |
| 1611 |
$ ,'Angle (Y view)' |
| 1612 |
$ ,nbin,from,to,0.) |
| 1613 |
call HBOOK1(id_hi_angle+200+ip |
| 1614 |
$ ,'Angle (X view)' |
| 1615 |
$ ,nbin,from,to,0.) |
| 1616 |
|
| 1617 |
c$$$ call HBOOK1(id_hi_angle_stat+100+ip |
| 1618 |
c$$$ $ ,'Angle (Y view)' |
| 1619 |
c$$$ $ ,nangbin,from,to,0.) |
| 1620 |
c$$$ call HBOOK1(id_hi_angle_stat+200+ip |
| 1621 |
c$$$ $ ,'Angle (X view)' |
| 1622 |
c$$$ $ ,nangbin,from,to,0.) |
| 1623 |
|
| 1624 |
xbins(1)=angL(1) |
| 1625 |
do i=1,nangbin |
| 1626 |
xbins(i+1)=angR(i) |
| 1627 |
enddo |
| 1628 |
call HBOOKB(id_hi_angle_stat+100+ip |
| 1629 |
$ ,'Angle (Y view)' |
| 1630 |
$ ,nangbin,XBINS,0.) |
| 1631 |
call HBOOKB(id_hi_angle_stat+200+ip |
| 1632 |
$ ,'Angle (X view)' |
| 1633 |
$ ,nangbin,XBINS,0.) |
| 1634 |
enddo |
| 1635 |
|
| 1636 |
|
| 1637 |
* * * * * * * * * * * * * * * * * * * * * * * * |
| 1638 |
* MULTIPLICITY |
| 1639 |
511 format('Multiplicity (view ',i2,', angle bin ',i2,' )') |
| 1640 |
512 format('Multiplicity vs angle (view ',i2,' )') |
| 1641 |
|
| 1642 |
do i=1,nviews |
| 1643 |
do ia=1,nangbin |
| 1644 |
write(title,511)i,ia |
| 1645 |
call HBOOK1(id_hi_mult+i*100+ia |
| 1646 |
$ ,title |
| 1647 |
$ ,10,0.5,10.5,0.) |
| 1648 |
enddo |
| 1649 |
write(title,512)i |
| 1650 |
call HBOOK2(id_hi_mult+i*100 |
| 1651 |
$ ,title |
| 1652 |
$ ,50,0.,25.,10,0.5,10.5,0.) |
| 1653 |
call HBOOK2(80000+id_hi_mult+i*100 |
| 1654 |
$ ,title |
| 1655 |
$ ,50,0.,25.,10,0.5,10.5,0.) |
| 1656 |
enddo |
| 1657 |
* * * * * * * * * * * * * * * * * * * * * * * * |
| 1658 |
* AVERAGE CLUSTER |
| 1659 |
510 format('Average cluster (view ',i2,', angle bin ',i2,' )') |
| 1660 |
|
| 1661 |
do i=1,nviews |
| 1662 |
do ia=1,nangbin |
| 1663 |
write(title,510)i,ia |
| 1664 |
call HBOOK1(id_hi_mult+i*100+ia+50000 |
| 1665 |
$ ,title |
| 1666 |
$ ,int(2*nhside+1) |
| 1667 |
$ ,-(float(nhside)+0.5),(float(nhside)+0.5),0.) |
| 1668 |
enddo |
| 1669 |
enddo |
| 1670 |
|
| 1671 |
|
| 1672 |
* * * * * * * * * * * * * * * * * * * * * * * * |
| 1673 |
* CHARGE |
| 1674 |
|
| 1675 |
5111 format('dE/dX scaled to 300 micron (view ',i2,' )') |
| 1676 |
5112 format('dE/dX vs R (view ',i2,' )') |
| 1677 |
5113 format('dE/dX vs multiplicity (view ',i2,' )') |
| 1678 |
do i=1,nviews |
| 1679 |
write(title,5111)i |
| 1680 |
call HBOOK1(id_hi_charge+i*100 |
| 1681 |
$ ,title |
| 1682 |
$ ,500,0.,1000.,0.) |
| 1683 |
write(title,5112)i |
| 1684 |
call HBPROF(id_hi_charge+i*100+1 |
| 1685 |
$ ,title |
| 1686 |
$ ,1000,0.,100.,0.,1000.,0.) |
| 1687 |
write(title,5113)i |
| 1688 |
call HBOOK2(id_hi_charge+i*100+2 |
| 1689 |
$ ,title |
| 1690 |
$ ,100,0.,1000.,10,0.5,10.5,0.) |
| 1691 |
call HBOOK2(80000+id_hi_charge+i*100+2 |
| 1692 |
$ ,title |
| 1693 |
$ ,100,0.,1000.,10,0.5,10.5,0.) |
| 1694 |
enddo |
| 1695 |
call HBOOK1(id_hi_charge |
| 1696 |
$ ,'dE/dX scaled to 300 micron AVERAGE' |
| 1697 |
$ ,500,0.,1000.,0.) |
| 1698 |
|
| 1699 |
* * * * * * * * * * * * * * * * * * * * * * * * |
| 1700 |
* EVENT SPATIAL DISTRIBUTION |
| 1701 |
|
| 1702 |
51111 format('Strip Y (plane ',i2,', sensor ',i1,')') |
| 1703 |
51112 format('Strip X (plane ',i2,', sensor ',i1,')') |
| 1704 |
do ip=1,nplanes |
| 1705 |
do is=1,2 |
| 1706 |
write(title,51111)ip,is |
| 1707 |
call HBOOK1(id_hi_beam+ip*100+is+100000 |
| 1708 |
$ ,title |
| 1709 |
$ ,256,0.5,3072.5,0.) |
| 1710 |
write(title,51112)ip,is |
| 1711 |
call HBOOK1(id_hi_beam+ip*100+is+200000 |
| 1712 |
$ ,title |
| 1713 |
$ ,256,0.5,3072.5,0.) |
| 1714 |
enddo |
| 1715 |
enddo |
| 1716 |
|
| 1717 |
* * * * * * * * * * * * * * * * * * * * * * * * |
| 1718 |
* MAGNETIC FIELD |
| 1719 |
|
| 1720 |
61 format('Bx (plane ',i2,')') |
| 1721 |
62 format('By (plane ',i2,')') |
| 1722 |
63 format('Bz (plane ',i2,')') |
| 1723 |
do ip=1,nplanes |
| 1724 |
write(title,61)ip |
| 1725 |
call HBOOK1(id_hi_bx+ip |
| 1726 |
$ ,title |
| 1727 |
$ ,1000,-5.,5.,0.) |
| 1728 |
write(title,62)ip |
| 1729 |
call HBOOK1(id_hi_by+ip |
| 1730 |
$ ,title |
| 1731 |
$ ,1000,-5.,5.,0.) |
| 1732 |
write(title,63)ip |
| 1733 |
call HBOOK1(id_hi_bz+ip |
| 1734 |
$ ,title |
| 1735 |
$ ,1000,-5.,5.,0.) |
| 1736 |
enddo |
| 1737 |
|
| 1738 |
|
| 1739 |
|
| 1740 |
return |
| 1741 |
end |
| 1742 |
C............................................................ |
| 1743 |
subroutine fill_check_histos |
| 1744 |
|
| 1745 |
include '../common/commontracker.f' |
| 1746 |
include '../common/calib.f' |
| 1747 |
include '../common/common_calibration.f' |
| 1748 |
include '../common/common_mech.f' |
| 1749 |
include '../common/level1.f' |
| 1750 |
include '../common/level2.f' |
| 1751 |
|
| 1752 |
real v(3),b(3) |
| 1753 |
character*10 rzdir |
| 1754 |
|
| 1755 |
real clusterx(2*nhside+1) |
| 1756 |
real clustery(2*nhside+1) |
| 1757 |
logical DRAY,LMSX,LMSY |
| 1758 |
|
| 1759 |
*--------------------------------------------------------- |
| 1760 |
* track info parameters |
| 1761 |
*--------------------------------------------------------- |
| 1762 |
rig = abs(1./al_nt(5,1)) |
| 1763 |
chi = chi2_nt(1) |
| 1764 |
|
| 1765 |
|
| 1766 |
|
| 1767 |
1010 format('LADDER',i1) |
| 1768 |
|
| 1769 |
* i plot della rigidita` del chi2 |
| 1770 |
* non dipendono dal ladder... |
| 1771 |
do il=1,3 |
| 1772 |
write(rzdir,1010)il |
| 1773 |
call HCDIR('//'//rzdir,'') |
| 1774 |
call HCDIR('//PAWC/'//rzdir,'') |
| 1775 |
call HFILL(id_hi_chi,chi,0.,1.) |
| 1776 |
call HFILL(id_hi_rig,rig,0.,1.) |
| 1777 |
thetax = axv_nt(1,1) |
| 1778 |
thetay = ayv_nt(1,1) |
| 1779 |
call HFILL(id_hi_chi+100,chi,rig,1.) |
| 1780 |
call HFILL(id_hi_chi+200,chi,thetax,1.) |
| 1781 |
|
| 1782 |
enddo |
| 1783 |
|
| 1784 |
dedxscalav = 0 |
| 1785 |
npdedxscalav = 0 |
| 1786 |
* start loop on planes |
| 1787 |
do ip=1,nplanes !TOP-->BOTTOM |
| 1788 |
*--------------------------------------------------------- |
| 1789 |
* main x-y cluster parameters |
| 1790 |
*--------------------------------------------------------- |
| 1791 |
iviewx = nviewx(7-ip) |
| 1792 |
iviewy = nviewy(7-ip) |
| 1793 |
xxx = xm_nt(ip,1) |
| 1794 |
yyy = ym_nt(ip,1) |
| 1795 |
call whichsensor((7-ip),xxx,yyy,iladder,isensor) |
| 1796 |
if(iladder.eq.0.or.isensor.eq.0)then |
| 1797 |
print*,'*** out *** ',ip,xxx,yyy,iladder,isensor |
| 1798 |
goto 10101 |
| 1799 |
endif |
| 1800 |
multx = mult(which_cl_x(ip)) |
| 1801 |
multy = mult(which_cl_y(ip)) |
| 1802 |
thetax = axv_nt(ip,1) |
| 1803 |
thetay = ayv_nt(ip,1) |
| 1804 |
* =============================================== |
| 1805 |
* looser multiplicity cut ... after LANDI warning |
| 1806 |
* (I also require that both clx anc cly are low-multiplicity) |
| 1807 |
LMSX=.false. |
| 1808 |
LMSY=.false. |
| 1809 |
if( |
| 1810 |
$ (abs(thetax).le.10..and.multx.le.3).or. |
| 1811 |
$ ((abs(thetax).gt.10..and.abs(thetax).le.15.).and.multx.le.4).or. |
| 1812 |
$ (abs(thetax).gt.15..and.multx.le.5).or. |
| 1813 |
$ .false.)LMSX=.true. |
| 1814 |
if( |
| 1815 |
$ (multy.le.3).or. |
| 1816 |
$ .false.)LMSY=.true. |
| 1817 |
|
| 1818 |
* charge scaling factor (to vertical trajectory) |
| 1819 |
pi=acos(-1.) |
| 1820 |
fact= |
| 1821 |
$ sqrt( |
| 1822 |
$ tan(pi*thetax/180.)**2 + |
| 1823 |
$ tan(pi*thetay/180.)**2 + |
| 1824 |
$ 1. |
| 1825 |
$ ) |
| 1826 |
dedxscalx = dedx(which_cl_x(ip)) / fact |
| 1827 |
dedxscaly = dedx(which_cl_y(ip)) / fact |
| 1828 |
dedxscal = (dedxscalx+dedxscaly)/2 |
| 1829 |
DRAY=.false. |
| 1830 |
if(dedxscal.gt.dedxmax)DRAY=.true. |
| 1831 |
|
| 1832 |
if(.not.DRAY.and.LMSX.and.LMSX)then |
| 1833 |
dedxscalav = dedxscalav + dedxscal |
| 1834 |
npdedxscalav = npdedxscalav + 1 |
| 1835 |
endif |
| 1836 |
|
| 1837 |
|
| 1838 |
|
| 1839 |
ista = indstart(which_cl_x(ip)) |
| 1840 |
ice = indmax(which_cl_x(ip)) |
| 1841 |
ilast = totCLlength |
| 1842 |
if(which_cl_x(ip).ne.nclstr1) |
| 1843 |
$ ilast = indstart(which_cl_x(ip)+1)-1 |
| 1844 |
do i = -nhside,nhside |
| 1845 |
ii = ice + i |
| 1846 |
iii=i+nhside+1 |
| 1847 |
clusterx(iii) = 0 |
| 1848 |
if(ii.le.ilast.and.ii.ge.ista) |
| 1849 |
$ clusterx(iii) = clsignal(ii) / dedx(which_cl_x(ip)) |
| 1850 |
enddo |
| 1851 |
|
| 1852 |
ista = indstart(which_cl_y(ip)) |
| 1853 |
ice = indmax(which_cl_y(ip)) |
| 1854 |
ilast = totCLlength |
| 1855 |
if(which_cl_y(ip).ne.nclstr1) |
| 1856 |
$ ilast = indstart(which_cl_y(ip)+1)-1 |
| 1857 |
do i = -nhside,nhside |
| 1858 |
ii = ice + i |
| 1859 |
iii=i+nhside+1 |
| 1860 |
clustery(iii) = 0 |
| 1861 |
if(ii.le.ilast.and.ii.ge.ista) |
| 1862 |
$ clustery(iii) = clsignal(ii) / dedx(which_cl_y(ip)) |
| 1863 |
enddo |
| 1864 |
|
| 1865 |
*--------------------------------------------------------- |
| 1866 |
|
| 1867 |
write(rzdir,1010)iladder |
| 1868 |
call HCDIR('//'//rzdir,'') |
| 1869 |
call HCDIR('//PAWC/'//rzdir,'') |
| 1870 |
|
| 1871 |
|
| 1872 |
|
| 1873 |
*--------------------------------------------------------- |
| 1874 |
* find the angular bins (ibin_angy ibin_angx) |
| 1875 |
* associated to the track intersection point with plane ip |
| 1876 |
*--------------------------------------------------------- |
| 1877 |
ia=0 |
| 1878 |
if( thetax.lt.angL(1).or. |
| 1879 |
$ thetax.ge.angR(nangbin))goto 111 |
| 1880 |
do ia=1,nangbin |
| 1881 |
if (thetay.lt.angR(ia))goto 111 |
| 1882 |
enddo |
| 1883 |
111 ibin_angy = ia |
| 1884 |
|
| 1885 |
ia=0 |
| 1886 |
if( thetay.lt.angL(1).or. |
| 1887 |
$ thetay.ge.angR(nangbin))goto 112 |
| 1888 |
do ia=1,nangbin |
| 1889 |
if (thetax.lt.angR(ia))goto 112 |
| 1890 |
enddo |
| 1891 |
112 ibin_angx = ia |
| 1892 |
|
| 1893 |
|
| 1894 |
|
| 1895 |
*--------------------------------------------------------- |
| 1896 |
* fill histos |
| 1897 |
*--------------------------------------------------------- |
| 1898 |
|
| 1899 |
************************************************ ANGLE |
| 1900 |
|
| 1901 |
call HFILL(id_hi_angle+200+ip,thetax,0.,1.) |
| 1902 |
call HFILL(id_hi_angle+100+ip,thetay,0.,1.) |
| 1903 |
|
| 1904 |
call HFILL(id_hi_angle_stat+200+ip,thetax,0.,1.) |
| 1905 |
call HFILL(id_hi_angle_stat+100+ip,thetay,0.,1.) |
| 1906 |
|
| 1907 |
************************************************ multiplicity & cluster |
| 1908 |
if(ibin_angx.ne.0)then |
| 1909 |
call HFILL(id_hi_mult+iviewx*100+ibin_angx, |
| 1910 |
$ float(mult(which_cl_x(ip))),0.,1.) |
| 1911 |
do i=-nhside,nhside |
| 1912 |
ii = i + nhside +1 |
| 1913 |
xxxx=clusterx(ii) |
| 1914 |
call HFILL(id_hi_mult+iviewx*100+ibin_angx+50000, |
| 1915 |
$ float(i),0.,xxxx) |
| 1916 |
enddo |
| 1917 |
endif |
| 1918 |
|
| 1919 |
if(ibin_angy.ne.0)then |
| 1920 |
call HFILL(id_hi_mult+iviewy*100+ibin_angy, |
| 1921 |
$ float(mult(which_cl_y(ip))),0.,1.) |
| 1922 |
do i=-nhside,nhside |
| 1923 |
ii = i + nhside +1 |
| 1924 |
xxxx=clustery(ii) |
| 1925 |
call HFILL(id_hi_mult+iviewy*100+ibin_angy+50000, |
| 1926 |
$ float(i),0.,xxxx) |
| 1927 |
enddo |
| 1928 |
endif |
| 1929 |
|
| 1930 |
************************************************ charge |
| 1931 |
call HFILL(id_hi_charge+iviewx*100, |
| 1932 |
$ dedxscalx,0.,1.) |
| 1933 |
call HFILL(id_hi_charge+iviewy*100, |
| 1934 |
$ dedxscaly,0.,1.) |
| 1935 |
|
| 1936 |
call HFILL(id_hi_charge+iviewx*100+1, |
| 1937 |
$ rig,dedxscalx,1.) |
| 1938 |
call HFILL(id_hi_charge+iviewy*100+1, |
| 1939 |
$ rig,dedxscaly,1.) |
| 1940 |
|
| 1941 |
* -------------------------------------------- |
| 1942 |
C this are the condition applied to select the |
| 1943 |
C particle sample for pfa calibration |
| 1944 |
* -------------------------------------------- |
| 1945 |
if(.not.DRAY.and.LMSX.and.LMSY)then |
| 1946 |
call HFILL(id_hi_charge+iviewx*100+2, |
| 1947 |
$ dedxscalx,float(multx),1.) |
| 1948 |
call HFILL(id_hi_charge+iviewy*100+2, |
| 1949 |
$ dedxscaly,float(multy),1.) |
| 1950 |
call HFILL(id_hi_mult+iviewx*100, |
| 1951 |
$ abs(thetax),float(multx),1.) |
| 1952 |
call HFILL(id_hi_mult+iviewy*100, |
| 1953 |
$ abs(thetay),float(multy),1.) |
| 1954 |
endif |
| 1955 |
call HFILL(80000+id_hi_charge+iviewx*100+2, |
| 1956 |
$ dedxscalx,float(multx),1.) |
| 1957 |
call HFILL(80000+id_hi_charge+iviewy*100+2, |
| 1958 |
$ dedxscaly,float(multy),1.) |
| 1959 |
call HFILL(80000+id_hi_mult+iviewx*100, |
| 1960 |
$ abs(thetax),float(multx),1.) |
| 1961 |
call HFILL(80000+id_hi_mult+iviewy*100, |
| 1962 |
$ abs(thetay),float(multy),1.) |
| 1963 |
|
| 1964 |
************************************************ beam profle |
| 1965 |
call hfill(id_hi_beam+ip*100+isensor+100000, |
| 1966 |
$ float(MAXS(which_cl_y(ip))),0.,1.) |
| 1967 |
call hfill(id_hi_beam+ip*100+isensor+200000, |
| 1968 |
$ float(MAXS(which_cl_x(ip))),0.,1.) |
| 1969 |
|
| 1970 |
|
| 1971 |
************************************************ magnetic field |
| 1972 |
|
| 1973 |
v(1) = xv_nt(ip,1) |
| 1974 |
v(2) = yv_nt(ip,1) |
| 1975 |
v(3) = zv_nt(ip,1) |
| 1976 |
|
| 1977 |
call gufld(v,b) |
| 1978 |
|
| 1979 |
bx=b(1) |
| 1980 |
by=b(2) |
| 1981 |
bz=b(3) |
| 1982 |
call hfill(id_hi_bx+ip,bx,0.,1.) |
| 1983 |
call hfill(id_hi_by+ip,by,0.,1.) |
| 1984 |
call hfill(id_hi_bz+ip,bz,0.,1.) |
| 1985 |
|
| 1986 |
enddo |
| 1987 |
|
| 1988 |
dedxscalav = dedxscalav/npdedxscalav |
| 1989 |
call HFILL(id_hi_charge, |
| 1990 |
$ dedxscalav,0.,1.) |
| 1991 |
|
| 1992 |
|
| 1993 |
|
| 1994 |
10101 continue |
| 1995 |
call HCDIR('//PAWC','') |
| 1996 |
|
| 1997 |
return |
| 1998 |
end |
| 1999 |
|
| 2000 |
*************************************************************************** |
| 2001 |
* |
| 2002 |
* |
| 2003 |
* |
| 2004 |
* |
| 2005 |
* |
| 2006 |
* |
| 2007 |
* |
| 2008 |
* |
| 2009 |
* |
| 2010 |
*************************************************************************** |
| 2011 |
|
| 2012 |
|
| 2013 |
subroutine FCNcharge(npar,grad,fval,xval,iflag,futil) |
| 2014 |
|
| 2015 |
include '../common/commontracker.f' |
| 2016 |
include '../common/calib.f' |
| 2017 |
include '../common/common_calibration.f' |
| 2018 |
|
| 2019 |
double precision grad(*),xval(*),fval |
| 2020 |
integer npar,iflag |
| 2021 |
double precision KKK,CCC |
| 2022 |
|
| 2023 |
|
| 2024 |
c parameter (ncorrmax=20000) |
| 2025 |
c real chargexy(nplanes,ncorrmax,2) |
| 2026 |
c integer whichplane,goodtr |
| 2027 |
c common/chargecorr/chargexy,whichplane,goodtr |
| 2028 |
|
| 2029 |
KKK=xval(1) |
| 2030 |
CCC=xval(2) |
| 2031 |
|
| 2032 |
|
| 2033 |
fval = 0. |
| 2034 |
|
| 2035 |
|
| 2036 |
sum2=0 |
| 2037 |
do i=1,min(goodtr,ncorrmax) |
| 2038 |
sx=chargexy(whichplane,whichladder,i,1) |
| 2039 |
sy=chargexy(whichplane,whichladder,i,2) |
| 2040 |
* if(sx.gt.140..and.sy.gt.140)then |
| 2041 |
d=(sy-KKK*sx-CCC)**2 |
| 2042 |
d=d/(KKK**2+1) |
| 2043 |
d=d/(1.9*4.**2+1.2*8.**2) !errors |
| 2044 |
sum2=sum2+d |
| 2045 |
c$$$ PRINT*,'xx ',i,d,' sx-sy ' |
| 2046 |
c$$$ $ ,chargexy(whichplane,whichladder,i,1) |
| 2047 |
c$$$ $ ,chargexy(whichplane,whichladder,i,2) |
| 2048 |
* endif |
| 2049 |
enddo |
| 2050 |
|
| 2051 |
|
| 2052 |
fval = sum2 |
| 2053 |
|
| 2054 |
|
| 2055 |
|
| 2056 |
return |
| 2057 |
|
| 2058 |
end |
| 2059 |
|
| 2060 |
|
| 2061 |
|
| 2062 |
|
| 2063 |
************************************************** |
| 2064 |
************************************************** |
| 2065 |
* CHARGE CORRELATION CALIBRATION |
| 2066 |
************************************************** |
| 2067 |
************************************************** |
| 2068 |
subroutine charge_calibration(iladder) |
| 2069 |
|
| 2070 |
include '../common/commontracker.f' |
| 2071 |
include '../common/calib.f' |
| 2072 |
include '../common/common_calibration.f' |
| 2073 |
c include '../common/common_mech.f' |
| 2074 |
c include '../common/level1.f' |
| 2075 |
c include '../common/level2.f' |
| 2076 |
external FCNcharge |
| 2077 |
|
| 2078 |
character*12 chnam |
| 2079 |
double precision val,error |
| 2080 |
|
| 2081 |
|
| 2082 |
do ip=1,nplanes |
| 2083 |
print*,' ' |
| 2084 |
print*,' ' |
| 2085 |
print*,' ' |
| 2086 |
print*, |
| 2087 |
$ ' -----------------------------------' |
| 2088 |
print*, |
| 2089 |
$ '////////////// MINUIT for plane ',ip,' - ladder ',iladder |
| 2090 |
print*, |
| 2091 |
$ ' -----------------------------------' |
| 2092 |
open(UNIT=40,FILE= |
| 2093 |
$ 'bin-aux/data-card.dat' |
| 2094 |
$ ,STATUS='OLD') |
| 2095 |
call mintio(40,6,7) |
| 2096 |
whichplane = ip |
| 2097 |
whichladder = iladder |
| 2098 |
call minuit(FCNcharge,0) |
| 2099 |
c call MNERRS(1,eplus,eminus,eparab,globcc) |
| 2100 |
call MNPOUT(1,chnam,val,error,bnd1,bnd2,ivarbl) |
| 2101 |
syksx(ip)=val |
| 2102 |
esyksx(ip)=error |
| 2103 |
call MNPOUT(2,chnam,val,error,bnd1,bnd2,ivarbl) |
| 2104 |
sycsx(ip)=val |
| 2105 |
esycsx(ip)=error |
| 2106 |
close(40) |
| 2107 |
enddo |
| 2108 |
|
| 2109 |
return |
| 2110 |
end |
| 2111 |
|
| 2112 |
************************************************** |
| 2113 |
************************************************** |
| 2114 |
* ETA CALIBRATION |
| 2115 |
************************************************** |
| 2116 |
************************************************** |
| 2117 |
subroutine eta_calibration(iladder,ia) |
| 2118 |
include '../common/commontracker.f' |
| 2119 |
include '../common/calib.f' |
| 2120 |
include '../common/common_calibration.f' |
| 2121 |
|
| 2122 |
* -------------------------------------------- |
| 2123 |
* some local variables for eta correction |
| 2124 |
real histo2(nhbin2) |
| 2125 |
c real eta2corr(nhbin2+1,nviews) !eta2 correction |
| 2126 |
c real eta2abs(nhbin2+1) !eta2 abscissa |
| 2127 |
* -------------------------------------------- |
| 2128 |
character*64 ctemp |
| 2129 |
character*10 rzdir |
| 2130 |
|
| 2131 |
1010 format('LADDER',i1) |
| 2132 |
write(rzdir,1010)iladder |
| 2133 |
call HCDIR('//'//rzdir,'') |
| 2134 |
call HCDIR('//PAWC/'//rzdir,'') |
| 2135 |
|
| 2136 |
do i=1,nviews |
| 2137 |
id=id_hi_eta2+i*100+ia |
| 2138 |
* eta2 histogram |
| 2139 |
call HUNPAK(id,histo2,' ',0) !eta2 histogram |
| 2140 |
call HGIVE(id,ctemp,nx,xmi,xma,ny,ymi,yma,nwt,loc) |
| 2141 |
sum=HSUM(id) |
| 2142 |
c print*,ctemp,nx,xmi,xma,ny,ymi,yma,nwt,loc |
| 2143 |
if(i.eq.1)then |
| 2144 |
do ie=0,nx |
| 2145 |
eta2abs(ie+1)=xmi+ie*(xma-xmi)/nx |
| 2146 |
enddo |
| 2147 |
endif |
| 2148 |
if(sum.eq.0)then |
| 2149 |
do ie=0,nx |
| 2150 |
eta2corr(ie+1,i)=eta2abs(ie+1) |
| 2151 |
enddo |
| 2152 |
goto 3 |
| 2153 |
endif |
| 2154 |
* ------------------------------------------------------ |
| 2155 |
* if the histogram contains less than NMIN4ETA2 entries, |
| 2156 |
* the correction factor is not calculated |
| 2157 |
* CHANGED!!! |
| 2158 |
* it is calculated anyway.... |
| 2159 |
* ------------------------------------------------------ |
| 2160 |
nmin4eta2 = 200 |
| 2161 |
if(sum.lt.nmin4eta2)then |
| 2162 |
c$$$ do ie=0,nx |
| 2163 |
c$$$ eta2corr(ie+1,i)=eta2abs(ie+1) |
| 2164 |
c$$$ enddo |
| 2165 |
print*,'*** WARNING ***' |
| 2166 |
print*,'ETA2 parameters: less than ',nmin4eta2 |
| 2167 |
$ ,' entries ----> ',rzdir,' view ',i,' ang. bin ',ia |
| 2168 |
c$$$ goto 3 |
| 2169 |
endif |
| 2170 |
* integral of the eta2 histogram |
| 2171 |
eta2corr(1,i)=0. |
| 2172 |
do ie=1,nx |
| 2173 |
eta2corr(ie+1,i)=eta2corr(ie,i)+histo2(ie) |
| 2174 |
enddo |
| 2175 |
* eta2 normalization |
| 2176 |
angle=avangbin(ia,i) |
| 2177 |
x0=eta2fix(i,angle) !point of reference for correction function |
| 2178 |
do ifix=1,nx+1 |
| 2179 |
if(eta2abs(ifix).ge.x0)goto 2 |
| 2180 |
enddo |
| 2181 |
2 y0=eta2corr(ifix,i)/sum |
| 2182 |
c shift=eta2shift(i,angle) |
| 2183 |
shift = 0. |
| 2184 |
do ie=0,nx |
| 2185 |
eta2corr(ie+1,i)= |
| 2186 |
$ x0-y0+eta2corr(ie+1,i)/sum-shift |
| 2187 |
enddo |
| 2188 |
3 continue |
| 2189 |
enddo !end loop on views |
| 2190 |
|
| 2191 |
|
| 2192 |
|
| 2193 |
return |
| 2194 |
end |
| 2195 |
*---***---***---***---***---***---***---***---*** |
| 2196 |
* |
| 2197 |
* functions to define the point where to "fix" |
| 2198 |
* the eta correction |
| 2199 |
* and the shift |
| 2200 |
* |
| 2201 |
*---***---***---***---***---***---***---***---*** |
| 2202 |
function eta2fix(iview,angle) |
| 2203 |
|
| 2204 |
if(mod(iview,2).eq.0)then !>>>>>>>> X |
| 2205 |
|
| 2206 |
eta2fix = -0.5 |
| 2207 |
|
| 2208 |
else !>>>>>>>> Y |
| 2209 |
|
| 2210 |
eta2fix = -0.5 |
| 2211 |
|
| 2212 |
endif |
| 2213 |
|
| 2214 |
|
| 2215 |
return |
| 2216 |
end |
| 2217 |
|
| 2218 |
|
| 2219 |
c$$$ function eta2shift(iview,angle) |
| 2220 |
c$$$ |
| 2221 |
c$$$ parameter (hallm_e=1650.) !cm**2 / Vs |
| 2222 |
c$$$ parameter (hallm_h=310.) |
| 2223 |
c$$$ |
| 2224 |
c$$$ include '../common/commontracker.f' |
| 2225 |
c$$$ include '../common/common_preanalysis.f' |
| 2226 |
c$$$ character*32 ctemp |
| 2227 |
c$$$ real hallm |
| 2228 |
c$$$ |
| 2229 |
c$$$ pi=acos(-1.) |
| 2230 |
c$$$ |
| 2231 |
c$$$ np=npl(iview) |
| 2232 |
c$$$ if(mod(iview,2).eq.0)then !>>>>>>>> X |
| 2233 |
c$$$ id=id_hi_by+np |
| 2234 |
c$$$ hallm=hallm_h*1e-4 |
| 2235 |
c$$$ else |
| 2236 |
c$$$ id=id_hi_bx+np |
| 2237 |
c$$$ hallm=hallm_e*1e-4 |
| 2238 |
c$$$ endif |
| 2239 |
c$$$ bfield=HSTATI(id,1,' ',0) |
| 2240 |
c$$$ bfield=bfield*1e-1 |
| 2241 |
c$$$ tgth_h=-hallm*bfield |
| 2242 |
c$$$ |
| 2243 |
c$$$ tgth_eff=tan(pi*angle/180.)-tgth_h |
| 2244 |
c$$$ theta_eff=atan(tgth_eff)*180./pi |
| 2245 |
c$$$ if(iview.ne.12)theta_eff=-theta_eff |
| 2246 |
c$$$ |
| 2247 |
c$$$ if(mod(iview,2).eq.0)then !>>>>>>>> X |
| 2248 |
c$$$ |
| 2249 |
c$$$ eta2shift = -1.549*sin(2*pi*theta_eff/13.59) |
| 2250 |
c$$$ eta2shift = eta2shift/51. |
| 2251 |
c$$$ |
| 2252 |
c$$$ else !>>>>>>>> Y |
| 2253 |
c$$$ |
| 2254 |
c$$$ eta2shift = -1.348*sin(2*pi*theta_eff/27.73) |
| 2255 |
c$$$ eta2shift = eta2shift/67. |
| 2256 |
c$$$* NBNBNBNBNB ... da risolvere il problema della simulazione!!!! |
| 2257 |
c$$$ |
| 2258 |
c$$$ endif |
| 2259 |
c$$$ |
| 2260 |
c$$$c print*,iview,np,bfield,angle,theta_eff,'>>>>>>> ',eta2shift |
| 2261 |
c$$$ |
| 2262 |
c$$$ return |
| 2263 |
c$$$ end |
Parent Directory
| 
Revision Log