git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixed Warnings in Reconstructioner and Prob
[u/mrichter/AliRoot.git] / TOF / AliTOFReconstructioner.cxx
CommitLineData
db9ba97f 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//_________________________________________________________________________
17// Manager class for TOF reconstruction.
18//
19//
20//-- Authors: Bologna-ITEP-Salerno Group
21//
22// Description: Manager class for TOF reconstruction (derived from TTask)
23// Summary of the main methods:
24// - extraction of the TPC (assumed to be) reconstructed tracks
25// comment: it has to me moved as soon as possible into a separate
26// class AliTOFTrackReader (K. Safarik suggestion)
27// - geometrical propagation of the above tracks till TOF detector
28// - matching of the tracks with the TOF signals
29//
30// Remark: the GEANT3.21 geometry is used during the geometrical propagation
31// of the tracks in order to know the current volume reached by the track.
32//
33//////////////////////////////////////////////////////////////////////////////
34
35
db9ba97f 36#include "AliConst.h"
37#include "AliRun.h"
38#include "AliTOFConstants.h"
39#include "AliTOFHitMap.h"
40#include "AliTOFSDigit.h"
41#include "AliTOFhit.h"
42#include "AliTOFRecHit.h"
43#include "AliTOFPad.h"
44#include "AliTOFTrack.h"
45#include "AliTOF.h"
46#include "AliTOFv1.h"
47#include "AliTOFv2.h"
48#include "AliTOFv2FHoles.h"
49#include "AliTOFv3.h"
50#include "AliTOFv4.h"
51#include "AliTOFv4T0.h"
52#include "AliTOFReconstructioner.h"
53// this line has to be commented till TPC will provide fPx fPy fPz and fL in
54// AliTPChit class or somewhere
55// #include "../TPC/AliTPC.h"
56#include "AliRun.h"
57#include "AliDetector.h"
58#include "AliMC.h"
59
b213b8bd 60#include "TTask.h"
61#include "TBenchmark.h"
62#include "TTree.h"
63#include "TSystem.h"
64#include "TFile.h"
65#include "TParticle.h"
db9ba97f 66#include <TClonesArray.h>
b9d0a01d 67#include "TGeant3.h"
68#include "TVirtualMC.h"
db9ba97f 69#include <TF1.h>
70#include <TF2.h>
db9ba97f 71#include "TROOT.h"
72#include "TFolder.h"
73#include "TNtuple.h"
74#include <stdlib.h>
75#include <iostream.h>
76#include <fstream.h>
77
78ClassImp(AliTOFReconstructioner)
79
80//____________________________________________________________________________
81 AliTOFReconstructioner::AliTOFReconstructioner():TTask("AliTOFReconstructioner","")
82{
83 // default ctor
84 fNevents = 0 ;
db9ba97f 85 foutputfile = 0;
86 foutputntuple= 0;
87 fZnoise = 0;
88 ftail = 0;
89}
90
91//____________________________________________________________________________
92 AliTOFReconstructioner::AliTOFReconstructioner(char* headerFile, Option_t* opt, char *RecFile ):TTask("AliTOFReconstructioner","")
93{
94 //
95 // ctor
96 //
97 fNevents = 0 ; // Number of events to reconstruct, 0 means all evens in current file
db9ba97f 98 foutputfile = 0;
99 foutputntuple= 0;
100 fZnoise = 0;
101 ftail = 0;
102
103 Init(opt);
104
105 // create output file
106 if (RecFile){
107 foutputfile= new TFile(RecFile,"RECREATE","root file for matching");
108 } else {
109 char outFileName[100];
110 strcpy(outFileName,"match");
111 strcat(outFileName,headerFile);
112 foutputfile= new TFile(outFileName,"RECREATE","root file for matching");
113 }
114
115 // initialize the ALIROOT geometry
116 gAlice->Init();
117 gAlice->Print();
118
db9ba97f 119 CreateNTuple();
120
121 // add Task to //root/Tasks folder
122 TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
123 roottasks->Add(this) ;
124}
125//____________________________________________________________________________
126void AliTOFReconstructioner::Init(Option_t* opt)
127{
128 // Initialize the AliTOFReconstructioner setting parameters for
129 // reconstruction.
130 // Option values: Pb-Pb for Pb-Pb events
131 // pp for pp events
132
133 // set common parameters
134 fdbg=1;
135 fNevents = 1;
136 fFirstEvent = 1;
137 fLastEvent = 1;
138 fTimeResolution =0.120;
139 fpadefficiency =0.99 ;
140 fEdgeEffect = 2 ;
141 fEdgeTails = 0 ;
142 fHparameter = 0.4 ;
143 fH2parameter = 0.15;
144 fKparameter = 0.5 ;
145 fK2parameter = 0.35;
146 fEffCenter = fpadefficiency;
147 fEffBoundary = 0.65;
148 fEff2Boundary = 0.90;
149 fEff3Boundary = 0.08;
150 fResCenter = 50. ;
151 fResBoundary = 70. ;
152 fResSlope = 40. ;
153 fTimeWalkCenter = 0. ;
154 fTimeWalkBoundary=0. ;
155 fTimeWalkSlope = 0. ;
156 fTimeDelayFlag = 1 ;
157 fPulseHeightSlope=2.0 ;
158 fTimeDelaySlope =0.060;
159 // was fMinimumCharge = TMath::Exp(fPulseHeightSlope*fKparameter/2.);
160 fMinimumCharge = TMath::Exp(-fPulseHeightSlope*fHparameter);
161 fChargeSmearing=0.0 ;
162 fLogChargeSmearing=0.13;
163 fTimeSmearing =0.022;
164 fAverageTimeFlag=0 ;
165 fChargeFactorForMatching=1;
166 fTrackingEfficiency=1.0; // 100% TPC tracking efficiency assumed
167 fSigmavsp = 1. ;
168 fSigmaZ = 0. ;
169 fSigmarphi= 0. ;
170 fSigmap = 0. ;
171 fSigmaPhi = 0. ;
172 fSigmaTheta=0. ;
173 fField = 0.2 ;
174 // fRadLenTPC : 0.2 includes TRD / 0.03 TPC only
175 fRadLenTPC=0.06 ; // last value
176 fCorrectionTRD=0. ;
177 fLastTPCRow=111 ;
178 fRadiusvtxBound=50. ; // expressed in [cm]
179 fStep = 0.1 ; // expressed in [cm] step during propagation of the
180 // track inside TOF volumes
181 fMatchingStyle=2 ;
182 /* previous values default
183 fMaxPixels=70000 ;
184 fMaxAllTracks=70000 ;
185 fMaxTracks=15000 ;
186 */
187 fMaxPixels=165000 ;
188 fMaxAllTracks=500000 ;
189 fMaxTracks=15000 ;
190
191 fMaxTOFHits=35000 ;
192 fPBound =0.0 ; // bending effect: P_t=0.3*z*B*R , z particle charge
193 fNoiseSlope=20. ;
194 // set parameters as specified in opt
195 //pp case
196 if(strstr(opt,"pp")){
197 fMaxTestTracks=500 ;
198 fNoise = 26. ;
199 fNoiseMeanTof= 26.4 ; // to check
200 }
201 //Pb-Pb case
202 if(strstr(opt,"Pb-Pb")){
203 fMaxTestTracks=20 ;
204 fNoise = 9400. ;
205 fNoiseMeanTof= 26.4 ;
206 }
207}
208
209//____________________________________________________________________________
210 AliTOFReconstructioner::~AliTOFReconstructioner()
211{
212 //
213 // dtor
214 //
b9d0a01d 215
db9ba97f 216 if (foutputfile)
217 {
218 delete foutputfile;
219 foutputfile = 0;
220 }
221 if (foutputntuple)
222 {
223 delete foutputntuple;
224 foutputntuple = 0;
225 }
226
227 if (fZnoise)
228 {
229 delete fZnoise;
230 fZnoise = 0;
231 }
232
233 if (ftail)
234 {
235 delete ftail;
236 ftail = 0;
237 }
238}
239
240//____________________________________________________________________________
241void AliTOFReconstructioner::CreateNTuple()
242{
243 //
244 // Create a Ntuple where information about reconstructed charged particles
245 // (both primaries and secondaries) are stored
246 // Variables: event ipart imam xvtx yvtx zvtx pxvtx pyvtx pzvtx time leng matc text mext
247 // Meaning:
248 // event - event number (0, 1, ...)
249 // ipart - PDG code of particles
250 // imam - PDG code for the parent
251 // =0 for primary particle
252 // xvtx - x-coordinate of the vertex (cm)
253 // yvtx - y-coordinate of the vertex (cm)
254 // zvtx - z-coordinate of the vertex (cm)
255 // pxvtx - x-coordinate of the momentum in the vertex (GeV)
256 // pyvtx - y-coordinate of the momentum in the vertex (GeV)
257 // pzvtx - z-coordinate of the momentum in the vertex (GeV)
258 // time - time of flight from TOF for given track (ps) - TOF time for the
259 // first TOF hit of the track
260 // leng - track length to the TOF pixel (cm), evaluate as a sum of the
261 // track length from the track vertex to TPC and the average
262 // length of the extrapolated track from TPC to TOF.
263 // for the track without TOF hits leng=-abs(leng)
264 // matc - index of the (TPC track) - (TOF pixel) matching
265 // =0 for tracks which are not tracks for matching, i.e.
266 // there is not hit on the TPC or Rvxt>200 cm
267 // >0 for tracks with positive matching procedure:
268 // =1 or 2 for non-identified tracks:
269 // =1, if the corresponding pixel is not fired,
270 // =2, if the corresponding pixel is also matched to the
271 // other track,
272 // =3 or 4 for identified tracks:
273 // =3, if identified with true time,
274 // =4, if identified with wrong time.
275 // <0 for tracks with negative mathing procedure:
276 // =-1, if track do not reach the pixel plate (curved in the
277 // magnetic field),
278 // =-2, if track is out of z-size of the TOF,
279 // =-3, if track is or into the RICH hole, or into the PHOS hole, or in the space between the plates,
280 // =-4, if track is into the dead space of the TOF.
281 // text - time of fligth from the matching procedure = time of the
282 // pixel corresponding to the track (ps)
283 // =0 for the tracks with matc<=1
284 // mext - mass of the track from the matching procedure
285 // =p*sqrt(900*(text/leng)**2-1), if 900*(text/leng)**2-1>=0
286 // =-p*sqrt(abs(900*(text/leng)**2-1)), if 900*(text/leng)**2-1<0
287
288 foutputntuple= new TNtuple("Ntuple","matching","event:ipart:imam:xvtx:yvtx:zvtx:pxvtx:pyvtx:pzvtx:time:leng:matc:text:mext",2000000); // buffersize set for 25 Pb-Pb events
289}
290
291//__________________________________________________________________
292Double_t TimeWithTailR(Double_t* x, Double_t* par)
293{
294 // sigma - par[0], alpha - par[1], part - par[2]
295 // at x<part*sigma - gauss
296 // at x>part*sigma - TMath::Exp(-x/alpha)
297 Float_t xx =x[0];
298 Double_t f;
299 if(xx<par[0]*par[2]) {
300 f = TMath::Exp(-xx*xx/(2*par[0]*par[0]));
301 } else {
302 f = TMath::Exp(-(xx-par[0]*par[2])/par[1]-0.5*par[2]*par[2]);
303 }
304 return f;
305}
306
307//____________________________________________________________________________
308void AliTOFReconstructioner::Exec(const char* datafile, Option_t *option)
309{
310 //
311 // Performs reconstruction for TOF detector
312 //
313 gBenchmark->Start("TOFReconstruction");
314
315 TFile *file = TFile::Open(datafile);
316
317 // Get AliRun object from file or create it if not on file
318 gAlice = (AliRun*)file->Get("gAlice");
319
320 AliTOF* TOF = (AliTOF *) gAlice->GetDetector ("TOF");
321 AliDetector* TPC = gAlice->GetDetector("TPC");
322
323 if (!TOF) {
324 Error("AliTOFReconstructioner","TOF not found");
325 return;
326 }
327 if (!TPC) {
328 Error("AliTOFReconstructioner","TPC Detector not found");
329 return;
330 }
331
332 if (fEdgeTails) ftail = new TF1("tail",TimeWithTailR,-2,2,3);
333
334 if (fNevents == 0) fNevents = (Int_t) gAlice->TreeE()->GetEntries();
335 // You have to set the number of event with the ad hoc setter
336 // see testrecon.C
337
338 for (Int_t ievent = 0; ievent < fNevents; ievent++) { // start loop on events
339
340 Int_t nparticles=gAlice->GetEvent(ievent);
341 if (nparticles <= 0) return;
342
343 TClonesArray* tofhits=0;
344 TClonesArray* tpchits=0;
345
346 if (TOF) tofhits = TOF->Hits();
347 if (TPC) tpchits = TPC->Hits();
348
349 TTree *TH = gAlice->TreeH();
350 if (!TH) return;
351 Int_t ntracks = (Int_t) (TH->GetEntries()); // primary tracks
352 cout << "number of primary tracked tracks in current event " << ntracks << endl; // number of primary tracked tracks
353 // array declaration and initialization
354 // TOF arrays
355 // Int_t mapPixels[AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates][AliTOFConstants::fgkNStripC][AliTOFConstants::fgkNpadZ*AliTOFConstants::fgkNpadX];
356
357 Int_t *** mapPixels = new Int_t**[AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates];
358 for (Int_t i=0; i<AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates; i++) mapPixels[i] = new Int_t*[AliTOFConstants::fgkNStripC];
359 for (Int_t i=0; i<AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates; i++) {
360 for (Int_t j=0; j<AliTOFConstants::fgkNStripC; j++) {
361 mapPixels[i][j]= new Int_t[AliTOFConstants::fgkNpadZ*AliTOFConstants::fgkNpadX];
362 }
363 }
364
365
366 // initializing the previous array
367 for (Int_t i=0;i<AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates;i++) {
368 for (Int_t j=0;j<AliTOFConstants::fgkNStripC;j++) {
369 for (Int_t l=0;l<AliTOFConstants::fgkNpadZ*AliTOFConstants::fgkNpadX;l++) {
370 mapPixels[i][j][l]=0;
371 }
372 }
373 }
374
a020d84f 375 Float_t * toftime = new Float_t[fMaxAllTracks];
f9a28264 376 InitArray(toftime, fMaxAllTracks);
db9ba97f 377 AliTOFPad* pixelArray = new AliTOFPad[fMaxPixels];
f9a28264 378 Int_t* iTOFpixel = new Int_t[fMaxAllTracks];
379 InitArray(iTOFpixel , fMaxAllTracks);
380 Int_t* kTOFhitFirst = new Int_t[fMaxAllTracks];
381 InitArray(kTOFhitFirst, fMaxAllTracks);
db9ba97f 382 AliTOFRecHit* hitArray = new AliTOFRecHit[fMaxTOFHits];
383 Int_t isHitOnFiredPad=0; // index used to fill hitArray (array used to store informations
384 // about pads that contains an hit)
385 Int_t ntotFiredPads=0; // index used to fill array -> total number of fired pads (at least one time)
386
387 // TPC arrays
388 AliTOFTrack* trackArray = new AliTOFTrack[fMaxTracks];
f9a28264 389 Int_t * iparticle = new Int_t[fMaxAllTracks];
390 InitArray(iparticle,fMaxAllTracks);
391 Int_t * iTrackPt = new Int_t[fMaxTracks];
392 InitArray(iTrackPt, fMaxTracks); // array
393 Float_t * ptTrack = new Float_t[fMaxTracks];
394 InitArray( ptTrack, fMaxTracks); // array for selected track pt
db9ba97f 395 Int_t ntotTPCtracks=0; // total number of selected TPC tracks
396
397
398 // reading TOF hits
399 if(TOF) ReadTOFHits(ntracks, TH, tofhits, mapPixels, kTOFhitFirst, pixelArray, iTOFpixel, toftime, hitArray,isHitOnFiredPad,ntotFiredPads);
400 cout << "isHitOnFiredPad " << isHitOnFiredPad << " for event " << ievent << endl;
401
402 // start debug for adding noise
403 // adding noise
404 Int_t nHitsNoNoise=isHitOnFiredPad;
405
406
407 if(fNoise) AddNoiseFromOuter(option,mapPixels,pixelArray,hitArray,isHitOnFiredPad,ntotFiredPads);
408 cout << "ntotFiredPads after adding noise " << ntotFiredPads << " for event " << ievent << endl;
409 // set the hitArray distance to nearest hit
410 SetMinDistance(hitArray,nHitsNoNoise);
411
412 // these lines has to be commented till TPC will provide fPx fPy fPz
413 // and fL in AliTPChit class
414 // reading TPC hits
415 /*
416 if(TPC) ReadTPCHits(ntracks, TH, tpchits, iTrackPt, iparticle, ptTrack, trackArray,ntotTPCtracks);
417 */
418
419 // geometrical matching
420 if(TOF && TPC) Matching(trackArray,hitArray,mapPixels,pixelArray,kTOFhitFirst,ntotFiredPads,iTrackPt,iTOFpixel,ntotTPCtracks);
421
422 // fill ntuple with reconstructed particles from current event
423 FillNtuple(ntracks,trackArray,hitArray,pixelArray,iTOFpixel,iparticle,toftime,ntotFiredPads,ntotTPCtracks);
424
425
426 // free used memory
f9a28264 427 delete [] toftime;
428 delete [] pixelArray;
429 delete [] iTOFpixel;
430 delete [] kTOFhitFirst;
431 delete [] hitArray;
432 delete [] trackArray;
433 delete [] iparticle;
434 delete [] iTrackPt;
435 delete [] ptTrack;
db9ba97f 436
437 for (Int_t i=0; i<AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates; i++) {
438 for (Int_t j=0; j<AliTOFConstants::fgkNStripC; j++) {
439 delete [] mapPixels[i][j];
440 }
441 }
442 for (Int_t i=0; i<AliTOFConstants::fgkNSectors*AliTOFConstants::fgkNPlates; i++) delete [] mapPixels[i];
443
444 delete [] mapPixels;
445
446 }//event loop
447
5fff655e 448 // free used memory for ftail
449 if (ftail)
450 {
451 delete ftail;
452 ftail = 0;
453 }
db9ba97f 454
455 // writing ntuple on output file
456 foutputfile->cd();
457 //foutputntuple->Write(0,TObject::kOverwrite);
458 foutputntuple->Write();
459 foutputfile->Write();
460 foutputfile->Close();
461
462 gBenchmark->Stop("TOFReconstruction");
463 cout << "AliTOFReconstructioner:" << endl ;
464 cout << " took " << gBenchmark->GetCpuTime("TOFReconstruction") << " seconds in order to make the reconstruction for " << fNevents << " events " << endl;
465 cout << gBenchmark->GetCpuTime("TOFReconstruction")/fNevents << " seconds per event " << endl ;
466 cout << endl ;
467
468}
469
470//__________________________________________________________________
471void AliTOFReconstructioner::SetRecFile(char * file )
472{
473 //
474 // Set the file name for reconstruction output
475 //
476 if(!fRecFile.IsNull())
477 cout << "Changing destination file for TOF reconstruction from " <<(char *)fRecFile.Data() << " to " << file << endl ;
478 fRecFile=file ;
479}
480//__________________________________________________________________
481void AliTOFReconstructioner::Print(Option_t* option)const
482{
483 //
484 // Print reconstruction output file name
485 //
486 cout << "------------------- "<< GetName() << " -------------" << endl ;
487 if(fRecFile.IsNull())
488 cout << " Writing reconstructed particles to file galice.root "<< endl ;
489 else
490 cout << " Writing reconstructed particle to file " << (char*) fRecFile.Data() << endl ;
491
492}
493
494//__________________________________________________________________
495void AliTOFReconstructioner::PrintParameters()const
496{
497 //
498 // Print parameters used for reconstruction
499 //
500 cout << " ------------------- "<< GetName() << " -------------" << endl ;
501 cout << " Parameters used for TOF reconstruction " << endl ;
502 // Printing the parameters
503
504 cout << " Number of events: " << fNevents << endl;
505 cout << " Recostruction from event "<< fFirstEvent << " to event "<< fLastEvent << endl;
506 cout << " TOF geometry parameters " << endl;
507 cout << " Min. radius of the TOF (cm) "<< AliTOFConstants::fgkrmin << endl;
508 cout << " Max. radius of the TOF (cm) "<< AliTOFConstants::fgkrmax << endl;
509 cout << " Number of TOF geom. levels "<< AliTOFConstants::fgkmaxtoftree<< endl;
510 cout << " Number of TOF sectors "<< AliTOFConstants::fgkNSectors << endl;
511 cout << " Number of TOF modules "<< AliTOFConstants::fgkNPlates << endl;
512 cout << " Max. Number of strips in a module "<< AliTOFConstants::fgkNStripC << endl;
513 cout << " Number of pads per strip "<< AliTOFConstants::fgkNpadX*AliTOFConstants::fgkNpadZ << endl;
514 cout << " Number of strips in central module "<< AliTOFConstants::fgkNStripA << endl;
515 cout << " Number of strips in intermediate modules "<< AliTOFConstants::fgkNStripB << endl;
516 cout << " Number of strips in outer modules "<< AliTOFConstants::fgkNStripC << endl;
517 cout << " Number of MRPC in x strip direction "<< AliTOFConstants::fgkNpadX<< endl;
518 cout << " Size of MRPC (cm) along X "<< AliTOFConstants::fgkXPad<< endl;
519 cout << " Number of MRPC in z strip direction "<< AliTOFConstants::fgkNpadZ<<endl;
520 cout << " Size of MRPC (cm) along Z "<< AliTOFConstants::fgkZPad<<endl;
521 cout << " Module Lengths (cm)" << endl;
522 cout << " A Module: "<< AliTOFConstants::fgkzlenA<< " B Modules: "<< AliTOFConstants::fgkzlenB<< " C Modules: "<< AliTOFConstants::fgkzlenC<< endl;
523 cout << " Inner radius of the TOF detector (cm): "<<AliTOFConstants::fgkrmin << endl;
524 cout << " Outer radius of the TOF detector (cm): "<<AliTOFConstants::fgkrmax << endl;
525 cout << " Max. half z-size of TOF (cm) : "<<AliTOFConstants::fgkMaxhZtof << endl;
526 cout << " TOF Pad parameters " << endl;
527 cout << " Time Resolution (ns) "<< fTimeResolution <<" Pad Efficiency: "<< fpadefficiency << endl;
528 cout << " Edge Effect option: "<< fEdgeEffect<< endl;
529
530 cout << " Boundary Effect Simulation Parameters " << endl;
531 cout << " Hparameter: "<< fHparameter<<" H2parameter:"<< fH2parameter <<" Kparameter:"<< fKparameter<<" K2parameter: "<< fK2parameter << endl;
532 cout << " Efficiency in the central region of the pad: "<< fEffCenter << endl;
533 cout << " Efficiency at the boundary region of the pad: "<< fEffBoundary << endl;
534 cout << " Efficiency value at H2parameter "<< fEff2Boundary << endl;
535 cout << " Efficiency value at K2parameter "<< fEff3Boundary << endl;
536 cout << " Resolution (ps) in the central region of the pad: "<< fResCenter << endl;
537 cout << " Resolution (ps) at the boundary of the pad : "<< fResBoundary << endl;
538 cout << " Slope (ps/K) for neighbouring pad : "<< fResSlope <<endl;
539 cout << " Time walk (ps) in the central region of the pad : "<< fTimeWalkCenter << endl;
540 cout << " Time walk (ps) at the boundary of the pad : "<< fTimeWalkBoundary<< endl;
541 cout << " Slope (ps/K) for neighbouring pad : "<< fTimeWalkSlope<<endl;
542 cout << " Pulse Heigth Simulation Parameters " << endl;
543 cout << " Flag for delay due to the PulseHeightEffect: "<< fTimeDelayFlag <<endl;
544 cout << " Pulse Height Slope : "<< fPulseHeightSlope<<endl;
545 cout << " Time Delay Slope : "<< fTimeDelaySlope<<endl;
546 cout << " Minimum charge amount which could be induced : "<< fMinimumCharge<<endl;
547 cout << " Smearing in charge in (q1/q2) vs x plot : "<< fChargeSmearing<<endl;
548 cout << " Smearing in log of charge ratio : "<< fLogChargeSmearing<<endl;
549 cout << " Smearing in time in time vs log(q1/q2) plot : "<< fTimeSmearing<<endl;
550 cout << " Flag for average time : "<< fAverageTimeFlag<<endl;
551 cout << " Charge factor flag for matching : "<< fChargeFactorForMatching<<endl;
552 cout << " Edge tails option : "<< fEdgeTails << endl;
553 cout << " TPC tracking parameters " << endl;
554 cout << " TPC tracking efficiency : "<< fTrackingEfficiency<< endl;
555 cout << " Sigma vs momentum dependency flag : "<< fSigmavsp << endl;
556 cout << " Space uncertainties (cm). sigma(z) (cm): "<< fSigmaZ << " sigma(R(phi)) (cm): "<< fSigmarphi << endl;
557 cout << " Momentum uncertainties. sigma(delta(P)/P): "<< fSigmap <<" sigma(phi) (rad): "<< fSigmaPhi <<" sigma(theta) (rad): "<< fSigmaTheta << endl;
558 cout << " Parameters for additional noise hits " << endl;
559 cout << " Number of noise hits : " << fNoise <<" Slope parameter (ns) in the time distribution: " << fNoiseSlope << endl;
560 cout << " Mean TOF for noise from outer regions (ns)" << fNoiseMeanTof << endl;
561 cout << " Physical parameters " << endl;
562 cout << " Magnetic Field (tesla) : "<< fField <<endl;
563 cout << " Radiation length of the outer wall of TPC: "<< fRadLenTPC << endl;
564 cout << " (TPC tracks)-(TOF pads) matching parameters " << endl;
565 cout << " TRD Correction flag : "<< fCorrectionTRD <<endl;
566 cout << " Number of the last TPC row: "<< fLastTPCRow <<" Vertex radius (cm) for selected tracks: "<<fRadiusvtxBound<<endl;
567 cout << " Max. number of test tracks: "<<fMaxTestTracks << endl;
568 cout << " Space step (cm) : "<< fStep <<endl;
569 cout << " Matching style option : "<< fMatchingStyle <<endl;
570 cout << " Array parameters " << endl;
571 cout << " Max.number of pads involved in the matching procedure: "<< fMaxPixels << endl;
572 cout << " Max.number of TOF hits per event : "<< fMaxTOFHits<< endl;
573 cout << " Max.number of tracks selected for matching : "<< fMaxTracks << endl;
574 cout << " Max.number of all tracks including the neutral ones : "<< fMaxAllTracks<< endl;
575 cout << " Debug Flag : "<< fdbg << endl;
576 cout << " Cut on momentum for selecting tracks : "<< fPBound << endl;
577
578}
579
580//__________________________________________________________________
b9d0a01d 581void AliTOFReconstructioner::IsInsideThePad(TVirtualMC *vmc, Float_t x, Float_t y, Float_t z, Int_t *nGeom, Float_t& zPad, Float_t& xPad)
db9ba97f 582{
583 // input: x,y,z - coordinates of a hit
584 // output: array nGeom[]
585 // nGeom[0] - the TOF sector number, 1,2,...,18 along azimuthal direction starting from -90 deg.!!!
586 // nGeom[1] - the TOF module number, 1,2,3,4,5=C,B,A,B,C along z-direction
587 // nGeom[2] - the TOF strip number, 1,2,... along z-direction
588 // nGeom[3] - the TOF padz number, 1,2=NPZ across a strip
589 // nGeom[4] - the TOF padx number, 1,2,...,48=NPX along a strip
590 // zPad, xPad - coordinates of the hit in the pad frame
591 // numbering is adopted for the version 3.05 of AliRoot
592 // example:
593 // from Hits: sec,pla,str,padz,padx=4,2,14,2,35
594 // Vol. n.0: ALIC, copy number 1
595 // Vol. n.1: B077, copy number 1
596 // Vol. n.2: B074, copy number 5
597 // Vol. n.3: BTO2, copy number 1
598 // Vol. n.4: FTOB, copy number 2
599 // Vol. n.5: FLTB, copy number 0
600 // Vol. n.6: FSTR, copy number 14
601 // Vol. n.7: FSEN, copy number 0
602 // Vol. n.8: FSEZ, copy number 2
603 // Vol. n.9: FSEX, copy number 35
604 // Vol. n.10: FPAD, copy number 0
605
606
607 Float_t xTOF[3];
608 Int_t sector=0,module=0,strip=0,padz=0,padx=0;
609 Int_t i,numed,nLevel,copyNumber;
610 Gcvolu_t* gcvolu;
611 char name[5];
612 name[4]=0;
613
614 for (i=0; i<AliTOFConstants::fgkmaxtoftree; i++) nGeom[i]=0;
615 zPad=100.;
616 xPad=100.;
617
618 xTOF[0]=x;
619 xTOF[1]=y;
620 xTOF[2]=z;
621
b9d0a01d 622 TGeant3 * g3 = (TGeant3*) vmc;
623
db9ba97f 624 g3->Gmedia(xTOF, numed);
625 gcvolu=g3->Gcvolu();
626 nLevel=gcvolu->nlevel;
627 if(fdbg) {
628 for (Int_t i=0; i<nLevel; i++) {
629 strncpy(name,(char*) (&gcvolu->names[i]),4);
630 cout<<"Vol. n."<<i<<": "<<name<<", copy number "<<gcvolu->number[i]<<endl;
631 }
632 }
633 if(nLevel>=2) {
634 // sector type name: B071(1,2,...,10),B074(1,2,3,4,5-PHOS),B075(1,2,3-RICH)
635 strncpy(name,(char*) (&gcvolu->names[2]),4);
636 // volume copy: 1,2,...,10 for B071, 1,2,3,4,5 for B074, 1,2,3 for B075
637 copyNumber=gcvolu->number[2];
638 if(!strcmp(name,"B071")) {
639 if (copyNumber>=6 && copyNumber<=8) {
640 sector=copyNumber+10;
641 } else if (copyNumber>=1 && copyNumber<=5){
642 sector=copyNumber+7;
643 } else {
644 sector=copyNumber-8;
645 }
646 } else if(!strcmp(name,"B075")) {
647 sector=copyNumber+12;
648 } else if(!strcmp(name,"B074")) {
649 if (copyNumber>=1 && copyNumber<=3){
650 sector=copyNumber+4;
651 } else {
652 sector=copyNumber-1;
653 }
654 }
655 }
656 if(sector) {
657 nGeom[0]=sector;
658 if(nLevel>=4) {
659 // we'll use the module value in z-direction:
660 // 1 2 3 4 5
661 // the module order in z-direction: FTOC,FTOB,FTOA,FTOB,FTOC
662 // the module copy: 2 2 0 1 1
663 // module type name: FTOA, FTOB, FTOC
664 strncpy(name,(char*) (&gcvolu->names[4]),4);
665 // module copy:
666 copyNumber=gcvolu->number[4];
667 if(!strcmp(name,"FTOC")) {
668 if (copyNumber==2) {
669 module=1;
670 } else {
671 module=5;
672 }
673 } else if(!strcmp(name,"FTOB")) {
674 if (copyNumber==2) {
675 module=2;
676 } else {
677 module=4;
678 }
679 } else if(!strcmp(name,"FTOA")) {
680 module=3;
681 }
682 }
683 }
684
685 if(module) {
686 nGeom[1]=module;
687 if(nLevel>=6) {
688 // strip type name: FSTR
689 strncpy(name,(char*) (&gcvolu->names[6]),4);
690 // strip copy:
691 copyNumber=gcvolu->number[6];
692 if(!strcmp(name,"FSTR")) strip=copyNumber;
693 }
694 }
695
696 if(strip) {
697 nGeom[2]=strip;
698 if(nLevel>=8) {
699 // padz type name: FSEZ
700 strncpy(name,(char*) (&gcvolu->names[8]),4);
701 // padz copy:
702 copyNumber=gcvolu->number[8];
703 if(!strcmp(name,"FSEZ")) padz=copyNumber;
704 }
705 }
706 if(padz) {
707 nGeom[3]=padz;
708 if(nLevel>=9) {
709 // padx type name: FSEX
710 strncpy(name,(char*) (&gcvolu->names[9]),4);
711 // padx copy:
712 copyNumber=gcvolu->number[9];
713 if(!strcmp(name,"FSEX")) padx=copyNumber;
714 }
715 }
716
717 if(padx) {
718 nGeom[4]=padx;
719 zPad=gcvolu->glx[2]; // check here
720 xPad=gcvolu->glx[0]; // check here
721 }
722
723 // printf(" nGeom[0,1,2,3,4]=%i,%i,%i,%i,%i\n",nGeom[0],nGeom[1],nGeom[2],nGeom[3],nGeom[4]);
724}
725
726//__________________________________________________________________
727void AliTOFReconstructioner::EpMulScatt(Float_t& px, Float_t& py, Float_t& pz, Float_t& p, Float_t& theta)
728{
729 // Momentum p - before mult.scat.
730 // Momentum p2 - after mult.scat.
731 // THE0 - r.m.s. of deviation angle in plane
732 // (see RPP'96: Phys.Rev.D54 (1996) 134)
733
734 Float_t pt,thex,they,tantx,tanty,p2px,p2py,p2pz,costhe,sinthe,cospsi,sinpsi,p2x,p2y,p2z,p2,g;
735
736 pt=TMath::Sqrt(px*px+py*py);
737 // angles for p in the ' frame with Z'along p
738 if(fMatchingStyle==1) {
739 thex=theta*gRandom->Gaus();
740 they=theta*gRandom->Gaus();
741 } else {
742 thex=3*(-theta+2*theta*gRandom->Rndm());
743 they=3*(-theta+2*theta*gRandom->Rndm());
744 }
745 tantx=TMath::Tan(thex);
746 tanty=TMath::Tan(they);
747
748 // p2p - p2 in the ' frame
749 p2pz=p/TMath::Sqrt(1.+tantx*tantx+tanty*tanty);
750 p2py=p2pz*tanty;
751 p2px=p2pz*tantx;
752 // choose X'so that PHI=0 (see Il'in, Pozdnyak Analiticheskaya geometriya, 1968, c.88
753 // for Euler angles PSI, THETA (PHI=0)
754 costhe=pz/p;
755 sinthe=pt/p;
756 cospsi=-py/pt;
757 sinpsi=px/pt;
758 //
759 g=p2py*costhe-p2pz*sinthe;
760 p2x=p2px*cospsi-g*sinpsi;
761 p2y=p2px*sinpsi+g*cospsi;
762 p2z=p2py*sinthe+p2pz*costhe;
763 p2=TMath::Sqrt(p2x*p2x+p2y*p2y+p2z*p2z);
764
765 // Test angle
766 g=(px*p2x+py*p2y+pz*p2z)/(p*p2);
767 if(g>1) g=1;
768 theta=TMath::ACos(g);
769 px=p2x;
770 py=p2y;
771 pz=p2z;
772 p=p2;
773
774}
775
776// std border effect algorithm
777//__________________________________________________________________
778void AliTOFReconstructioner::BorderEffect(Float_t z0, Float_t x0, Float_t geantTime, Int_t& nActivatedPads, Int_t& nFiredPads, Bool_t* isFired, Int_t* nPlace, Float_t* qInduced, Float_t* tofTime, Float_t& averageTime)
779{
780 // Input: z0, x0 - hit position in the strip system (0,0 - center of the strip), cm
781 // geantTime - time generated by Geant, ns
782 // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4)
783 // nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads)
784 // qInduced[iPad]- charge induced on pad, arb. units
785 // this array is initialized at zero by the caller
786 // tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns
787 // this array is initialized at zero by the caller
788 // averageTime - time given by pad hited by the Geant track taking into account the times (weighted) given by the pads fired for edge effect also.
789 // The weight is given by the qInduced[iPad]/qCenterPad
790 // this variable is initialized at zero by the caller
791 // nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3
792 // this variable is initialized at zero by the caller
793 //
794 // Description of used variables:
795 // eff[iPad] - efficiency of the pad
796 // res[iPad] - resolution of the pad, ns
797 // timeWalk[iPad] - time walk of the pad, ns
798 // timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns
799 // PadId[iPad] - Pad Identifier
800 // E | F --> PadId[iPad] = 5 | 6
801 // A | B --> PadId[iPad] = 1 | 2
802 // C | D --> PadId[iPad] = 3 | 4
803 // nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB
804 // qCenterPad - charge extimated for each pad, arb. units
805 // weightsSum - sum of weights extimated for each pad fired, arb. units
806
807 const Float_t kSigmaForTail[2] = {AliTOFConstants::fgkSigmaForTail1,AliTOFConstants::fgkSigmaForTail2}; //for tail
808 Int_t iz = 0, ix = 0;
809 Float_t dX = 0., dZ = 0., x = 0., z = 0.;
810 Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter;
811 Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.;
812 Float_t logOfqInd = 0.;
813 Float_t weightsSum = 0.;
814 Int_t nTail[4] = {0,0,0,0};
815 Int_t padId[4] = {0,0,0,0};
816 Float_t eff[4] = {0.,0.,0.,0.};
817 Float_t res[4] = {0.,0.,0.,0.};
818 // Float_t qCenterPad = fMinimumCharge * fMinimumCharge;
819 Float_t qCenterPad = 1.;
820 Float_t timeWalk[4] = {0.,0.,0.,0.};
821 Float_t timeDelay[4] = {0.,0.,0.,0.};
822
823 nActivatedPads = 0;
824 nFiredPads = 0;
825
826 (z0 <= 0) ? iz = 0 : iz = 1;
827 dZ = z0 + (0.5 * AliTOFConstants::fgkNpadZ - iz - 0.5) * AliTOFConstants::fgkZPad; // hit position in the pad frame, (0,0) - center of the pad
828 z = 0.5 * AliTOFConstants::fgkZPad - TMath::Abs(dZ); // variable for eff., res. and timeWalk. functions
829 iz++; // z row: 1, ..., AliTOFConstants::fgkNpadZ = 2
830 ix = (Int_t)((x0 + 0.5 * AliTOFConstants::fgkNpadX * AliTOFConstants::fgkXPad) / AliTOFConstants::fgkXPad);
831 dX = x0 + (0.5 * AliTOFConstants::fgkNpadX - ix - 0.5) * AliTOFConstants::fgkXPad; // hit position in the pad frame, (0,0) - center of the pad
832 x = 0.5 * AliTOFConstants::fgkXPad - TMath::Abs(dX); // variable for eff., res. and timeWalk. functions;
833 ix++; // x row: 1, ..., AliTOFConstants::fgkNpadX = 48
834
835 ////// Pad A:
836 nActivatedPads++;
837 nPlace[nActivatedPads-1] = (iz - 1) * AliTOFConstants::fgkNpadX + ix;
838 qInduced[nActivatedPads-1] = qCenterPad;
839 padId[nActivatedPads-1] = 1;
840
841 if (fEdgeEffect == 0) {
842 eff[nActivatedPads-1] = fEffCenter;
843 if (gRandom->Rndm() < eff[nActivatedPads-1]) {
844 nFiredPads = 1;
845 res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + fResCenter * fResCenter); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns;
846 isFired[nActivatedPads-1] = kTRUE;
847 tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]);
848 averageTime = tofTime[nActivatedPads-1];
849 }
850 } else {
851
852 if(z < h) {
853 if(z < h2) {
854 effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2;
855 } else {
856 effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2);
857 }
858 resZ = fResBoundary + (fResCenter - fResBoundary) * z / h;
859 timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h;
860 nTail[nActivatedPads-1] = 1;
861 } else {
862 effZ = fEffCenter;
863 resZ = fResCenter;
864 timeWalkZ = fTimeWalkCenter;
865 }
866
867 if(x < h) {
868 if(x < h2) {
869 effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2;
870 } else {
871 effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2);
872 }
873 resX = fResBoundary + (fResCenter - fResBoundary) * x / h;
874 timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h;
875 nTail[nActivatedPads-1] = 1;
876 } else {
877 effX = fEffCenter;
878 resX = fResCenter;
879 timeWalkX = fTimeWalkCenter;
880 }
881
882 (effZ<effX) ? eff[nActivatedPads-1] = effZ : eff[nActivatedPads-1] = effX;
883 (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
884 (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
885
886
887 ////// Pad B:
888 if(z < k2) {
889 effZ = fEffBoundary - (fEffBoundary - fEff3Boundary) * (z / k2);
890 } else {
891 effZ = fEff3Boundary * (k - z) / (k - k2);
892 }
893 resZ = fResBoundary + fResSlope * z / k;
894 timeWalkZ = fTimeWalkBoundary + fTimeWalkSlope * z / k;
895
896 if(z < k && z > 0) {
897 if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
898 nActivatedPads++;
899 nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX;
900 eff[nActivatedPads-1] = effZ;
901 res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
902 timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns
903 nTail[nActivatedPads-1] = 2;
904 if (fTimeDelayFlag) {
905 // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
906 // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
907 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
908 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
909 timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
910 } else {
911 timeDelay[nActivatedPads-1] = 0.;
912 }
913 padId[nActivatedPads-1] = 2;
914 }
915 }
916
917
918 ////// Pad C, D, E, F:
919 if(x < k2) {
920 effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2);
921 } else {
922 effX = fEff3Boundary * (k - x) / (k - k2);
923 }
924 resX = fResBoundary + fResSlope*x/k;
925 timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k;
926
927 if(x < k && x > 0) {
928 // C:
929 if(ix > 1 && dX < 0) {
930 nActivatedPads++;
931 nPlace[nActivatedPads-1] = nPlace[0] - 1;
932 eff[nActivatedPads-1] = effX;
933 res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
934 timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
935 nTail[nActivatedPads-1] = 2;
936 if (fTimeDelayFlag) {
937 // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
938 // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
939 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
940 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
941 timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
942 } else {
943 timeDelay[nActivatedPads-1] = 0.;
944 }
945 padId[nActivatedPads-1] = 3;
946
947 // D:
948 if(z < k && z > 0) {
949 if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
950 nActivatedPads++;
951 nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX - 1;
952 eff[nActivatedPads-1] = effX * effZ;
953 (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
954 (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
955
956 nTail[nActivatedPads-1] = 2;
957 if (fTimeDelayFlag) {
958 if (TMath::Abs(x) < TMath::Abs(z)) {
959 // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
960 // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
961 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
962 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
963 } else {
964 // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
965 // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
966 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
967 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
968 }
969 timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
970 } else {
971 timeDelay[nActivatedPads-1] = 0.;
972 }
973 padId[nActivatedPads-1] = 4;
974 }
975 } // end D
976 } // end C
977
978 // E:
979 if(ix < AliTOFConstants::fgkNpadX && dX > 0) {
980 nActivatedPads++;
981 nPlace[nActivatedPads-1] = nPlace[0] + 1;
982 eff[nActivatedPads-1] = effX;
983 res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(10400 + resX * resX)); // ns
984 timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
985 nTail[nActivatedPads-1] = 2;
986 if (fTimeDelayFlag) {
987 // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
988 // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
989 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
990 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
991 timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
992 } else {
993 timeDelay[nActivatedPads-1] = 0.;
994 }
995 padId[nActivatedPads-1] = 5;
996
997
998 // F:
999 if(z < k && z > 0) {
1000 if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1001 nActivatedPads++;
1002 nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX + 1;
1003 eff[nActivatedPads - 1] = effX * effZ;
1004 (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1005 (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001*timeWalkX; // ns
1006 nTail[nActivatedPads-1] = 2;
1007 if (fTimeDelayFlag) {
1008 if (TMath::Abs(x) < TMath::Abs(z)) {
1009 // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1010 // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1011 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * z);
1012 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1013 } else {
1014 // qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1015 // qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1016 qInduced[nActivatedPads-1] = TMath::Exp(-fPulseHeightSlope * x);
1017 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1018 }
1019 timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1020 } else {
1021 timeDelay[nActivatedPads-1] = 0.;
1022 }
1023 padId[nActivatedPads-1] = 6;
1024 }
1025 } // end F
1026 } // end E
1027 } // end if(x < k)
1028
1029
1030 for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) {
1031 if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution;
1032 if(gRandom->Rndm() < eff[iPad]) {
1033 isFired[iPad] = kTRUE;
1034 nFiredPads++;
1035 if(fEdgeTails) {
1036 if(nTail[iPad] == 0) {
1037 tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1038 } else {
1039 ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]);
1040 Double_t timeAB = ftail->GetRandom();
1041 tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB;
1042 }
1043 } else {
1044 tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1045 }
1046 if (fAverageTimeFlag) {
1047 averageTime += tofTime[iPad] * qInduced[iPad];
1048 weightsSum += qInduced[iPad];
1049 } else {
1050 averageTime += tofTime[iPad];
1051 weightsSum += 1.;
1052 }
1053 }
1054 }
1055 if (weightsSum!=0) averageTime /= weightsSum;
1056 } // end else (fEdgeEffect != 0)
1057}
1058
1059
1060/* new algorithm (to be checked)
1061//__________________________________________________________________
1062void AliTOFReconstructioner::BorderEffect(Float_t z0, Float_t x0, Float_t geantTime, Int_t& nActivatedPads, Int_t& nFiredPads, Bool_t* isFired, Int_t* nPlace, Float_t* qInduced, Float_t* tofTime, Float_t& averageTime)
1063{
1064 // Input: z0, x0 - hit position in the strip system (0,0 - center of the strip), cm
1065 // geantTime - time generated by Geant, ns
1066 // Output: nActivatedPads - the number of pads activated by the hit (1 || 2 || 4)
1067 // nFiredPads - the number of pads fired (really activated) by the hit (nFiredPads <= nActivatedPads)
1068 // qInduced[iPad]- charge induced on pad, arb. units
1069 // this array is initialized at zero by the caller
1070 // tofAfterSimul[iPad] - time calculated with edge effect algorithm, ns
1071 // this array is initialized at zero by the caller
1072 // averageTime - time given by pad hited by the Geant track taking into account the times (weighted) given by the pads fired for edge effect also.
1073 // The weight is given by the qInduced[iPad]/qCenterPad
1074 // this variable is initialized at zero by the caller
1075 // nPlace[iPad] - the number of the pad place, iPad = 0, 1, 2, 3
1076 // this variable is initialized at zero by the caller
1077 //
1078 // Description of used variables:
1079 // eff[iPad] - efficiency of the pad
1080 // res[iPad] - resolution of the pad, ns
1081 // timeWalk[iPad] - time walk of the pad, ns
1082 // timeDelay[iPad] - time delay for neighbouring pad to hited pad, ns
1083 // PadId[iPad] - Pad Identifier
1084 // E | F --> PadId[iPad] = 5 | 6
1085 // A | B --> PadId[iPad] = 1 | 2
1086 // C | D --> PadId[iPad] = 3 | 4
1087 // nTail[iPad] - the tail number, = 1 for tailA, = 2 for tailB
1088 // qCenterPad - charge extimated for each pad, arb. units
1089 // weightsSum - sum of weights extimated for each pad fired, arb. units
1090
1091 const Float_t kSigmaForTail[2] = {AliTOFConstants::fgkSigmaForTail1,AliTOFConstants::fgkSigmaForTail2}; //for tail
1092 Int_t iz = 0, ix = 0;
1093 Float_t dX = 0., dZ = 0., x = 0., z = 0.;
1094 Float_t h = fHparameter, h2 = fH2parameter, k = fKparameter, k2 = fK2parameter;
1095 Float_t effX = 0., effZ = 0., resX = 0., resZ = 0., timeWalkX = 0., timeWalkZ = 0.;
1096 Float_t logOfqInd = 0.;
1097 Float_t weightsSum = 0.;
1098 Int_t nTail[4] = {0,0,0,0};
1099 Int_t padId[4] = {0,0,0,0};
1100 Float_t eff[4] = {0.,0.,0.,0.};
1101 Float_t res[4] = {0.,0.,0.,0.};
1102 Float_t qCenterPad = fMinimumCharge * fMinimumCharge;
1103 Float_t timeWalk[4] = {0.,0.,0.,0.};
1104 Float_t timeDelay[4] = {0.,0.,0.,0.};
1105
1106 nActivatedPads = 0;
1107 nFiredPads = 0;
1108
1109 (z0 <= 0) ? iz = 0 : iz = 1;
1110 dZ = z0 + (0.5 * AliTOFConstants::fgkNpadZ - iz - 0.5) * AliTOFConstants::fgkZPad; // hit position in the pad frame, (0,0) - center of the pad
1111 z = 0.5 * AliTOFConstants::fgkZPad - TMath::Abs(dZ); // variable for eff., res. and timeWalk. functions
1112 iz++; // z row: 1, ..., AliTOFConstants::fgkNpadZ = 2
1113 ix = (Int_t)((x0 + 0.5 * AliTOFConstants::fgkNpadX * AliTOFConstants::fgkXPad) / AliTOFConstants::fgkXPad);
1114 dX = x0 + (0.5 * AliTOFConstants::fgkNpadX - ix - 0.5) * AliTOFConstants::fgkXPad; // hit position in the pad frame, (0,0) - center of the pad
1115 x = 0.5 * AliTOFConstants::fgkXPad - TMath::Abs(dX); // variable for eff., res. and timeWalk. functions;
1116 ix++; // x row: 1, ..., AliTOFConstants::fgkNpadX = 48
1117
1118 ////// Pad A:
1119 nActivatedPads++;
1120 nPlace[nActivatedPads-1] = (iz - 1) * AliTOFConstants::fgkNpadX + ix;
1121 qInduced[nActivatedPads-1] = qCenterPad;
1122 padId[nActivatedPads-1] = 1;
1123
1124 if (fEdgeEffect == 0) {
1125 eff[nActivatedPads-1] = fEffCenter;
1126 if (gRandom->Rndm() < eff[nActivatedPads-1]) {
1127 nFiredPads = 1;
1128 res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + fResCenter * fResCenter); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns;
1129 isFired[nActivatedPads-1] = kTRUE;
1130 tofTime[nActivatedPads-1] = gRandom->Gaus(geantTime + fTimeWalkCenter, res[0]);
1131 averageTime = tofTime[nActivatedPads-1];
1132 }
1133 } else {
1134
1135 if(z < h) {
1136 if(z < h2) {
1137 effZ = fEffBoundary + (fEff2Boundary - fEffBoundary) * z / h2;
1138 } else {
1139 effZ = fEff2Boundary + (fEffCenter - fEff2Boundary) * (z - h2) / (h - h2);
1140 }
1141 resZ = fResBoundary + (fResCenter - fResBoundary) * z / h;
1142 timeWalkZ = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * z / h;
1143 nTail[nActivatedPads-1] = 1;
1144 } else {
1145 effZ = fEffCenter;
1146 resZ = fResCenter;
1147 timeWalkZ = fTimeWalkCenter;
1148 }
1149
1150 if(x < h) {
1151 if(x < h2) {
1152 effX = fEffBoundary + (fEff2Boundary - fEffBoundary) * x / h2;
1153 } else {
1154 effX = fEff2Boundary + (fEffCenter - fEff2Boundary) * (x - h2) / (h - h2);
1155 }
1156 resX = fResBoundary + (fResCenter - fResBoundary) * x / h;
1157 timeWalkX = fTimeWalkBoundary + (fTimeWalkCenter - fTimeWalkBoundary) * x / h;
1158 nTail[nActivatedPads-1] = 1;
1159 } else {
1160 effX = fEffCenter;
1161 resX = fResCenter;
1162 timeWalkX = fTimeWalkCenter;
1163 }
1164
1165 (effZ<effX) ? eff[nActivatedPads-1] = effZ : eff[nActivatedPads-1] = effX;
1166 (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1167 (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1168
1169
1170 ////// Pad B:
1171 if(z < k2) {
1172 effZ = fEffBoundary - (fEffBoundary - fEff3Boundary) * (z / k2);
1173 } else {
1174 effZ = fEff3Boundary * (k - z) / (k - k2);
1175 }
1176 resZ = fResBoundary + fResSlope * z / k;
1177 timeWalkZ = fTimeWalkBoundary + fTimeWalkSlope * z / k;
1178
1179 if(z < k && z > 0) {
1180 if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1181 nActivatedPads++;
1182 nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX;
1183 eff[nActivatedPads-1] = effZ;
1184 res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1185 timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ; // ns
1186 nTail[nActivatedPads-1] = 2;
1187 if (fTimeDelayFlag) {
1188 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1189 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1190 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1191 timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1192 } else {
1193 timeDelay[nActivatedPads-1] = 0.;
1194 }
1195 padId[nActivatedPads-1] = 2;
1196 }
1197 }
1198
1199
1200 ////// Pad C, D, E, F:
1201 if(x < k2) {
1202 effX = fEffBoundary - (fEffBoundary - fEff3Boundary) * (x / k2);
1203 } else {
1204 effX = fEff3Boundary * (k - x) / (k - k2);
1205 }
1206 resX = fResBoundary + fResSlope*x/k;
1207 timeWalkX = fTimeWalkBoundary + fTimeWalkSlope*x/k;
1208
1209 if(x < k && x > 0) {
1210 // C:
1211 if(ix > 1 && dX < 0) {
1212 nActivatedPads++;
1213 nPlace[nActivatedPads-1] = nPlace[0] - 1;
1214 eff[nActivatedPads-1] = effX;
1215 res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1216 timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1217 nTail[nActivatedPads-1] = 2;
1218 if (fTimeDelayFlag) {
1219 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1220 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1221 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1222 timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1223 } else {
1224 timeDelay[nActivatedPads-1] = 0.;
1225 }
1226 padId[nActivatedPads-1] = 3;
1227
1228 // D:
1229 if(z < k && z > 0) {
1230 if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1231 nActivatedPads++;
1232 nPlace[nActivatedPads-1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX - 1;
1233 eff[nActivatedPads-1] = effX * effZ;
1234 (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1235 (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1236
1237 nTail[nActivatedPads-1] = 2;
1238 if (fTimeDelayFlag) {
1239 if (TMath::Abs(x) < TMath::Abs(z)) {
1240 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1241 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1242 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1243 } else {
1244 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1245 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1246 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1247 }
1248 timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1249 } else {
1250 timeDelay[nActivatedPads-1] = 0.;
1251 }
1252 padId[nActivatedPads-1] = 4;
1253 }
1254 } // end D
1255 } // end C
1256
1257 // E:
1258 if(ix < AliTOFConstants::fgkNpadX && dX > 0) {
1259 nActivatedPads++;
1260 nPlace[nActivatedPads-1] = nPlace[0] + 1;
1261 eff[nActivatedPads-1] = effX;
1262 res[nActivatedPads-1] = 0.001 * (TMath::Sqrt(10400 + resX * resX)); // ns
1263 timeWalk[nActivatedPads-1] = 0.001 * timeWalkX; // ns
1264 nTail[nActivatedPads-1] = 2;
1265 if (fTimeDelayFlag) {
1266 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1267 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1268 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1269 timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1270 } else {
1271 timeDelay[nActivatedPads-1] = 0.;
1272 }
1273 padId[nActivatedPads-1] = 5;
1274
1275
1276 // F:
1277 if(z < k && z > 0) {
1278 if( (iz == 1 && dZ > 0) || (iz == 2 && dZ < 0) ) {
1279 nActivatedPads++;
1280 nPlace[nActivatedPads - 1] = nPlace[0] + (3 - 2 * iz) * AliTOFConstants::fgkNpadX + 1;
1281 eff[nActivatedPads - 1] = effX * effZ;
1282 (resZ<resX) ? res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resX * resX) : res[nActivatedPads-1] = 0.001 * TMath::Sqrt(10400 + resZ * resZ); // 10400=30^2+20^2+40^2+50^2+50^2+50^2 ns
1283 (timeWalkZ<timeWalkX) ? timeWalk[nActivatedPads-1] = 0.001 * timeWalkZ : timeWalk[nActivatedPads-1] = 0.001*timeWalkX; // ns
1284 nTail[nActivatedPads-1] = 2;
1285 if (fTimeDelayFlag) {
1286 if (TMath::Abs(x) < TMath::Abs(z)) {
1287 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * z / 2.);
1288 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * z / 2.);
1289 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * z, fLogChargeSmearing);
1290 } else {
1291 qInduced[0] = fMinimumCharge * TMath::Exp(fPulseHeightSlope * x / 2.);
1292 qInduced[nActivatedPads-1] = fMinimumCharge * TMath::Exp(-fPulseHeightSlope * x / 2.);
1293 logOfqInd = gRandom->Gaus(-fPulseHeightSlope * x, fLogChargeSmearing);
1294 }
1295 timeDelay[nActivatedPads-1] = gRandom->Gaus(-fTimeDelaySlope * logOfqInd, fTimeSmearing);
1296 } else {
1297 timeDelay[nActivatedPads-1] = 0.;
1298 }
1299 padId[nActivatedPads-1] = 6;
1300 }
1301 } // end F
1302 } // end E
1303 } // end if(x < k)
1304
1305
1306 for (Int_t iPad = 0; iPad < nActivatedPads; iPad++) {
1307 if (res[iPad] < fTimeResolution) res[iPad] = fTimeResolution;
1308 if(gRandom->Rndm() < eff[iPad]) {
1309 isFired[iPad] = kTRUE;
1310 nFiredPads++;
1311 if(fEdgeTails) {
1312 if(nTail[iPad] == 0) {
1313 tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1314 } else {
1315 ftail->SetParameters(res[iPad], 2. * res[iPad], kSigmaForTail[nTail[iPad]-1]);
1316 Double_t timeAB = ftail->GetRandom();
1317 tofTime[iPad] = geantTime + timeWalk[iPad] + timeDelay[iPad] + timeAB;
1318 }
1319 } else {
1320 tofTime[iPad] = gRandom->Gaus(geantTime + timeWalk[iPad] + timeDelay[iPad], res[iPad]);
1321 }
1322 if (fAverageTimeFlag) {
1323 averageTime += tofTime[iPad] * qInduced[iPad];
1324 weightsSum += qInduced[iPad];
1325 } else {
1326 averageTime += tofTime[iPad];
1327 weightsSum += 1.;
1328 }
1329 }
1330 }
1331 if (weightsSum!=0) averageTime /= weightsSum;
1332
1333 } // end else (fEdgeEffect != 0)
1334
1335 //cout << "timedelay " << timeDelay[0] << endl;
1336 //cout << "timedelay " << timeDelay[1] << endl;
1337 //cout << "timedelay " << timeDelay[2] << endl;
1338 //cout << "timedelay " << timeDelay[3] << endl;
1339
1340}
1341*/
1342
1343
1344//__________________________________________________________________
1345Int_t AliTOFReconstructioner::PDGtoGeantCode(Int_t pdgcode)
1346{
1347 //
1348 // Gives the GEANT code from KF code of LUND JETSET
1349 //
1350 Int_t geantCode=0; // default value
1351 switch (pdgcode) {
1352 case 22:
1353 geantCode=1; // GAMMA
1354 break ;
1355 case -11:
1356 geantCode=2; // E+
1357 break ;
1358 case 11:
1359 geantCode=3; // E-
1360 break ;
1361 case 12:
1362 geantCode=4; // NUE
1363 break ;
1364 case 14:
1365 geantCode=4; // NUMU
1366 break ;
1367 case -13:
1368 geantCode=5; // MU+
1369 break ;
1370 case 13:
1371 geantCode=6; // MU-
1372 break ;
1373 case 111:
1374 geantCode=7; // PI0
1375 break ;
1376 case 211:
1377 geantCode=8; // PI+
1378 break ;
1379 case -211:
1380 geantCode=9; // PI-
1381 break ;
1382 case 130:
1383 geantCode=10; // K_L0
1384 break ;
1385 case 321:
1386 geantCode=11; // K+
1387 break ;
1388 case -321:
1389 geantCode=12; // K-
1390 break ;
1391 case 2112:
1392 geantCode=13; // N0
1393 break ;
1394 case 2212:
1395 geantCode=14; // P+
1396 break ;
1397 case -2212:
1398 geantCode=15; // P~-
1399 break ;
1400 case 310:
1401 geantCode=16; // K_S0
1402 break ;
1403 case 221:
1404 geantCode=17; // ETA
1405 break ;
1406 case 3122:
1407 geantCode=18; // LAMBDA0
1408 break ;
1409 case 3222:
1410 geantCode=19; // SIGMA+
1411 break ;
1412 case 3212:
1413 geantCode=20; // SIGMA0
1414 break ;
1415 case 3112:
1416 geantCode=21; // SIGMA-
1417 break ;
1418 case 3322:
1419 geantCode=22; // XI0
1420 break ;
1421 case 3312:
1422 geantCode=23; // XI-
1423 break ;
1424 case 3334:
1425 geantCode=24; // OMEGA-
1426 break ;
1427 case -2112:
1428 geantCode=25; // N~0
1429 break ;
1430 case -3122:
1431 geantCode=26; // LAMBDA~0
1432 break ;
1433 case -3112:
1434 geantCode=27; // SIGMA~+
1435 break ;
1436 case -3212:
1437 geantCode=28; // SIGMA~0
1438 break ;
1439 case -3222:
1440 geantCode=29; // SIGMA~-
1441 break ;
1442 case -3322:
1443 geantCode=30; // XI~0
1444 break ;
1445 case -3312:
1446 geantCode=31; // XI~+
1447 break ;
1448 case -3334:
1449 geantCode=32; // OMEGA~+
1450 break ;
1451 case 223:
1452 geantCode=33; // OMEGA(782)
1453 break ;
1454 case 333:
1455 geantCode=34; // PHI(1020)
1456 break ;
1457 case 411:
1458 geantCode=35; // D+
1459 break ;
1460 case -411:
1461 geantCode=36; // D-
1462 break ;
1463 case 421:
1464 geantCode=37; // D0
1465 break ;
1466 case -421:
1467 geantCode=38; // D~0
1468 break ;
1469 case 431:
1470 geantCode=39; // D_S+
1471 break ;
1472 case -431:
1473 geantCode=40; // D_S~-
1474 break ;
1475 case 4122:
1476 geantCode=41; // LAMBDA_C+
1477 break ;
1478 case 213:
1479 geantCode=42; // RHP(770)+
1480 break ;
1481 case -213:
1482 geantCode=43; // RHO(770)-
1483 break ;
1484 case 113:
1485 geantCode=44; // RHO(770)0
1486 break ;
1487 default:
1488 geantCode=45;
1489 break;
1490 }
1491
1492 return geantCode;
1493}
1494
1495//__________________________________________________________________
1496Bool_t AliTOFReconstructioner::operator==( AliTOFReconstructioner const & tofrec)const
1497{
1498 // Equal operator.
1499 // Reconstructioners are equal if their parameters are equal
1500
1501 // split the member variables in analogous categories
1502
1503 // time resolution and edge effect parameters
be167bf4 1504 Bool_t dummy0=(fTimeResolution==tofrec.fTimeResolution)&&(fpadefficiency==tofrec.fpadefficiency)&&(fEdgeEffect==tofrec.fEdgeEffect)&&(fEdgeTails==tofrec.fEdgeTails)&&(fHparameter==tofrec.fHparameter)&&(fH2parameter==tofrec.fH2parameter)&&(fKparameter==t
1505ofrec.fKparameter)&&(fK2parameter==tofrec.fK2parameter);
db9ba97f 1506
1507 // pad efficiency parameters
be167bf4 1508 Bool_t dummy1=(fEffCenter==tofrec.fEffCenter)&&(fEffBoundary==tofrec.fEffBoundary)&&(fEff2Boundary==tofrec.fEff2Boundary)&&(fEff3Boundary==tofrec.fEff3Boundary)&&(fResCenter==tofrec.fResCenter)&&(fResBoundary==tofrec.fResBoundary)&&(fResSlope==tofrec.fR
1509esSlope);
db9ba97f 1510
1511 // time walk parameters
be167bf4 1512 Bool_t dummy2=(fTimeWalkCenter==tofrec.fTimeWalkCenter)&&(fTimeWalkBoundary==tofrec.fTimeWalkBoundary)&&(fTimeWalkSlope==tofrec.fTimeWalkSlope)&&(fTimeDelayFlag==tofrec.fTimeDelayFlag)&&(fPulseHeightSlope==tofrec.fPulseHeightSlope)&&(fTimeDelaySlope==to
1513frec.fTimeDelaySlope);
db9ba97f 1514
1515 // ADC-TDC correlation parameters
be167bf4 1516 Bool_t dummy3=(fMinimumCharge==tofrec.fMinimumCharge)&&(fChargeSmearing==tofrec.fChargeSmearing )&&(fLogChargeSmearing==tofrec.fLogChargeSmearing )&&(fTimeSmearing==tofrec.fTimeSmearing )&&(fAverageTimeFlag==tofrec.fAverageTimeFlag)&&(fChargeFactorForMa
1517tching==tofrec.fChargeFactorForMatching)&&(fMatchingStyle==tofrec.fMatchingStyle);
db9ba97f 1518
be167bf4 1519 Bool_t dummy4=(fTrackingEfficiency==tofrec.fTrackingEfficiency)&&(fSigmavsp==tofrec.fSigmavsp)&&(fSigmaZ==tofrec.fSigmaZ)&&(fSigmarphi==tofrec.fSigmarphi)&&(fSigmap==tofrec.fSigmap)&&(fSigmaPhi==tofrec.fSigmaPhi)&&(fSigmaTheta==tofrec.fSigmaTheta)&&(fNo
1520ise==tofrec.fNoise)&&(fNoiseSlope==tofrec.fNoiseSlope)&&(fField==tofrec.fField)&&(fRadLenTPC==tofrec.fRadLenTPC)&&(fCorrectionTRD==tofrec.fCorrectionTRD)&&(fLastTPCRow==tofrec.fLastTPCRow)&&(fRadiusvtxBound==tofrec.fRadiusvtxBound)&&(fMaxTestTracks==tofre
1521c.fMaxTestTracks)&&(fStep==tofrec.fStep)&&(fMaxPixels==tofrec.fMaxPixels)&&(fMaxAllTracks==tofrec.fMaxAllTracks)&&(fMaxTracks==tofrec.fMaxTracks)&&(fMaxTOFHits==tofrec.fMaxTOFHits)&&(fPBound==tofrec.fPBound);
db9ba97f 1522
1523 if( dummy0 && dummy1 && dummy2 && dummy3 && dummy4)
1524 return kTRUE ;
1525 else
1526 return kFALSE ;
1527
1528}
1529//____________________________________________________________________________
1530void AliTOFReconstructioner::UseHitsFrom(const char * filename)
1531{
1532 SetTitle(filename) ;
1533}
1534
1535//____________________________________________________________________________
1536void AliTOFReconstructioner::InitArray(Float_t array[], Int_t nlocations)
1537{
1538 //
1539 // Initialize the array of Float_t
1540 //
1541 for (Int_t i = 0; i < nlocations; i++) {
1542 array[i]=0.;
1543 } // end loop
1544
1545}
1546
1547//____________________________________________________________________________
1548void AliTOFReconstructioner::InitArray(Int_t array[], Int_t nlocations)
1549{
1550 //
1551 // Initialize the array of Int_t
1552 //
1553 for (Int_t i = 0; i < nlocations; i++) {
1554 array[i]=0;
1555 } // end loop
1556
1557}
1558
1559
1560//____________________________________________________________________________
be167bf4 1561void AliTOFReconstructioner::ReadTOFHits(Int_t ntracks, TTree* treehits, TClonesArray* tofhits, Int_t ***MapPixels, Int_t* kTOFhitFirst, AliTOFPad* pixelArray , Int_t* iTOFpixel, Float_t* toftime, AliTOFRecHit* hitArray, Int_t& isHitOnFiredPad, Int_t& ipi
1562xel)
db9ba97f 1563{
1564 //
1565 // Read TOF hits for the current event and fill arrays
1566 //
1567 // Start loop on primary tracks in the hits containers
1568 //
1569 // Noise meaning in ReadTOFHits: we use the word 'noise' in the
1570 // following cases
1571 // - signals produced by secondary particles
1572 // - signals produced by the next hits (out of the first) of a given track
1573 // (both primary and secondary)
1574 // - signals produced by edge effect
1575
1576
1577 TParticle *particle;
1578 Int_t nHitOutofTofVolumes; // number of hits out of TOF GEANT volumes (it happens in very
1579 // few cases)
f9a28264 1580 Int_t * npixel = new Int_t[AliTOFConstants::fgkmaxtoftree]; // array used by TOFRecon for check on TOF geometry
db9ba97f 1581 Int_t npions=0; // number of pions for the current event
1582 Int_t nkaons=0; // number of kaons for the current event
1583 Int_t nprotons=0; // number of protons for the current event
1584 Int_t nelectrons=0;// number of electrons for the current event
1585 Int_t nmuons=0; // number of muons for the current event
1586 Float_t tofpos[3]; // TOF hit position and GEANT time
1587 Float_t zPad,xPad;
1588 Int_t nbytes = 0;
1589 Int_t ipart, nhits=0, nHitsFromPrimaries=0;
1590 Int_t ntotalTOFhits=0; // total number of TOF hits for the current event
1591 Int_t ipartLast=-1; // last track identifier
1592 Int_t iFirstHit; // flag to check if the current hit is the first hit on TOF for the
1593 // current track
1594 Int_t iNoiseHit=0; // flag used to tag noise hits (the noise meaning is reported in the
1595 // header of the ReadTOFHits method)
1596 Int_t nhitWithoutNoise;// number of hits not due to noise
1597 Int_t inoise=0,inoise2=0;
1598 Int_t nMultipleSignOnSamePad=0; // number of cases where a pad is fired more than one time
1599 Int_t nPixEdge=0; // additional pads fired due to edge effect in ReadTOFHits (local var)
1600 // array used for counting different types of primary particles
1601 Int_t particleTypeGEANT[50]={0,4,4,0,5,5,0,3,3,0,
1602 2,2,0,1,1,0,0,0,0,0,
1603 0,0,0,0,0,0,0,0,0,0,
1604 0,0,0,0,0,0,0,0,0,0,
1605 0,0,0,0,0,0,0,0,0,0};
1606 Int_t particleType,particleInTOFtype[6][3];
1607 for (Int_t i=0;i<6;i++) {
1608 for (Int_t j=0;j<3;j++) {
1609 particleInTOFtype[i][j]=0;
1610 }
1611 }
1612
5fff655e 1613 // speed-up the code
1614 treehits->SetBranchStatus("*",0); // switch off all branches
1615 treehits->SetBranchStatus("TOF*",1); // switch on only TOF
db9ba97f 1616
1617 for (Int_t track=0; track<ntracks;track++) { // starting loop on primary tracks for the current event
1618
1619 gAlice->ResetHits();
1620 nbytes += treehits->GetEvent(track);
1621 nhits = tofhits->GetEntriesFast();
1622
1623 ntotalTOFhits+=nhits;
1624
1625 // Start loop on hits connected to the current primary tracked particle
1626 // (including hits produced by secondary particles generaterd by the
1627 // current ptimary tracked particle)
1628 for (Int_t hit=0;hit<nhits;hit++) {
1629 AliTOFhit* tofHit = (AliTOFhit*)tofhits->UncheckedAt(hit);
1630 ipart = tofHit->GetTrack();
1631 if(ipart>=fMaxAllTracks) break;
1632 Float_t geantTime= tofHit->GetTof(); // it is given in [s]
1633 particle = (TParticle*)gAlice->Particle(ipart);
1634
1635 Int_t pdgCode=particle->GetPdgCode();
1636 // Only high momentum tracks (see fPBound value)
1637 // momentum components at vertex
1638 Float_t pxvtx = particle->Px();
1639 Float_t pyvtx = particle->Py();
1640 Float_t pzvtx = particle->Pz();
1641 Float_t pvtx = TMath::Sqrt(pxvtx*pxvtx+pyvtx*pyvtx+pzvtx*pzvtx);
1642 if(pvtx>fPBound) {
1643
1644 if(particle->GetFirstMother() < 0) nHitsFromPrimaries++; // count primaries
1645
1646 // x and y coordinates of the particle production vertex
1647 Float_t vx = particle->Vx();
1648 Float_t vy = particle->Vy();
1649 Float_t vr = TMath::Sqrt(vx*vx+vy*vy); // cylindrical radius of the particle production vertex
1650
1651 Float_t x = tofHit->X(); tofpos[0]=x;
1652 Float_t y = tofHit->Y(); tofpos[1]=y;
1653 Float_t z = tofHit->Z(); tofpos[2]=z;
b213b8bd 1654 /* var used for QA
db9ba97f 1655 Float_t tofradius = TMath::Sqrt(x*x+y*y); // radius cilindrical coordinate of the TOF hit
b213b8bd 1656 */
db9ba97f 1657 // momentum components (cosine) when striking the TOF
1658 Float_t pxtof = tofHit->GetPx();
1659 Float_t pytof = tofHit->GetPy();
1660 Float_t pztof = tofHit->GetPz();
1661 // scalar product indicating the direction of the particle when striking the TOF
b213b8bd 1662 /* var used for QA
db9ba97f 1663 // (>0 for outgoing particles)
1664 Float_t isGoingOut = (x*pxtof+y*pytof+z*pztof)/TMath::Sqrt(x*x+y*y+z*z);
b213b8bd 1665 */
db9ba97f 1666 Float_t momtof = tofHit->GetMom();
1667 // now momentum components when striking the TOF
1668 pxtof *= momtof;
1669 pytof *= momtof;
1670 pztof *= momtof;
1671 particleType=particleTypeGEANT[PDGtoGeantCode(pdgCode)-1];
1672 if(particleType) {
1673 particleInTOFtype[5][2]++;
1674 particleInTOFtype[particleType-1][2]++;
1675 }
1676 iFirstHit=0;
1677 // without noise hits
1678
1679 if(ipart!=ipartLast) {
1680 iFirstHit=1;
1681 toftime[ipart]=geantTime; //time [s]
1682 // tofMom[ipart]=momtof;
1683 ipartLast=ipart;
1684 if(particle->GetFirstMother() < 0) {
1685 Int_t abspdgCode=TMath::Abs(pdgCode);
1686 switch (abspdgCode) {
1687 case 211:
1688 npions++;
1689 break ;
1690 case 321:
1691 nkaons++;
1692 break ;
1693 case 2212:
1694 nprotons++;
1695 break ;
1696 case 11:
1697 nelectrons++;
1698 break ;
1699 case 13:
1700 nmuons++;
1701 break ;
1702 }
1703 }
1704 if(vr>fRadiusvtxBound) {
1705 if(particleType) {
1706 particleInTOFtype[5][1]++;
1707 particleInTOFtype[particleType-1][1]++;
1708 }
1709 inoise++;
1710 inoise2++;
1711 } else {
1712 if(particleType) {
1713 particleInTOFtype[5][0]++;
1714 particleInTOFtype[particleType-1][0]++;
1715 }
1716 }
1717 } else {
1718 inoise++;
1719 if(particleType) {
1720 particleInTOFtype[5][1]++;
1721 particleInTOFtype[particleType-1][1]++;
1722 }
1723 } //end if(ipart!=ipartLast)
1724
b9d0a01d 1725 IsInsideThePad(gMC,x,y,z,npixel,zPad,xPad);
db9ba97f 1726
1727 Int_t sec = tofHit->GetSector();
1728 Int_t pla = tofHit->GetPlate();
1729 Int_t str = tofHit->GetStrip();
1730 if(sec!=npixel[0] || pla!=npixel[1] || str!=npixel[2]){// check on volume
1731 cout << "sector" << sec << " npixel[0] " << npixel[0] << endl;
1732 cout << "plate " << pla << " npixel[1] " << npixel[1] << endl;
1733 cout << "strip " << str << " npixel[2] " << npixel[2] << endl;
1734 } // close check on volume
1735
1736 Int_t padz = tofHit->GetPadz();
1737 Int_t padx = tofHit->GetPadx();
1738 Float_t Zpad = tofHit->GetDz();
1739 Float_t Xpad = tofHit->GetDx();
1740
1741
1742 if (npixel[4]==0){
b9d0a01d 1743 IsInsideThePad(gMC,x,y,z,npixel,zPad,xPad);
db9ba97f 1744 if (npixel[4]==0){
1745 nHitOutofTofVolumes++;
1746 }
1747 } else {
1748 Float_t zStrip=AliTOFConstants::fgkZPad*(padz-0.5-0.5*AliTOFConstants::fgkNpadZ)+Zpad;
1749 if(padz!=npixel[3]) printf(" : Zpad=%f, padz=%i, npixel[3]=%i, zStrip=%f\n",Zpad,padz,npixel[3],zStrip);
1750 Float_t xStrip=AliTOFConstants::fgkXPad*(padx-0.5-0.5*AliTOFConstants::fgkNpadX)+Xpad;
1751
1752 Int_t nPlace[4]={0,0,0,0};
1753 nPlace[0]=(padz-1)*AliTOFConstants::fgkNpadX+padx;
1754
1755 Int_t nActivatedPads=0;
1756 Int_t nFiredPads=0;
1757 Bool_t isFired[4]={kFALSE,kFALSE,kFALSE,kFALSE};
1758 Float_t tofAfterSimul[4]={0.,0.,0.,0.};
1759 Float_t qInduced[4]={0.,0.,0.,0.};
1760 Float_t averageTime=0.;
1761
1762
1763 BorderEffect(zStrip,xStrip,geantTime*1.0e+09,nActivatedPads,nFiredPads,isFired,nPlace,qInduced,tofAfterSimul,averageTime); // simulate edge effect
1764
1765
1766 if(nFiredPads) {
1767 for(Int_t indexOfPad=0; indexOfPad<nActivatedPads; indexOfPad++) {
1768 if(isFired[indexOfPad]){// the pad has fired
1769 if(indexOfPad==0) {// the hit belongs to a fired pad
1770 isHitOnFiredPad++;
1771 hitArray[isHitOnFiredPad-1].SetHit(ipart,pdgCode,tofpos,momtof,vr,iFirstHit);
1772 iNoiseHit=0;
1773
1774 if(vr>fRadiusvtxBound || iFirstHit==0) iNoiseHit=1;
1775
1776 hitArray[isHitOnFiredPad-1].SetNoise(iNoiseHit);
1777 if(iFirstHit) kTOFhitFirst[ipart]=isHitOnFiredPad;
1778
1779 }// close - the hit belongs to a fired pad
1780
1781 Int_t iMapFirstIndex=AliTOFConstants::fgkNSectors*(npixel[1]-1)+npixel[0]-1;
1782 Int_t iMapValue=MapPixels[iMapFirstIndex][npixel[2]-1][nPlace[indexOfPad]-1];
1783
1784 if(iMapValue==0) {
1785 ipixel++;
1786 if(indexOfPad) {
1787 iNoiseHit=1;
1788 nPixEdge++;
1789 } else {
1790 iTOFpixel[ipart]=ipixel;
1791 }
1792
1793 if(ipixel>fMaxPixels){ // check on the total number of activated pads
1794 cout << "ipixel=" << ipixel << " > fMaxPixels=" << fMaxPixels << endl;
1795 return;
1796 } // close check on the number of activated pads
1797
1798 MapPixels[iMapFirstIndex][npixel[2]-1][nPlace[indexOfPad]-1]=ipixel;
1799 pixelArray[ipixel-1].SetGeom(npixel[0],npixel[1],npixel[2],nPlace[indexOfPad]);
1800 pixelArray[ipixel-1].SetTrack(ipart);
1801 if(iNoiseHit) {
1802 pixelArray[ipixel-1].AddState(1);
1803 } else {
1804 if(tofAfterSimul[indexOfPad]<0) cout << "Time of Flight after detector simulation is negative" << endl;
1805 pixelArray[ipixel-1].AddState(10);
1806 }
1807
1808 pixelArray[ipixel-1].SetTofChargeHit(tofAfterSimul[indexOfPad],qInduced[indexOfPad],geantTime*1.0e+09,isHitOnFiredPad);
1809 } else { //else if(iMapValue==0)
1810 if(indexOfPad==0) iTOFpixel[ipart]=iMapValue;
1811 nMultipleSignOnSamePad++;
1812
1813 if(tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetRealTime() ) {
1814 pixelArray[iMapValue-1].SetTrack(ipart);
1815 // if(indexOfPad==0) pixelArray[iMapValue-1].SetTrack(ipart);
1816 if(indexOfPad) iNoiseHit=1;
1817 if(iNoiseHit) {
1818 pixelArray[iMapValue-1].AddState(1);
1819 } else {
1820 pixelArray[iMapValue-1].AddState(10);
1821 }
1822 pixelArray[iMapValue-1].SetRealTime(tofAfterSimul[indexOfPad]);
1823 pixelArray[iMapValue-1].SetGeantTime(geantTime*1.0e+09);
1824 pixelArray[iMapValue-1].SetHit(isHitOnFiredPad);
1825 } // close if(tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetTime() )
1826 } //end of Pixel filling
1827 } // close if(isFired[indexOfPad])
1828 } //end loop on activated pads indexOfPad
1829 } // close if(nFiredPads)
1830 } //end of hit with npixel[3]!=0
1831 } //high momentum tracks
1832 } //end on TOF hits
1833 } //end on primary tracks
1834
1835
1836 if(fdbg) {
1837 cout << ntotalTOFhits << " - total number of TOF hits " << nHitsFromPrimaries << " - primary " << endl;
1838 cout << inoise << " - noise hits, " << inoise2<< " - first crossing of a track with Rvtx>" << fRadiusvtxBound << endl;
1839 // cout << inoise << " - noise hits (" << 100.*inoise/ihit << " %), " << inoise2
1840 //<< " - first crossing of a track with Rvtx>" << RVTXBOUND << endl;
1841 nhitWithoutNoise=isHitOnFiredPad;
1842
1843 cout << ipixel << " fired pixels (" << nMultipleSignOnSamePad << " multiple fired pads, " << endl;
1844 //j << " fired by noise, " << j1 << " noise+track)" << endl;
1845 printf(" %i additional pads are fired due to edge effect\n",nPixEdge);
1846 cout << npions << " primary pions reached TOF" << endl;
1847 cout << nkaons << " primary kaons reached TOF" << endl;
1848 cout << nprotons << " primary protons reached TOF" << endl;
1849 cout << nelectrons<<" primary electrons reached TOF" << endl;
1850 cout << nmuons << " primary muons reached TOF" << endl;
1851 cout << "number of TOF hits for different species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl;
1852 cout << " first number - track hits, second - noise ones, third - all" << endl;
1853 for (Int_t i=0;i<6;i++) cout << i+1 << " " << particleInTOFtype[i][0] << " " << particleInTOFtype[i][1] << " " << particleInTOFtype[i][2] << endl;
1854
1855 Int_t primaryReachedTOF[6];
1856 primaryReachedTOF[0]=npions;
1857 primaryReachedTOF[1]=nkaons;
1858 primaryReachedTOF[2]=nprotons;
1859 primaryReachedTOF[3]=nelectrons;
1860 primaryReachedTOF[4]=nmuons;
1861 primaryReachedTOF[5]=npions+nkaons+nprotons+nelectrons+nmuons;
1862
1863 cout << " Reading TOF hits done" << endl;
1864 }
1865
f9a28264 1866 delete [] npixel;
db9ba97f 1867}
1868
1869//____________________________________________________________________________
1870void AliTOFReconstructioner::AddNoiseFromOuter(Option_t *option, Int_t ***MapPixels, AliTOFPad* pixelArray , AliTOFRecHit* hitArray, Int_t& isHitOnFiredPad, Int_t& ipixel)
1871{
1872 //
1873 // Add noise hits from outer regions (forward and backward) according
1874 // to parameterized fZNoise distribution (to be used with events
1875 // generated in the barrel region)
1876
f9a28264 1877 Float_t * zLen = new Float_t[AliTOFConstants::fgkNPlates+1];
1878 Float_t * zStrips = new Float_t[AliTOFConstants::fgkNPlates];
db9ba97f 1879 zStrips[0]=(Float_t) (AliTOFConstants::fgkNStripC);
1880 zStrips[1]=(Float_t) (AliTOFConstants::fgkNStripB);
1881 zStrips[2]=(Float_t) (AliTOFConstants::fgkNStripA);
1882 zStrips[3]=(Float_t) (AliTOFConstants::fgkNStripB);
1883 zStrips[4]=(Float_t) (AliTOFConstants::fgkNStripC);
1884
1885 zLen[5]=AliTOFConstants::fgkzlenA*0.5+AliTOFConstants::fgkzlenB+AliTOFConstants::fgkzlenC;
1886 zLen[4]=zLen[5]-AliTOFConstants::fgkzlenC;
1887 zLen[3]=zLen[4]-AliTOFConstants::fgkzlenB;
1888 zLen[2]=zLen[3]-AliTOFConstants::fgkzlenA;
1889 zLen[1]=zLen[2]-AliTOFConstants::fgkzlenB;
1890 zLen[0]=zLen[1]-AliTOFConstants::fgkzlenC;
1891
1892
1893 Int_t isector; // random sector number
1894 Int_t iplate; // random plate number
1895 Int_t istrip; // random strip number in the plate
1896 Int_t ipadAlongX; // random pad number along x direction
1897 Int_t ipadAlongZ; // random pad number along z direction
1898 Int_t ipad;
1899 Int_t nPixEdge=0; // additional pads fired due to edge effect when adding noise from outer
1900 // regions
1901
1902 // x -> time of flight given in ns
1903 TF1 *noiseTof = new TF1("noiseTof","exp(-x/20)",0,100);
1904
1905 if(strstr(option,"pp")){
1906 fZnoise = new TF1("fZnoise","257.8-0.178*x-0.000457*x*x",-AliTOFConstants::fgkMaxhZtof,AliTOFConstants::fgkMaxhZtof);
1907 }
1908 if(strstr(option,"Pb-Pb")){
1909 fZnoise = new TF1("fZnoise","182.2-0.09179*x-0.0001931*x*x",-AliTOFConstants::fgkMaxhZtof,AliTOFConstants::fgkMaxhZtof);
1910 }
1911
1912 if(fNoise) {
1913 if(fdbg) cout << " Start adding additional noise hits from outer regions" << endl;
1914
1915 for(Int_t i=0;i<fNoise;i++) {
1916
1917 isector=(Int_t) (AliTOFConstants::fgkNSectors*gRandom->Rndm())+1; //the sector number
1918 // non-flat z-distribution of additional hits
1919 Float_t zNoise=fZnoise->GetRandom();
1920
1921 // holes for PHOS and HMPID
1922 if(((AliTOF *) gAlice->GetDetector("TOF"))->IsVersion()==2) {
1923 // to be checked the holes case
1924 if(isector>12 && isector<16) { // sectors 13,14,15 - RICH
1925 do {
1926 iplate=(Int_t) (AliTOFConstants::fgkNPlates*gRandom->Rndm())+1;
1927 } while (iplate==2 || iplate==3 || iplate==4);
1928 // } else if(isector>11 && isector<17) { // sectors 12,13,14,15,16 - PHOS
1929 } else if(isector>2 && isector<8) { // sectors 3,4,5,6,7 - PHOS
1930 do {
1931 iplate=(Int_t) (AliTOFConstants::fgkNPlates*gRandom->Rndm())+1;
1932 } while (iplate==3);
1933 } else {
1934 iplate=(Int_t) (AliTOFConstants::fgkNPlates*gRandom->Rndm())+1;
1935 }
1936 } else {
1937 iplate=0;
1938 do {
1939 iplate++;
1940 } while(zNoise>zLen[iplate]);
1941 }
1942 // end of holes
1943
1944 if(iplate<1 || iplate>5) {
1945 printf(" iplate<1 or iplate>5, iplate=%i\n",iplate);
1946 return;
1947 }
1948
1949 Float_t nStripes=0;
1950 if(iplate>1) {
1951 for (Int_t i=0;i<iplate-1;i++) {
1952 nStripes += zStrips[i];
1953 }
1954 }
1955
b213b8bd 1956 istrip=(Int_t)((zNoise-zLen[iplate-1])/((zLen[iplate]-zLen[iplate-1])/zStrips[iplate-1])); //the strip number in the plate
db9ba97f 1957 istrip++;
1958
1959 ipadAlongX = (Int_t)(AliTOFConstants::fgkNpadX*gRandom->Rndm())+1;
1960 ipadAlongZ = (Int_t)(AliTOFConstants::fgkNpadZ*gRandom->Rndm())+1;
1961 ipad=(Int_t)(ipadAlongZ-1)*AliTOFConstants::fgkNpadX+ipadAlongX; //the pad number
1962
1963 Float_t xStrip=(ipadAlongX-1)*AliTOFConstants::fgkXPad+AliTOFConstants::fgkXPad*gRandom->Rndm()-0.5*AliTOFConstants::fgkNpadX*AliTOFConstants::fgkXPad;//x-coor.in the strip frame
1964 Float_t zStrip=(ipadAlongZ-1)*AliTOFConstants::fgkZPad+AliTOFConstants::fgkZPad*gRandom->Rndm()-0.5*AliTOFConstants::fgkNpadZ*AliTOFConstants::fgkZPad;//z-coor.in the strip frame
1965
1966 Int_t nPlace[4]={0,0,0,0};
1967 nPlace[0]=ipad;
1968
1969 Int_t nActivatedPads=0;
1970 Int_t nFiredPads=0;
1971 Bool_t isFired[4]={kFALSE,kFALSE,kFALSE,kFALSE};
1972 Float_t tofAfterSimul[4]={0.,0.,0.,0.};
1973 Float_t qInduced[4]={0.,0.,0.,0.};
1974 Float_t averageTime=0.;
1975 Float_t toffornoise=10.+noiseTof->GetRandom(); // 10 ns offset + parameterization [ns]
1976
1977 BorderEffect(zStrip,xStrip,toffornoise,nActivatedPads,nFiredPads,isFired,nPlace,qInduced,tofAfterSimul,averageTime); // simulate edge effect
1978
1979 if(nFiredPads) {
1980 for(Int_t indexOfPad=0; indexOfPad<nActivatedPads; indexOfPad++) {
1981 if(isFired[indexOfPad]){// the pad has fired
1982
1983 if(indexOfPad==0) {// the hit belongs to a fired pad
1984 isHitOnFiredPad++;
1985 hitArray[isHitOnFiredPad-1].SetX(0.);
1986 hitArray[isHitOnFiredPad-1].SetY(0.);
1987 hitArray[isHitOnFiredPad-1].SetZ(zNoise);
1988 hitArray[isHitOnFiredPad-1].SetNoise(1);
1989 } // close if(indexOfPad==0)
1990
1991 ipad = nPlace[indexOfPad];
1992
1993 Int_t iMapValue=MapPixels[AliTOFConstants::fgkNSectors*(iplate-1)+isector-1][istrip-1][ipad-1];
1994
1995 if(iMapValue==0) {
1996 ipixel++;
1997 if(indexOfPad) nPixEdge++;
1998 MapPixels[AliTOFConstants::fgkNSectors*(iplate-1)+isector-1][istrip-1][ipad-1]=ipixel;
1999 pixelArray[ipixel-1].SetGeom(isector,iplate,istrip,ipad);
2000 pixelArray[ipixel-1].AddState(1);
2001 pixelArray[ipixel-1].SetRealTime(tofAfterSimul[indexOfPad]);
2002 pixelArray[ipixel-1].SetHit(isHitOnFiredPad);
2003 } else if( tofAfterSimul[indexOfPad] < pixelArray[iMapValue-1].GetRealTime() ) {
2004 pixelArray[iMapValue-1].SetTrack(-1);
2005 pixelArray[iMapValue-1].AddState(1);
2006 pixelArray[iMapValue-1].SetRealTime(tofAfterSimul[indexOfPad]);
2007 pixelArray[iMapValue-1].SetHit(isHitOnFiredPad);
2008 } //end of if(iMapValue==0)
2009
2010 }// close if(isFired[indexOfPad])
2011 } //end loop on activated pads indexOfPad
2012 } // close if(nFiredPads)
2013 } //end of NOISE cycle
2014 }
2015
2016 // free used memory
2017 if (fZnoise)
2018 {
2019 delete fZnoise;
2020 fZnoise = 0;
2021 }
2022
2023 if (noiseTof)
2024 {
2025 delete noiseTof;
2026 noiseTof = 0;
2027 }
2028
2029 Int_t nNoiseSignals=0;
2030 Int_t nAll=0;
2031 for(Int_t idummy=1; idummy<ipixel+1; idummy++) {
2032 if(hitArray[pixelArray[idummy-1].GetHit()-1].GetNoise()==1) {
2033 nNoiseSignals++;
2034 if(pixelArray[idummy-1].GetState()>10) nAll++;
2035 }
2036 }
2037
2038 if(fdbg) {
2039 cout << " after adding " << fNoise << " noise hits: " << ipixel << " fired pixels (" << nNoiseSignals << " fired by noise, " << nAll << " noise+track)" << endl;
2040 printf(" %i additional pixels are fired by noise due to edge effect\n",nPixEdge);
2041 cout << " End of adding additional noise hits from outer regions" << endl;
2042 }
2043
2044 Float_t occupancy;
2045 // numberOfPads for AliTOFV4 (Full coverage)
2046 // - to be upgraded checking the used TOF version -
2047 Float_t numberOfPads=AliTOFConstants::fgkPadXSector*AliTOFConstants::fgkNSectors;
2048 occupancy=100.*ipixel/numberOfPads; // percentage of fired pads
2049 printf(" Overall TOF occupancy (percentage of fired pads after adding noise) = %f\n",occupancy);
f9a28264 2050 delete [] zLen;
2051 delete [] zStrips;
db9ba97f 2052
2053}
2054
2055
2056//____________________________________________________________________________
2057void AliTOFReconstructioner::SetMinDistance(AliTOFRecHit* hitArray, Int_t ilastEntry)
2058{
2059 //
2060 // Set the distance to the nearest hit for hitArray
2061 // ilastEntry is the index of the last entry of hitArray
2062
2063 // starting the setting for the distance to the nearest TOFhit (cm)
2064 for(Int_t i=0; i<ilastEntry; i++) {
2065
2066 if(hitArray[i].GetFirst()==1 && hitArray[i].GetNoise()==0) { // select the first hit of the track
2067 // hits are not due to noise
2068 Float_t minDistance=10000.,squareDistance; // current values of the (square) distance
2069 Int_t jAtMin=0; // index of the hit nearest to the i-th hit
2070 Float_t xhit=hitArray[i].X(); // x coordinate for i-th hit
2071 Float_t yhit=hitArray[i].Y(); // y coordinate for i-th hit
2072 Float_t zhit=hitArray[i].Z(); // z coordinate for i-th hit
2073 // was for(Int_t j=0; j<isHitOnFiredPad; j++) {
2074 for(Int_t j=0; j<ilastEntry; j++) {
2075 if(i!=j) {
2076 squareDistance=(hitArray[j].X()-xhit)*(hitArray[j].X()-xhit)+
2077 (hitArray[j].Y()-yhit)*(hitArray[j].Y()-yhit)+
2078 (hitArray[j].Z()-zhit)*(hitArray[j].Z()-zhit);
2079 if(squareDistance<minDistance) {
2080 minDistance=squareDistance;
2081 jAtMin=j;
2082 }
2083 }
2084 }
2085 minDistance=TMath::Sqrt(minDistance);
2086 hitArray[i].SetRmin(minDistance);
2087 if(minDistance==0.) printf(" Rmin=0, i=%i, j=%i, x=%f,y=%f,z=%f\n",i,jAtMin,xhit,yhit,zhit);// it cannot happen
2088 }
2089 }
2090
2091}
2092
2093// these lines has to be commented till TPC will provide fPx fPy fPz
2094// and fL in AliTPChit class
2095//____________________________________________________________________________
2096/*
2097void AliTOFReconstructioner::ReadTPCHits(Int_t ntracks, TTree* treehits, TClonesArray* tpchits, Int_t* iTrackPt, Int_t* iparticle, Float_t* ptTrack, AliTOFTrack* trackArray, Int_t& itrack)
2098{
2099 //
2100 // Read TPC hits for the current event
2101 //
2102 TParticle *particle=0;
2103 Int_t npions=0; // number of pions for the current event
2104 Int_t nkaons=0; // number of kaons for the current event
2105 Int_t nprotons=0; // number of protons for the current event
2106 Int_t nelectrons=0;// number of electrons for the current event
2107 Int_t nmuons=0; // number of muons for the current event
2108 Int_t ntotalTPChits=0; // total number of TPC hits for the current event
2109 Int_t idummy=-1; // dummy var used to count double hit TPC cases
2110 Int_t nTpcDoubleHitsLastRow=0; // number of double TPC hits in the last pad row
2111 Int_t nTpcHitsLastRow=0; // number of TPC hits in the last pad row
2112 Float_t trdpos[2]={0.,0.};
2113 Float_t pos[3]; // TPC hit position
2114 Float_t mom[3]; // momentum components in the last TPC row
2115 Float_t pt=0., tpclen; // pt: transverse momentum in the last TPC row
2116 Int_t nbytes = 0;
2117 Int_t ipart=0, nhits=0, iprim=0;
2118
2119 itrack=0; // itrack: total number of selected TPC tracks
2120
5fff655e 2121 // speed-up the code
2122 treehits->SetBranchStatus("*",0); // switch off all branches
2123 treehits->SetBranchStatus("TPC*",1); // switch on only TPC
2124
db9ba97f 2125 for (Int_t track=0; track<ntracks;track++) {
2126 gAlice->ResetHits();
2127 nbytes += treehits->GetEvent(track);
2128
2129
2130 nhits = tpchits->GetEntriesFast();
2131
2132 for (Int_t hit=0;hit<nhits;hit++) {
2133 ntotalTPChits++;
2134 AliTPChit* tpcHit = (AliTPChit*)tpchits->UncheckedAt(hit);
2135 Int_t row = tpcHit->fPadRow;
2136 ipart = tpcHit->GetTrack();
2137 if(ipart>=fMaxAllTracks) break;
2138 particle = (TParticle*)gAlice->Particle(ipart);
2139 Int_t pdgCode=particle->GetPdgCode();
2140 // only high momentum tracks
2141 // momentum components at production vertex
2142 Float_t pxvtx = particle->Px();
2143 Float_t pyvtx = particle->Py();
2144 Float_t pzvtx = particle->Pz();
2145 Float_t pvtx = TMath::Sqrt(pxvtx*pxvtx+pyvtx*pyvtx+pzvtx*pzvtx);
2146 if(pvtx>fPBound && row == fLastTPCRow) {
2147 Float_t vx = particle->Vx();
2148 Float_t vy = particle->Vy();
2149 Float_t vr = TMath::Sqrt(vx*vx+vy*vy);
2150 Float_t x = tpcHit->X();
2151 Float_t y = tpcHit->Y();
2152 Float_t z = tpcHit->Z();
2153 pos[0]=x; pos[1]=y; pos[2]=z;
2154
2155 Float_t pxtpc = tpcHit->fPx;
2156 Float_t pytpc = tpcHit->fPy;
2157 Float_t pztpc = tpcHit->fPz;
2158 mom[0]=pxtpc; mom[1]=pytpc; mom[2]=pztpc;
2159 Float_t momtpc = TMath::Sqrt(pxtpc*pxtpc+pytpc*pytpc+pztpc*pztpc);
2160
2161 if(x*pxtpc+y*pytpc>0) { // only tracks going out of TPC
2162
2163 Float_t isoutgoing = x*pxtpc+y*pytpc+z*pztpc;
2164 isoutgoing /= (momtpc*TMath::Sqrt(x*x+y*y+z*z));
2165 tpclen = tpcHit->fL;
2166
2167
2168 if(ipart!=idummy) {
2169 if(particle->GetFirstMother() < 0) {
2170 Int_t abspdgCode=TMath::Abs(pdgCode);
2171 switch (abspdgCode) {
2172 case 211:
2173 npions++;
2174 break ;
2175 case 321:
2176 nkaons++;
2177 break ;
2178 case 2212:
2179 nprotons++;
2180 break ;
2181 case 11:
2182 nelectrons++;
2183 break ;
2184 case 13:
2185 nmuons++;
2186 break ;
2187 }
2188 } // close if(particle->GetFirstMother() < 0)
2189 } // close if(ipart!=idummy)
2190
2191 if(gRandom->Rndm()<fTrackingEfficiency && vr<fRadiusvtxBound && ipart!=idummy) {
2192
2193 itrack++;
2194 if(particle->GetFirstMother() < 0) iprim++;
2195
2196 if(itrack>fMaxTracks) {
2197 cout << "itrack=" << itrack << " > MAXTRACKS=" << fMaxTracks << endl;
2198 return;
2199 } // close if(itrack>fMaxTracks)
2200
2201
2202 iparticle[ipart]=itrack;
2203
2204 trackArray[itrack-1].SetTrack(ipart,pvtx,pdgCode,tpclen,pos,mom,trdpos);
2205
2206 pt=TMath::Sqrt(pxtpc*pxtpc+pytpc*pytpc); // pt: transverse momentum at TPC
2207 // Filling iTrackPt[MAXTRACKS] by itrack ordering on Pt
2208 if(itrack==1) {
2209 iTrackPt[itrack-1]=itrack;
2210 ptTrack[itrack-1]=pt;
2211 } else {
2212 for (Int_t i=0; i<itrack-1; i++) {
2213 if(pt>ptTrack[i]) {
2214 for(Int_t j=i; j<itrack-1; j++) {
2215 Int_t k=itrack-1+i-j;
2216 iTrackPt[k]= iTrackPt[k-1];
2217 ptTrack[k] = ptTrack[k-1];
2218 }
2219 iTrackPt[i]=itrack;
2220 ptTrack[i]=pt;
2221 break;
2222 }
2223 if(i==itrack-2) {
2224 iTrackPt[itrack-1]=itrack;
2225 ptTrack[itrack-1]=pt;
2226 }
2227 }
2228 }
2229
2230 } //end of itrack
2231 if(vr>fRadiusvtxBound) nTpcHitsLastRow++;
2232 if(ipart==idummy) nTpcDoubleHitsLastRow++;
2233 idummy=ipart;
2234 } // close if(x*px+y*py>0)
2235 } // close if(pvtx>fPBound && row == fLastTPCRow)
2236 } //end of hits
2237 } // close loop on tracks
2238
2239
2240 if(fdbg) {
2241 cout << ntotalTPChits << " TPC hits in the last TPC row " << fLastTPCRow << endl;
2242 cout << " " << nTpcHitsLastRow << " - hits with Rvtx>fRadiusvtxBound=" << fRadiusvtxBound << endl;
2243 cout << " " << nTpcDoubleHitsLastRow << " double TPC hits" << endl;
2244 cout << itrack << " - extracted TPC tracks " << iprim << " - primary" << endl;
2245 cout << npions << " primary pions reached TPC" << endl;
2246 cout << nkaons << " primary kaons reached TPC" << endl;
2247 cout << nprotons << " primary protons reached TPC" << endl;
2248 cout << nelectrons<< " primary electrons reached TPC" << endl;
2249 cout << nmuons << " primary muons reached TPC" << endl;
2250 } // if(fdbg)
2251
2252 Int_t primaryInTPC[6]={0,0,0,0,0,0};
2253 primaryInTPC[0]=npions;
2254 primaryInTPC[1]=nkaons;
2255 primaryInTPC[2]=nprotons;
2256 primaryInTPC[3]=nelectrons;
2257 primaryInTPC[4]=nmuons;
2258 primaryInTPC[5]=npions+nkaons+nprotons+nelectrons+nmuons;
2259
2260 if(fdbg) {
2261 printf(" contents of iTrackPt[MAXTRACKS],PtTrack[MAXTRACKS]\n");
2262 for (Int_t i=0; i<itrack; i++) {
2263 printf(" %i : iTrackPt=%i, PtTrack=%f\n",i+1,iTrackPt[i],ptTrack[i]);
2264 }
2265 printf(" Check ordered transverse momentum array\n");
2266 for (Int_t i=itrack-1; i>=0; i--) {
2267 printf(" %i : iTrackPt=%i, PtTrack=%f\n",i+1,iTrackPt[i],ptTrack[i]);
2268 }
2269 }// if(fdbg)
2270
2271}
2272*/
2273//____________________________________________________________________________
2274void cylcor(Float_t& x, Float_t& y) {
2275 Float_t rho,phi;
2276
2277 rho=TMath::Sqrt(x*x+y*y);
2278 phi=0.;
2279 if(TMath::Abs(x)>0. || TMath::Abs(y)>0.) phi=TMath::ATan2(y,x);
2280 if(phi<0.) phi=phi+2.*TMath::Pi();
2281 x=rho;
2282 y=phi;
2283
2284}
2285
2286//____________________________________________________________________________
2287void AliTOFReconstructioner::Matching(AliTOFTrack* trackArray, AliTOFRecHit* hitArray, Int_t ***mapPixels, AliTOFPad* pixelArray, Int_t* kTOFhitFirst, Int_t& ipixel, Int_t* iTrackPt, Int_t* iTOFpixel, Int_t ntotTpcTracks)
2288{
f9a28264 2289 Int_t TestTracks,iTestTrack,itest,wPixel=0,itestc;
2290 Int_t * ntest = new Int_t[fMaxTestTracks];
2291 Int_t * testPixel = new Int_t[fMaxTestTracks];
2292 Float_t wLength=0.,wRho=0.,wZ=0.;
2293 Float_t * testLength = new Float_t[fMaxTestTracks];
2294 Float_t * testRho = new Float_t[fMaxTestTracks];
2295 Float_t * testZ = new Float_t[fMaxTestTracks];
2296 Float_t weight;
2297 Float_t * testWeight = new Float_t[fMaxTestTracks];
db9ba97f 2298 Float_t rotationFactor,phi0,coslam,sinlam,helixRadius,xHelixCenter,yHelixCenter,zHelixCenter,helixFactor;
2299 Int_t npixel[5],iMapValue,iwork1,iwork2,iwork3,iwork4,ihit=0;
2300 Int_t charge[48]={ 0, 1,-1, 0, 1,-1, 0, 1,-1, 0,
2301 1,-1, 0, 1,-1, 0, 0, 0, 1, 0,
2302 -1, 0,-1,-1, 0, 0,-1, 0, 1, 0,
2303 1, 1, 0, 0, 1,-1, 0, 0, 1,-1,
2304 1, 1,-1, 0, 1, 1, 2, 0};
2305 Float_t theta0,gpx,gpy,gpz,gp,gpt,gtheta,gx,gy,gz,gr,gxLast,gyLast,gzLast,chargeField;
2306 Float_t sumOfTheta=0.,weightTestTracksOutTof[4];
2307 Float_t s,ds,xRespectToHelixCenter,yRespectToHelixCenter,deltaRadius,fp,xp,yp,grho;
2308 Float_t mass,energy,g;
2309 Int_t itrack=0,itr,particleCharge,istep,iplate=0,iPadAlongX=0;
2310 Int_t itra,t34=0,t32=0,t44=0,t43=0,t42=0;
2311 Int_t wstate=0,m2state=0,wPix;
2312 Int_t idelR=0,idelR1=0,idelR2=0,iRmin=0,iRmin1=0,iRmin2=0;
2313 Float_t massArray[50] = {0.0,0.00051,0.00051,0.0,0.1057,0.1057,0.135,0.1396,0.1396,0.4977,
2314 0.4936,0.4936,0.9396,0.9383,0.9383,0.4977,0.5488,1.1156,1.1894,1.1926,1.1926,
2315 1.3149,1.3213,1.6724,0.9396,1.1156,1.1894,1.1926,1.1974,1.3149,
2316 0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
2317 Float_t delR;
2318 Float_t radius,area,normR,normS,cosAngl;
2319 Int_t iPlateFirst,iTestGmax=0;
2320 Int_t fstate,iPrintM1=0,iPrintM2=0;
2321 Float_t gxExtrap=0.,gyExtrap=0.,gzExtrap=0.;
2322 Float_t avSigZ=0,avSigRPHI=0,avSigP=0,avSigPHI=0,avSigTHETA=0;
2323
2324 Float_t gxW,gyW,gzW;
2325 Float_t length0;
2326 Float_t snr=0;
2327 Int_t indexOfTestTrack;
2328 Float_t zPad,xPad;
2329 Int_t istate=0,imax=0,match,iMaxTestTracksOutTof=0,matchw;
2330 Float_t w,wmax=0.,inverseOfParticleSpeed,w2,smat[9],largestWeightTracksOutTof,sw;
2331 Float_t sumWeightTracksOutTof,sGeomWeigth;
2332 Int_t imatched;
2333 Int_t m10=0,m20=0,m22=0,m23=0;
2334 Int_t PRINT=0;
2335 TParticle *particle;
2336
2337 Float_t time=0.;
2338 itr=ntotTpcTracks;
2339 printf(" itr=%i\n",itr);
2340 for (itra=1; itra<itr+1; itra++) {
2341
2342 Int_t itrack=iTrackPt[itra-1];
2343 if(itrack==0) printf(" iTrackPt[itra-1]=0 for itra=%i\n",itra);
2344 Int_t ipart=trackArray[itrack-1].GetTrack();
2345 Float_t pvtx=trackArray[itrack-1].GetP();
2346 Int_t pdgCode=trackArray[itrack-1].GetPdgCode();
2347 Float_t tpclength=trackArray[itrack-1].GetlTPC();
2348 Float_t x=trackArray[itrack-1].GetRxTPC();
2349 Float_t y=trackArray[itrack-1].GetRyTPC();
2350 Float_t z=trackArray[itrack-1].GetRzTPC();
b213b8bd 2351 /* vars used for QA
db9ba97f 2352 Float_t RxTPC=x;
2353 Float_t RyTPC=y;
2354 Float_t RzTPC=z;
b213b8bd 2355 */
db9ba97f 2356 Float_t Wx=x;
2357 Float_t Wy=y;
2358 Float_t Wz=z;
2359 Float_t px=trackArray[itrack-1].GetPxTPC();
2360 Float_t py=trackArray[itrack-1].GetPyTPC();
2361 Float_t pz=trackArray[itrack-1].GetPzTPC();
b213b8bd 2362 /* vars used for QA
db9ba97f 2363 Float_t pxTPC=px;
2364 Float_t pyTPC=py;
2365 Float_t pzTPC=pz;
b213b8bd 2366 */
db9ba97f 2367 Float_t p = TMath::Sqrt(px*px+py*py+pz*pz);
b213b8bd 2368 /* var used for QA
db9ba97f 2369 Float_t pTPC=p;
b213b8bd 2370 */
db9ba97f 2371 Float_t rho = TMath::Sqrt(x*x+y*y);
2372 Float_t phi=0.;
2373 if(TMath::Abs(x)>0. || TMath::Abs(y)>0.) phi=TMath::ATan2(y,x);
2374 if(phi<0.) phi=phi+2.*TMath::Pi();
b213b8bd 2375 /* var used for QA
db9ba97f 2376 Float_t phiTPC=phi*kRaddeg;
b213b8bd 2377 */
db9ba97f 2378 if(fSigmavsp) {
2379 if(p==0) printf(" p=%f in g=0.022/p\n",p);
2380 g=0.022/p;
2381 avSigRPHI += g; // (cm)
2382 if(rho==0) printf(" rho=%f in phi += g*gRandom->Gaus()/rho\n",rho);
2383 phi += g*gRandom->Gaus()/rho;
2384 } else {
2385 if(rho==0) printf(" rho=%f in phi += (SIGMARPHI*gRandom->Gaus()/rho\n",rho);
2386 phi += (fSigmarphi*gRandom->Gaus()/rho);
2387 }
2388 x=rho*TMath::Cos(phi);
2389 y=rho*TMath::Sin(phi);
b213b8bd 2390 /* var used for QA
db9ba97f 2391 Float_t zTPC=z;
b213b8bd 2392 */
db9ba97f 2393 if(fSigmavsp) {
2394 if(p==0) printf(" p=%f in g=0.0275/p\n",p);
2395 g=0.0275/p;
2396 avSigZ += g; // (cm)
2397 z += g*gRandom->Gaus();
2398 } else {
2399 z += fSigmaZ*gRandom->Gaus();
2400 }
2401
2402 // smearing on TPC momentum
2403
2404 {
2405 Float_t pmom,phi,theta,arg;
2406
2407 pmom=TMath::Sqrt(px*px+py*py+pz*pz);
2408 phi=0.;
2409 if(TMath::Abs(px)>0. || TMath::Abs(py)>0.) phi=TMath::ATan2(py,px);
2410 if(phi<0.) phi=phi+2*TMath::Pi();
2411 arg=1.;
2412 if(pmom>0.) arg=pz/pmom;
2413 theta=0.;
2414 if(TMath::Abs(arg)<=1.) theta=TMath::ACos(arg);
2415
2416 if(fSigmavsp) {
2417 if(pmom<=0) printf(" pmom=%f in g = TMath::Abs(TMath::Log(pmom)/TMath::Log(10)+0.5)/0.7\n",pmom);
2418 g = TMath::Abs(TMath::Log(pmom)/TMath::Log(10)+0.5)/0.7;
2419 g = 0.01*(g*g*g+1.5)*1.24;
2420 avSigP += g;
2421 pmom *= (1+g*gRandom->Gaus());
2422
2423 if(p<10) {
2424 if(pmom<=0) printf(" pmom=%f in g = 1-TMath::Log(pmom)/TMath::Log(10)\n",pmom);
2425 g = 1-TMath::Log(pmom)/TMath::Log(10);
2426 g = 0.001*(g*g*g+0.3)*0.65; // (radian)
2427 } else {
2428 g = 0.001*0.3*0.65;
2429 }
2430 avSigPHI += g;
2431 phi += g*gRandom->Gaus();
2432 avSigTHETA += g;
2433 theta += g*gRandom->Gaus();
2434
2435 } else {
2436 pmom *= (1+fSigmap*gRandom->Gaus());
2437 phi += fSigmaPhi*gRandom->Gaus();
2438 theta += fSigmaTheta*gRandom->Gaus();
2439 }
2440 gxW=px;
2441 gyW=py;
2442 gzW=pz;
2443
2444 px=pmom*TMath::Sin(theta)*TMath::Cos(phi);
2445 py=pmom*TMath::Sin(theta)*TMath::Sin(phi);
2446 pz=pmom*TMath::Cos(theta);
2447
2448
2449 if(x*px+y*py<=0) {
2450 x=Wx;
2451 y=Wy;
2452 z=Wz;
2453 px=gxW;
2454 py=gyW;
2455 pz=gzW;
2456 }// if(x*px+y*py<=0)
2457 }
2458
2459 p = TMath::Sqrt(px*px+py*py+pz*pz);
2460
2461 particleCharge=charge[PDGtoGeantCode(pdgCode)-1];
2462 mass=massArray[PDGtoGeantCode(pdgCode)-1];
2463 mass=massArray[8-1]; //we take pion mass for all tracks
2464 // mass=massArray[14-1]; //here we take proton mass for all tracks
2465 energy=TMath::Sqrt(p*p+mass*mass);
2466 chargeField=particleCharge*fField;
2467
2468 g=fRadLenTPC/( (x*px+y*py)/(rho*p) );
2469
2470 if(g<=0) printf(" error, g<=0: g=%f, itra=%i, x,y,px,py=%f, %f, %f, %f\n",g,itra,x,y,px,py);
2471
2472 theta0=13.6*0.001*TMath::Sqrt(g)*(1.+0.038*TMath::Log(g))*energy/(p*p);
2473
2474
2475 // start Loop on test tracks
2476 sumOfTheta=0.;
2477 for(Int_t i=0;i<4;i++) {
2478 weightTestTracksOutTof[i]=0.;
2479 }
2480
2481 itest=0;
2482 for(Int_t i=0;i<fMaxTestTracks;i++) {
2483 ntest[i]=0;
2484 testPixel[i]=0;
2485 testLength[i]=0.;
2486 testRho[i]=0.;
2487 testZ[i]=0.;
2488 testWeight[i]=0.;
2489 }
2490
2491 iPlateFirst=0;
2492 TestTracks=0;
2493 iTestTrack=0;
2494 iTestGmax=0;
2495
2496 length0=0;
2497
2498 for (indexOfTestTrack=0; indexOfTestTrack<fMaxTestTracks; indexOfTestTrack++) {
2499
2500 iTestTrack++;
2501 gpx=px;
2502 gpy=py;
2503 gpz=pz;
2504 gp=p;
2505 if(indexOfTestTrack) {
2506 gtheta=theta0;
2507 EpMulScatt(gpx,gpy,gpz,gp,gtheta);
2508
2509 } else {
2510 gtheta=0;
2511 }
2512
2513 weight=TMath::Exp(-gtheta*gtheta/(2*theta0*theta0));
2514 sumOfTheta += gtheta;
2515
2516 // ==========================================================
2517 // Calculate crossing of the track in magnetic field with cylidrical surface
2518 // of radius RTOFINNER
2519 // chargeField = qB, where q is a charge of a particle in units of e,
2520 // B is magnetic field in tesla
2521 // see 3.3.1.1. in the book "Data analysis techniques for
2522 // high-energy physics experiments", edited by M.Regler
2523 // in Russian: "Metody analiza dannykh v fizicheskom eksperimente"
2524 // Moskva, "Mir", 1993. ctr.306
2525
2526 // Initial constants
2527 rotationFactor=1.;
2528 if(chargeField<0.) rotationFactor=-1.;
2529 rotationFactor=-rotationFactor;
2530 gpt=gpx;
2531 phi0=gpy;
2532 cylcor(gpt,phi0);
2533 phi0 -= rotationFactor*TMath::Pi()*0.5;
2534 // phi0 -= h*PID2;
2535 coslam=gpt/gp;
2536 sinlam=gpz/gp;
2537 // helixRadius=100.*gpt/TMath::Abs(0.299792458*chargeField);
2538 helixRadius=100.*gpt/TMath::Abs(AliTOFConstants::fgkSpeedOfLight*chargeField);
2539 xHelixCenter=x-helixRadius*TMath::Cos(phi0);
2540 yHelixCenter=y-helixRadius*TMath::Sin(phi0);
2541 zHelixCenter=z;
2542 helixFactor=rotationFactor*coslam/helixRadius;
2543
2544 // Solves the equation f(s)=r(s)-RTOFINNER=0 by the Newton's method:
2545 // snew=s-f/f'
2546 istep=0;
2547 s=AliTOFConstants::fgkrmin-TMath::Sqrt(x*x+y*y);;
2548 do {
2549 istep++;
2550 xRespectToHelixCenter=helixRadius*TMath::Cos(phi0+s*helixFactor);
2551 yRespectToHelixCenter=helixRadius*TMath::Sin(phi0+s*helixFactor);
2552 gx=xHelixCenter+xRespectToHelixCenter;
2553 gy=yHelixCenter+yRespectToHelixCenter;
2554 gr=TMath::Sqrt(gx*gx+gy*gy);
2555 deltaRadius=gr-AliTOFConstants::fgkrmin;
2556 xp=-helixFactor*yRespectToHelixCenter;
2557 yp= helixFactor*xRespectToHelixCenter;
2558 fp=(gx*xp+gy*yp)/gr;
2559 ds=deltaRadius/fp;
2560 s -= ds;
2561 if(istep==20) {
2562 istep=0;
2563 break;
2564 }
2565 } while (TMath::Abs(ds)>0.01);
2566
2567
2568 if(istep==0) goto end;
2569
2570 // Steps along the circle till a pad
2571 wPixel=0;
2572 wLength=0.;
2573 iplate=0;
2574 iPadAlongX=0;
2575 grho=0.;
2576 ds=fStep;
2577 gxLast=xHelixCenter+helixRadius*TMath::Cos(phi0+s*helixFactor);
2578 gyLast=yHelixCenter+helixRadius*TMath::Sin(phi0+s*helixFactor);
2579 gzLast=zHelixCenter+s*sinlam;
2580
2581
2582 do {
2583 istep++;
2584 s += ds;
2585 gx=xHelixCenter+helixRadius*TMath::Cos(phi0+s*helixFactor);
2586 gy=yHelixCenter+helixRadius*TMath::Sin(phi0+s*helixFactor);
2587 gz=zHelixCenter+s*sinlam;
2588 rho=TMath::Sqrt(gx*gx+gy*gy);
2589
b9d0a01d 2590 IsInsideThePad(gMC,gx,gy,gz,npixel,zPad,xPad);
db9ba97f 2591
2592 iplate += npixel[1];
2593 iPadAlongX += npixel[4];
2594
2595 if(indexOfTestTrack==0 && iplate && iPlateFirst==0) {
2596 iPlateFirst=1;
2597 length0=s;
2598
2599 radius=s*3*theta0;
2600 area=TMath::Pi()*radius*radius;
2601 normR=TMath::Sqrt(gx*gx+gy*gy);
2602 normS=TMath::Sqrt((gx-gxLast)*(gx-gxLast)+
2603 (gy-gyLast)*(gy-gyLast)+
2604 (gz-gzLast)*(gz-gzLast));
2605
2606 cosAngl=(gx*(gx-gxLast)+gy*(gy-gyLast))/(normR*normS);
2607 if(cosAngl<0) printf(" cosAngl<0: gx=%f,gy=%f, gxLast=%f,gyLast=%f,gzLast=%f\n",gx,gy,gxLast,gyLast,gzLast);
2608
2609 area /= cosAngl;
2610 TestTracks=(Int_t) (2*area/(AliTOFConstants::fgkXPad * AliTOFConstants::fgkZPad));
2611
2612 if(TestTracks<12) TestTracks=12;
2613
2614 // Angles of entering into the TOF plate
2615
2616 Int_t iZ=0;
2617 if(TMath::Abs(gz)>300) {
2618 iZ=4;
2619 } else if(TMath::Abs(gz)>200) {
2620 iZ=3;
2621 } else if(TMath::Abs(gz)>100) {
2622 iZ=2;
2623 } else if(TMath::Abs(gz)>0) {
2624 iZ=1;
2625 }
2626
2627
2628 } // end of if(indexOfTestTrack==0 && iplate && iPlateFirst==0)
2629
2630
2631 if(npixel[4]>0) {
2632
2633 iwork1=npixel[0];
2634 iwork2=npixel[1];
2635 iwork3=npixel[2];
2636 // iwork4=npixel[3];
2637 iwork4=(npixel[3]-1)*AliTOFConstants::fgkNpadX+npixel[4];
2638
2639 Int_t ifirstindex=AliTOFConstants::fgkNSectors*(npixel[1]-1)+npixel[0];
2640 iMapValue=mapPixels[ifirstindex-1][iwork3-1][iwork4-1];
2641 if(iMapValue==0) {
2642 ipixel++;
2643 if(ipixel>fMaxPixels) {
2644 cout << "ipixel=" << ipixel << " > MAXPIXELS=" << fMaxPixels << endl;
2645 break;
2646 }
2647 mapPixels[ifirstindex-1][iwork3-1][iwork4-1]=ipixel;
2648 pixelArray[ipixel-1].SetGeom(iwork1,iwork2,iwork3,iwork4);
2649 iMapValue=ipixel;
2650 }
2651
2652 wPixel=iMapValue;
2653 wLength=tpclength+s;
2654 wRho=rho;
2655 wZ=gz;
2656
2657 ihit=kTOFhitFirst[ipart];
2658
2659 if(ihit) {
2660 if(indexOfTestTrack==0) {
2661 {
2662 idelR++;
2663 delR=TMath::Sqrt((gx-hitArray[ihit-1].X())*(gx-hitArray[ihit-1].X())+
2664 (gy-hitArray[ihit-1].Y())*(gy-hitArray[ihit-1].Y())+
2665 (gz-hitArray[ihit-1].Z())*(gz-hitArray[ihit-1].Z()));
2666
2667 }
2668
2669 if(delR>hitArray[ihit-1].GetRmin()) iRmin++;
2670 gxExtrap=gx;
2671 gyExtrap=gy;
2672 gzExtrap=gz;
2673 } else {
2674 delR=TMath::Sqrt((gx-gxExtrap)*(gx-gxExtrap)+
2675 (gy-gyExtrap)*(gy-gyExtrap)+
2676 (gz-gzExtrap)*(gz-gzExtrap));
2677 }
2678 } //end of if(ihit)
2679
2680 break;
2681
2682 } //end of npixel[4]
2683
2684 if(rho<grho) {
2685 istep=0;
2686 break;
2687 }
2688 grho=rho;
2689
2690 gxLast=gx;
2691 gyLast=gy;
2692 gzLast=gz;
2693
2694 } while(rho<AliTOFConstants::fgkrmax); //end of do
2695
2696
2697 if(istep>0) {
2698 if(iplate) {
2699 if(iPadAlongX==0) {
2700 istep=-3; // holes in TOF
2701 }
2702 } else {
2703 if(TMath::Abs(gz)<AliTOFConstants::fgkMaxhZtof) {
2704 // if(TMath::Abs(gz)<MAXZTOF2) {
2705 istep=-2; // PHOS and RICH holes or holes in between TOF plates
2706 } else {
2707 istep=-1; // out of TOF on z-size
2708 }
2709 }
2710 }
2711
2712 if(iPadAlongX>0) {
2713 if(itest==0) {
2714 itest=1;
2715 ntest[itest-1]=1;
2716 testPixel[itest-1]=wPixel;
2717 testLength[itest-1]=wLength;
2718 testRho[itest-1]=wRho;
2719 testZ[itest-1]=wZ;
2720 testWeight[itest-1]=weight;
2721 } else {
be167bf4 2722 Int_t k=0;
db9ba97f 2723 for(Int_t i=0;i<itest;i++) {
2724 k=0;
2725 if(testPixel[i]==wPixel) {
2726 k=1;
2727 ntest[i]++;
2728 testLength[i] += wLength;
2729 testRho[i] += wRho;
2730 testZ[i] += wZ;
2731 testWeight[i] += weight;
2732 break;
2733 }
2734 } //end for i
2735 if(k==0) {
2736 itest++;
2737 ntest[itest-1]=1;
2738 testPixel[itest-1]=wPixel;
2739 testLength[itest-1]=wLength;
2740 testRho[itest-1]=wRho;
2741 testZ[itest-1]=wZ;
2742 testWeight[itest-1]=weight;
2743 }
2744 }
2745 }
2746
2747 end: ;
2748 // Statistics
2749 if(fMatchingStyle==1) {
2750 if(istep>-4 && istep<1) weightTestTracksOutTof[-istep] ++;
2751 } else {
2752 if(istep>-4 && istep<1) weightTestTracksOutTof[-istep] += weight;
2753 }
2754
2755 if(fMatchingStyle==2) {
2756 if(indexOfTestTrack==0 && istep==0) break;
2757 if(indexOfTestTrack+1==TestTracks) break;
2758 }
2759
2760 } //end of indexOfTestTrack
2761
2762 snr += (Float_t) (indexOfTestTrack+1);
2763
2764 // Search for the "hole" with the largest weigth
2765 largestWeightTracksOutTof=0.;
2766 sumWeightTracksOutTof=0.;
2767 for(Int_t i=0;i<4;i++) {
2768 w=weightTestTracksOutTof[i];
2769 sumWeightTracksOutTof += w;
2770 if(w>largestWeightTracksOutTof) {
2771 largestWeightTracksOutTof=w;
2772 iMaxTestTracksOutTof=i;
2773 }
2774 }
2775
2776 itestc=itest;
2777 if(itest>0) {
2778 for(Int_t i=0;i<itest;i++) {
2779 testLength[i] /= ntest[i];
2780 testRho[i] /= ntest[i];
2781 testZ[i] /= ntest[i];
2782 }
2783 // Search for the pixel with the largest weigth
2784 wmax=0.;
2785 wstate=0;
2786 sw=0;
2787 sGeomWeigth=0;
2788 for(Int_t i=0;i<itest;i++) {
2789 istate=pixelArray[testPixel[i]-1].GetState();
2790 fstate=0;
2791 if(istate>0) {
2792 fstate=1;
2793 wstate++;
2794 }
2795 if(fMatchingStyle==1) {
2796 sGeomWeigth += ntest[i];
2797 w=(fpadefficiency*fstate+(1.-fpadefficiency)*(1-fstate))*ntest[i];
2798 if(pixelArray[testPixel[i]-1].GetTrackMatched()>0) w *= 0.1;
2799 } else {
2800 sGeomWeigth += testWeight[i];
2801 w=(fpadefficiency*fstate+(1.-fpadefficiency)*(1-fstate))*testWeight[i];
2802 if(pixelArray[testPixel[i]-1].GetTrackMatched()>0) w *= 0.1;
2803 }
2804
2805 // weighting according to the Pulse Height (we use the square of weight)
2806 // if (fChargeFactorForMatching) w *= (pixelArray[testPixel[i]-1].GetCharge())*(pixelArray[testPixel[i]-1].GetCharge());
2807 if (fChargeFactorForMatching && fstate==1) w *= (pixelArray[testPixel[i]-1].GetCharge())*(pixelArray[testPixel[i]-1].GetCharge());
2808
2809 if(w>wmax) {
2810 wmax=w;
2811 imax=i;
2812 }
2813 sw += w;
2814 }
2815 wPixel=testPixel[imax];
2816 wLength=testLength[imax];
2817 istate=pixelArray[wPixel-1].GetState();
2818
2819 //Choose the TOF dead space
2820 // if(istate==0 && largestWeightTracksOutTof>wmax) {
2821 // if(istate==0 && largestWeightTracksOutTof>=sw) {
2822 if(istate==0 && sumWeightTracksOutTof>sGeomWeigth) {
2823 itestc=itest;
2824 itest=0;
2825 }
2826 }
2827
2828 if(itest>0) {
2829
2830 // Set for MyTrack: Pixel
2831 trackArray[itrack-1].SetPixel(wPixel);
2832
2833 istate=pixelArray[wPixel-1].GetState();
2834
2835 if(istate) {
2836
2837 // Set for MyTrack: Pixel, Length, TOF, MassTOF
2838 //fp
2839 //time=pixelArray[wPixel-1].GetTime();
2840 time=pixelArray[wPixel-1].GetRealTime();
2841 trackArray[itrack-1].SetLength(wLength);
2842 trackArray[itrack-1].SetTof(time);
2843
2844 inverseOfParticleSpeed=time/wLength;
2845 //w=900.*inverseOfParticleSpeed*inverseOfParticleSpeed-1.;
2846 w=(100.*AliTOFConstants::fgkSpeedOfLight)*(100.*AliTOFConstants::fgkSpeedOfLight)*inverseOfParticleSpeed*inverseOfParticleSpeed-1.;
2847 w2=pvtx*pvtx;
2848 Float_t squareMass=w2*w;
2849 mass=TMath::Sqrt(TMath::Abs(squareMass));
2850 if(w<0.) mass=-mass;
2851
2852 trackArray[itrack-1].SetMassTOF(mass);
2853
2854 // Set for MyTrack: Matching
2855 match=4;
2856 // if(ipart==pixelArray[wPixel-1].GetTrack()) match=3;
2857 if( (ipart==pixelArray[wPixel-1].GetTrack()) && hitArray[pixelArray[wPixel-1].GetHit()-1].GetNoise()==0)match=3;
2858 imatched=pixelArray[wPixel-1].GetTrackMatched();
2859 // Set for TOFPixel the number of matched track
2860 pixelArray[wPixel-1].SetTrackMatched(itrack);
2861
2862 if(imatched>0) {
2863 matchw=trackArray[imatched-1].GetMatching();
2864 if(match==3 && matchw==4) t34++;
2865 if(match==3 && matchw==2) t32++;
2866 if(match==4 && matchw==4) t44++;
2867 if(match==4 && matchw==3) t43++;
2868 if(match==4 && matchw==2) t42++;
2869 if(iTOFpixel[ipart]==0 || iTOFpixel[trackArray[imatched-1].GetTrack()]==0) {
2870 m20++;
2871 } else if(iTOFpixel[ipart]==iTOFpixel[trackArray[imatched-1].GetTrack()]) {
2872 m22++;
2873 } else {
2874 m23++;
2875 wPix=iTOFpixel[ipart];
2876 if(PRINT && iPrintM1==10 && iPrintM2<10) {
2877 if(iPrintM2==0) {
2878 printf("*** test print for tracks matched with the pixel for with we had matched track\n");
2879 }
2880 iPrintM2++;
2881 printf(" m=2: ipart=%i, pdgCode=%i, p=%f, theta0=%f, %i Pixel(LP=%i,SP=%i,P=%i) \n",
2882 ipart,pdgCode,p,theta0,wPix,
2883 pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel());
2884 printf(" mat=%i, %i Pixel(LP=%i,SP=%i,P=%i), Test(n=%i,i=%i,w=%f,z=%f), wst=%i \n",
2885 match,wPixel,
2886 pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetPixel(),
2887 itest,imax,wmax,testZ[imax],wstate);
2888 Int_t fstat,istat;
2889 for(Int_t i=0;i<itest;i++) {
2890 wPix=testPixel[i];
2891 istat=pixelArray[wPix-1].GetState();
2892 fstat=0;
2893 if(istat>0) fstat=1;
2894 w=(fpadefficiency*fstat+(1.-fpadefficiency)*(1-fstat))*ntest[i];
2895 if(istat>0)
2896 printf(" %i: %i Pixel(LP=%i,SP=%i,P=%i), istat=%i, ntest=%i, w=%f\n",i+1,
2897 wPix,pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel(),
2898 istat,ntest[i],w);
2899 }
2900 printf(" mat=%i, %i Pixel \n",matchw,trackArray[imatched-1].GetPad());
2901 }
2902 }
2903 if(wstate>1) m2state++;
2904 smat[matchw+4]--;
2905 match=2;
2906 trackArray[imatched-1].SetMatching(match);
2907 smat[match+4]++;
2908
2909 } // if(imatched>0)
2910
2911 } else { //else if(istate)
2912
2913 match=1;
2914 if(iTOFpixel[ipart]==0) m10++;
2915 if(PRINT && iPrintM1<10) {
2916 Int_t wPix;
2917 wPix=iTOFpixel[ipart];
2918 if(wPix) {
2919 if(iPrintM1==0) {
2920 printf("*** test print for tracks fired a pixel but matched with non-fired pixel\n");
2921 }
2922 iPrintM1++;
2923 printf(" m=1: itra=%i,ipart=%i, pdgCode=%i, p=%f, theta0=%f, %i Pixel(LP=%i,SP=%i,P=%i) \n",
2924 itra,ipart,pdgCode,p,theta0,wPix,
2925 pixelArray[wPix-1].GetSector(),pixelArray[wPix-1].GetPlate(),pixelArray[wPix-1].GetPixel());
2926 printf(" mat=%i, %i Pixel(LP=%i,SP=%i,P=%i), Test(n=%i,i=%i,w=%f,z=%f), wst=%i \n",
2927 match,wPixel,
2928 pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetPixel(),
2929 itest,imax,wmax,testZ[imax],wstate);
2930
2931 }
2932 } //end if(PRINT && iPrintM1<10)
2933
2934 } //end if(istate)
2935
2936 } else {
2937 match=-1-iMaxTestTracksOutTof;
2938
2939 } //end itest
2940
2941 trackArray[itrack-1].SetMatching(match);
2942 // if(iTestGmax==1) hMTT->Fill(match);
2943 smat[match+4]++;
2944
2945 sumOfTheta /= iTestTrack;
2946
2947 itest=itestc;
2948
2949 //Test
2950 if(PRINT) {
2951 if(iTOFpixel[ipart] && match!=3) {
2952 particle = (TParticle*)gAlice->Particle(ipart); //for V3.05
2953
be167bf4 2954 printf(" ipixel=%i (Sector=%i, Plate=%i, Strip=%i, Pixel=%i), fired by %i track\n",iTOFpixel[ipart],pixelArray[iTOFpixel[ipart]-1].GetSector(),pixelArray[iTOFpixel[ipart]-1].GetPlate(),pixelArray[iTOFpixel[ipart]-1].GetStrip(),pixelArray[iTOFpix
2955el[ipart]-1].GetPixel(),pixelArray[iTOFpixel[ipart]-1].GetTrack());
2956 printf(" indexOfTestTrack=%i itest=%i weightTestTracksOutTof[4]=%f weightTestTracksOutTof[2]=%f weightTestTracksOutTof[1]=%f weightTestTracksOutTof[0]=%f\n",indexOfTestTrack,itest,weightTestTracksOutTof[3],weightTestTracksOutTof[2],weightTestTracksOu
2957tTof[1],weightTestTracksOutTof[0]);
db9ba97f 2958 if(itest) {
2959
be167bf4 2960 printf(" take ipixel=%i (Sector=%i, Plate=%i, Strip=%i, Pixel=%i), (fired by %i track), match=%i\n",wPixel,pixelArray[wPixel-1].GetSector(),pixelArray[wPixel-1].GetPlate(),pixelArray[wPixel-1].GetStrip(),pixelArray[wPixel-1].GetPixel(),pixelArray[w
2961Pixel-1].GetTrack(),match);
db9ba97f 2962 }
2963 }
2964 }
2965 if(PRINT && itra<10 ) {
2966
2967 if(itest) {
2968 cout << " number of pixels with test tracks=" << itest << endl;
2969 for(Int_t i=0;i<itest;i++) {
2970 cout << " " << i+1 << " tr.=" << ntest[i] << " w=" << testWeight[i] << " pix.= " << testPixel[i] << " (" <<
2971 pixelArray[testPixel[i]-1].GetSector() << " " << " " << pixelArray[testPixel[i]-1].GetPlate() << " " <<
2972 pixelArray[testPixel[i]-1].GetPixel() << " )" << " l= " << testLength[i] << " sig=" <<
2973 theta0*(testLength[i]-tpclength) << " rho= " << testRho[i] << " z= " << testZ[i] << endl;
2974 }
2975 cout << " pixel=" << wPixel << " state=" << istate << " l=" << wLength << " TOF=" << time << " m=" << mass << " match=" << match << endl;
2976 if(istate>0) cout << " fired by track " << pixelArray[wPixel-1].GetTrack() << endl;
2977 }
2978 }
2979 } //end of track
2980
2981
2982 if(itr) {
2983 printf(" %f probe tracks per 1 real track\n",snr/itr);
2984 itrack=itr;
2985 }
2986
2987
2988 cout << ipixel << " - total number of TOF pixels after matching" << endl;
2989 w=iRmin;
2990 if(idelR!=0) {
2991 w /= idelR;
2992 printf(" %i tracks with delR, %f of them have delR>Rmin \n",idelR,w);
2993 }
2994 w=iRmin1;
2995 if(idelR1!=0) {
2996 w /= idelR1;
2997 printf(" %i tracks with delR1 (|z|<175), %f of them have delR>Rmin \n",idelR1,w);
2998 }
2999 w=iRmin2;
3000 if(idelR2!=0) {
3001 w /= idelR2;
3002 printf(" %i tracks with delR2 (|z|>175), %f of them have delR>Rmin \n",idelR2,w);
3003 }
3004
3005 cout << " ******************** End of matching **********" << endl;
f9a28264 3006 delete [] ntest;
3007 delete [] testPixel;
3008 delete [] testLength;
3009 delete [] testRho;
3010 delete [] testZ;
3011 delete [] testWeight;
db9ba97f 3012}
3013
3014//____________________________________________________________________________
3015void AliTOFReconstructioner::FillNtuple(Int_t ntracks, AliTOFTrack* trackArray, AliTOFRecHit* hitArray, AliTOFPad* pixelArray, Int_t* iTOFpixel, Int_t* iparticle, Float_t* toftime, Int_t& ipixelLastEntry, Int_t itrack){
3016
3017 // itrack : total number of TPC selected tracks
3018 // for the caller is ntotTPCtracks
3019
3020 cout << " ******************** Start of searching non-matched fired pixels **********" << endl;
3021 const Int_t charge[48]={ 0, 1,-1, 0, 1,-1, 0, 1,-1, 0,
3022 1,-1, 0, 1,-1, 0, 0, 0, 1, 0,
3023 -1, 0,-1,-1, 0, 0,-1, 0, 1, 0,
3024 1, 1, 0, 0, 1,-1, 0, 0, 1,-1,
3025 1, 1,-1, 0, 1, 1, 2, 0};
3026
3027 Int_t macthm1=0;
3028 Int_t macthm2=0;
3029 Int_t macthm3=0;
3030 Int_t macthm4=0;
3031 Int_t macth0=0;
3032 Int_t macth1=0;
3033 Int_t macth2=0;
3034 Int_t macth3=0;
3035 Int_t macth4=0;
3036
3037
3038 Float_t smat[9],smat0[9],smat1[9];
3039 for(Int_t i=0;i<9;i++) {
3040 smat[i]=0.;
3041 smat0[i]=0.;
3042 smat1[i]=0.;
3043 }
3044
3045 Int_t nFiredPixelsNotMatchedWithTracks=0;
3046 Int_t istate;
3047 for (Int_t i=0; i<ipixelLastEntry; i++) {
3048 istate=pixelArray[i].GetState();
3049 if(istate==0) break;
3050 if(pixelArray[i].GetTrackMatched()==-1) nFiredPixelsNotMatchedWithTracks++;
3051 }
3052 printf(" %i fired pixels have not matched tracks\n",nFiredPixelsNotMatchedWithTracks);
3053 cout << " ******************** End of searching non-matched fired pixels **********" << endl;
3054
3055 Int_t nTPCHitMissing=0;
3056 for(Int_t i=0; i<ipixelLastEntry; i++) {
3057 if(pixelArray[i].GetHit()>0) {
3058 if(hitArray[pixelArray[i].GetHit()-1].GetNoise()==0) {
3059 if(iparticle[pixelArray[i].GetTrack()]==0) nTPCHitMissing++;
3060 }
3061 }
3062 }
3063 printf(" %i pixels fired by track hit without a hit on the last layer of TPC\n",nTPCHitMissing);
3064
3065
3066 Int_t icharge=0; // total number of charged particles
3067 Int_t iprim=0; // number of primaries
3068 Int_t ipions=0; // number of primary pions
3069 Int_t ikaons=0; // number of primary kaons
3070 Int_t iprotons=0; // number of primary protons
3071 Int_t ielectrons=0;// number of primary electrons
3072 Int_t imuons=0; // number of primary muons
3073 Float_t particleTypeArray[6][5][2];
3074
3075 for (Int_t index1=0;index1<6;index1++) {
3076 for (Int_t index2=0;index2<5;index2++) {
3077 for (Int_t index3=0;index3<2;index3++) {
3078 particleTypeArray[index1][index2][index3]=0.;
3079 }
3080 }
3081 }
3082
3083 Int_t nTOFhitsWithNoTPCTracks=0; // to be moved later when used
3084
3085 /*
3086 TObjArray *Particles = gAlice->Particles();
3087 Int_t numberOfParticles=Particles->GetEntries();
3088 cout << "numberOfParticles " << numberOfParticles << endl;
3089 // fpdbg
3090 if(numberOfParticles>fMaxAllTracks) numberOfParticles=fMaxAllTracks;
3091 */
3092
3093 for (Int_t i=0; i<ntracks; i++) { // starting loop on all primaries charged particles for current event)
3094
3095 /*
3096 cout << "particle " << i << endl;
3097 cout << "total " << numberOfParticles << endl;
3098 */
3099 TParticle *part = (TParticle *) gAlice->Particle(i);
3100 if(charge[PDGtoGeantCode(part->GetPdgCode())-1]) {
3101 icharge++;
3102 /*
3103 cout << "charged particles " << icharge << endl;
3104 */
3105 Int_t particleType=0;
3106 Int_t absPdgCode = TMath::Abs(part->GetPdgCode());
3107 switch (absPdgCode) {
3108 case 211:
3109 particleType=3;
3110 break ;
3111 case 321:
3112 particleType=2;
3113 break ;
3114 case 2212:
3115 particleType=1;
3116 break ;
3117 case 11:
3118 particleType=4;
3119 break ;
3120 case 13:
3121 particleType=5;
3122 break ;
3123 }
3124
3125 if(part->GetFirstMother() < 0) {
3126 iprim++;
3127 switch (particleType) {
3128 case 1:
3129 iprotons++;
3130 break ;
3131 case 2:
3132 ikaons++;
3133 break ;
3134 case 3:
3135 ipions++;
3136 break ;
3137 case 4:
3138 ielectrons++;
3139 break ;
3140 case 5:
3141 imuons++;
3142 break ;
3143 }
3144 }
3145
3146 Int_t match=0;
3147 Float_t wLength=-1.;
3148 Float_t time=-1.;
3149 Float_t mass=-1.;
3150
3151 Int_t itr=iparticle[i]; // get the track number for the current charged particle
3152
3153 if(iTOFpixel[i]>0 && itr==0) nTOFhitsWithNoTPCTracks++;
3154
3155 if(itr) {
3156 match=trackArray[itr-1].GetMatching();
3157 //cout << "match " << match << endl;
3158 wLength=trackArray[itr-1].GetLength();
3159 //cout << "wLength " << wLength << endl;
3160 time=trackArray[itr-1].GetTof();
3161 mass=trackArray[itr-1].GetMassTOF();
3162 //cout << "mext " << mass << endl;
3163 // if(PRINT && (i>789 && i<800) ) cout << i << " track: l=" << wLength << " TOF=" << time << " m=" << mass << " match=" << match << endl;
3164 if(iTOFpixel[i]==0) {
3165 smat0[match+4]++;
3166 wLength=-wLength;
3167 }
3168 }
3169 Int_t ikparen=part->GetFirstMother();
3170 Int_t imam;
3171 if(ikparen<0) {
3172 imam=0;
3173 } else {
3174 imam=part->GetPdgCode();
3175 }
3176
3177 Int_t evnumber=gAlice->GetEvNumber();
3178 if(match==-1) macthm1++;
3179 if(match==-2) macthm2++;
3180 if(match==-3) macthm3++;
3181 if(match==-4) macthm4++;
3182 if(match==0) macth0++;
3183 if(match==1) macth1++;
3184 if(match==2) macth2++;
3185 if(match==3) macth3++;
3186 if(match==4) macth4++;
3187 foutputntuple->Fill(evnumber,part->GetPdgCode(),imam,part->Vx(),part->Vy(),part->Vz(),part->Px(),part->Py(),part->Pz(),toftime[i],wLength,match,time,mass);
3188
3189
3190
3191 // -----------------------------------------------------------
3192 // Filling 2 dimensional Histograms true time vs matched time
3193 // Filling 1 dimensional Histogram true time - matched time
3194 //
3195 // time = time associated to the matched pad [ns]
3196 // it could be the average time of the cluster fired
3197 //
3198 // toftime[i] = real time (including pulse height delays) [s]
3199 //
3200 //
3201 // if (time>=0) {
3202 // if (imam==0) TimeTrueMatched->Fill(time, toftime[i]*1E+09);
3203 // if (imam==0) DeltaTrueTimeMatched->Fill(time-toftime[i]*1E+09);
3204 // }
3205 //
3206 //---------------------------------------------------------------
3207
3208 if(match==-4 || match>0) {
3209 Int_t matchW;
3210 matchW=match;
3211 if(match==-4) matchW=1;
3212 if(particleType) {
3213 particleTypeArray[particleType-1][matchW-1][1]++;
3214 particleTypeArray[5][matchW-1][1]++;
3215 particleTypeArray[particleType-1][4][1]++;
3216 particleTypeArray[5][4][1]++;
3217 if(part->GetFirstMother() < 0) {
3218 particleTypeArray[particleType-1][matchW-1][0]++;
3219 particleTypeArray[5][matchW-1][0]++;
3220 particleTypeArray[particleType-1][4][0]++;
3221 particleTypeArray[5][4][0]++;
3222
3223 // fill histos for QA
3224 //if(particleType==3 && matchW==3) hPiWithTrueTime->Fill(sqrt((part->Px())*(part->Px())+(part->Py())*(part->Py())+(part->Pz())*(part->Pz())));
3225 //if(particleType==2 && matchW==3) hKWithTrueTime->Fill(sqrt((part->Px())*(part->Px())+(part->Py())*(part->Py())+(part->Pz())*(part->Pz())));
3226 //if(particleType==1 && matchW==3) hPWithTrueTime->Fill(sqrt((part->Px())*(part->Px())+(part->Py())*(part->Py())+(part->Pz())*(part->Pz())));
3227 //
3228
3229 } // close if(part->GetFirstMother() < 0)
3230 } // close if(particleType)
3231 } // close if(match==-4 || match>0)
3232 } // close if(charge[PDGtoGeantCode(part->GetPdgCode())-1])
3233 } // close for (Int_t i=0; i<ntracks; i++) {
3234
3235 cout << " macthm1 " << macthm1 << endl;
3236 cout << " macthm2 " << macthm2 << endl;
3237 cout << " macthm3 " << macthm3 << endl;
3238 cout << " macthm4 " << macthm4 << endl;
3239 cout << " macth0 " << macth0 << endl;
3240 cout << " macth1 " << macth1 << endl;
3241 cout << " macth2 " << macth2 << endl;
3242 cout << " macth3 " << macth3 << endl;
3243 cout << " macth4 " << macth4 << endl;
3244
3245
3246 printf(" %i TOF hits have not TPC track\n",nTOFhitsWithNoTPCTracks);
3247 Int_t imatch=0;
3248 for(Int_t i=0;i<9;i++) {
3249 if(itrack) cout << " " << smat[i]*100./itrack << " % of them (="<<smat[i]<<") have match=" << i-4 << " " << smat0[i] << " have not TOF hits" << endl;
3250 if(i==0 || i>4) imatch += (Int_t) (smat[i]);
3251
3252 // cout << " " << smat[i]*100./itrack << " % of them (="<<smat[i]<<") have match=" << i-4 << " " << smat0[i] << " have not TOF hits" << " " << smat1[i] << " have (r.p)<0 for first hit" << endl;
3253 }
3254
3255 if(fdbg){
3256 /*
3257 cout << " nparticles = " << numberOfParticles << " charged = " << icharge << " prim.=" << iprim << endl;
3258 */
3259 cout << " nparticles = " << ntracks << " charged = " << icharge << " prim.=" << iprim << endl;
3260 cout << ipions << " - primary pions" << endl;
3261 cout << ikaons << " - primary kaons" << endl;
3262 cout << iprotons << " - primary protons" << endl;
3263 cout << ielectrons << " - primary electrons" << endl;
3264 cout << imuons << " - primary muons reached TPC" << endl;
3265 cout << " ********** " << imatch << " TPC tracks are matched with TOF pixels (incl.match=-4) **********" << endl;
3266 }
3267
3268 /*
3269 Float_t PrimaryInBarrel[6],Acceptance[6];
3270 PrimaryInBarrel[0]=ipions;
3271 PrimaryInBarrel[1]=ikaons;
3272 PrimaryInBarrel[2]=iprotons;
3273 PrimaryInBarrel[3]=ielectrons;
3274 PrimaryInBarrel[4]=imuons;
3275 PrimaryInBarrel[5]=ipions+ikaons+iprotons+ielectrons+imuons;
3276
3277 // cout << " TPC acceptance for the primary species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl;
3278 for(Int_t i=0; i<6; i++) {
3279 Acceptance[i]=0.;
3280 if(PrimaryInBarrel[i]) Acceptance[i]=100.*PrimaryReachedTPC[i]/PrimaryInBarrel[i];
3281 //hTPCacceptance[i]->Fill(Acceptance[i]);
3282 // printf(" species: %i %f\n",i+1,Acceptance[i]);
3283 }
3284
3285 // cout << " TOF acceptance for the primary species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl;
3286 for(Int_t i=0; i<6; i++) {
3287 Acceptance[i]=0.;
3288 if(PrimaryInBarrel[i]) Acceptance[i]=100.*PrimaryReachedTOF[i]/PrimaryInBarrel[i];
3289 //hTOFacceptance[i]->Fill(Acceptance[i]);
3290 // printf(" species: %i %f\n",i+1,Acceptance[i]);
3291 }
3292
3293 for (Int_t index1=0;index1<6;index1++) {
3294 for (Int_t index2=0;index2<4;index2++) {
3295 for (Int_t index3=0;index3<2;index3++) {
3296 if(particleTypeArray[index1][4][index3]) particleTypeArray[index1][index2][index3]=
3297 100.*particleTypeArray[index1][index2][index3]/particleTypeArray[index1][4][index3];
3298 }
3299 }
3300 }
3301
3302 cout << "species: 1-p, 2-K, 3-pi, 4-e, 5-mu, 6-all" << endl;
3303 cout << " matched pixels(%): 1-unfired 2-double 3-true 4-wrong 5-total number of tracks" << endl;
3304
3305 cout << " primary tracks:" << endl;
3306 for (Int_t i=0;i<6;i++) {
3307 cout << i+1 << " " << particleTypeArray[i][0][0] << " " << particleTypeArray[i][1][0] << " " << particleTypeArray[i][2][0] << " " << particleTypeArray[i][3][0] << " " << particleTypeArray[i][4][0] << endl;
3308 }
3309
3310 // cout<<" contam.for all prim.(%)="<<100*particleTypeArray[5][3][0]/(particleTypeArray[5][3][0]+particleTypeArray[5][2][0])<<endl;
3311
3312 cout << " all tracks:" << endl;
3313 for (Int_t i=0;i<6;i++) {
3314 cout << i+1 << " " << particleTypeArray[i][0][1] << " " << particleTypeArray[i][1][1] << " " << particleTypeArray[i][2][1] << " " << particleTypeArray[i][3][1] << " " << particleTypeArray[i][4][1] << endl;
3315 }
3316
3317 // cout<<" contam.for all (%)="<<100*particleTypeArray[5][3][1]/(particleTypeArray[5][3][1]+particleTypeArray[5][2][1])<<endl;
3318 // printf(" t34=%i, t32=%i, t44=%i, t43=%i, t42=%i\n",t34,t32,t44,t43,t42);
3319 // printf(" m10=%f, m20=%f, m22=%f, m23=%f, m2state=%i\n",m10,m20,m22,m23,m2state);
3320 */
3321}