Fixed Warnings in Reconstructioner and Prob
[u/mrichter/AliRoot.git] / TOF / AliTOFReconstructionerV2.cxx
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 // This is a TTask for reconstruction V2 in TOF
18 // Description of the algorithm
19 //-- Author: F. Pierella | pierella@bo.infn.it
20 //////////////////////////////////////////////////////////////////////////////
21
22 #include "TBenchmark.h"
23 #include "TFile.h"
24 #include "TFolder.h"
25 #include "TParticle.h"
26 #include "TTask.h"
27 #include "TTree.h"
28 #include "TClonesArray.h"
29 #include "TROOT.h"
30 #include "TSystem.h"
31 #include <stdlib.h>
32 #include <iostream.h>
33 #include <fstream.h>
34 #include "TGeant3.h"
35 #include "TVirtualMC.h"
36
37 #include "AliDetector.h"
38 #include "AliMC.h"
39 #include "AliRun.h"
40 #include "AliTOF.h"
41 #include "AliTOFDigitMap.h"
42 #include "AliTOFHitMap.h"
43 #include "AliTOFhitT0.h"
44 #include "AliTOFdigit.h"
45 #include "AliTOFConstants.h"
46 #include "AliTOFReconstructionerV2.h"
47 #include "AliTOFTrackV2.h"
48
49 #include "AliKalmanTrack.h"
50 #include "../TPC/AliTPCtrack.h"
51 #include "../TRD/AliTRDtrack.h"
52
53 ClassImp(AliTOFReconstructionerV2)
54
55 //____________________________________________________________________________ 
56   AliTOFReconstructionerV2::AliTOFReconstructionerV2():TTask("AliTOFReconstructionerV2","") 
57 {
58   //
59   // std ctor
60   // set all member vars to zero
61   //
62   fdbg      =0;
63   fDigitsMap=0x0;
64   fField    =0;
65   fNDummyTracks=0;
66   fScaleSigmaFactor=0.;
67   fStep     =0; 
68   fTOFDigits=0x0;
69   fTOFTracks=0x0;
70   fTOFDigitsFile    ="digits.root";
71   fTPCBackTracksFile="AliTPCBackTracks.root";
72   fKalmanTree       =0x0;
73   fBranchWithTracks =0x0;
74 }
75            
76 //____________________________________________________________________________ 
77 AliTOFReconstructionerV2::AliTOFReconstructionerV2(char* tpcBackTracks, char* tofDigits):TTask("AliTOFReconstructionerV2","") 
78 {
79   //
80   // par ctor
81   // default arguments are specified only in
82   // the header file
83   // Parameters:
84   // tpcBackTracks -> file name with backpropagated tracks in TPC
85   // tofDigits     -> file name with TOF digits
86   //
87
88   fdbg      =0;
89   fDigitsMap=0x0;
90   fField    =0.2;   // default value 0.2 [Tesla]
91   fNDummyTracks=20; // by default 20 test tracks used
92   fScaleSigmaFactor=1.;
93   fStep     =0.01;  // [cm] 
94   fTOFDigits=0x0;
95   fTOFTracks=0x0;
96   fTOFDigitsFile    =tofDigits;
97   fTPCBackTracksFile=tpcBackTracks;
98   fKalmanTree       =0x0;
99   fBranchWithTracks =0x0;
100
101   // initialize the G3 geometry 
102   gAlice->Init();
103   gAlice->Print(); 
104
105   // add Task to //root/Tasks folder
106   TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ; 
107   roottasks->Add(this) ; 
108 }
109
110 //____________________________________________________________________________
111 AliTOFReconstructionerV2::AliTOFReconstructionerV2(const AliTOFReconstructionerV2 & rec)
112 {
113   //
114   // Dummy copy constructor
115   // required by coding conventions
116   ;
117 }
118
119
120
121 //____________________________________________________________________________ 
122   AliTOFReconstructionerV2::~AliTOFReconstructionerV2()
123 {
124   //
125   // dtor
126   // some delete has to be moved
127   // 
128
129   if (fDigitsMap)
130     {
131       delete fDigitsMap;
132       fDigitsMap = 0;
133     }
134
135   if (fTOFDigits)
136     {
137       delete fTOFDigits;
138       fTOFDigits = 0;
139     }
140
141   if (fTOFTracks)
142     {
143       delete fTOFTracks;
144       fTOFTracks = 0;
145     }
146
147   if (fTOFDigitsFile)
148     {
149       delete fTOFDigitsFile;
150     }
151
152   if (fTPCBackTracksFile)
153     {
154       delete fTPCBackTracksFile;
155     }
156
157   if (fKalmanTree)
158     {
159       delete fKalmanTree;
160       fKalmanTree = 0;
161     }
162
163   if (fBranchWithTracks)
164     {
165       delete fBranchWithTracks;
166       fBranchWithTracks = 0;
167     }
168 }
169
170
171 //____________________________________________________________________________
172 void AliTOFReconstructionerV2::Exec(Option_t* option) 
173
174   //
175   // Description of the algorithm:
176   //
177   //
178   //
179   //
180   //
181   //
182
183   // load TOF digits and fill the digit map
184   Int_t tofDigitsLoading=LoadTOFDigits();
185
186   // load back-propagated tracks in TPC
187   Int_t tpcTracksLoading=LoadTPCTracks();
188
189   if(tofDigitsLoading || tpcTracksLoading) {
190     cout<<" Couldn't start reconstruction V2. Exit."<<endl;
191     exit(1);
192   }
193   // or load TRD reconstructed tracks
194   // Int_t trdTracksLoading=LoadTRDTracks();
195
196   // create a TObjArray to store reconstructed tracks
197   // and reject fake tracks
198   const Int_t maxRecTracks = 10000; // max number of reconstructed tracks
199                                     // per event 
200
201   const Float_t stripRegionHeight = 15.3; // [cm] height in radial direction 
202                                           // of the volume where strips are placed
203   TObjArray trackArray(maxRecTracks);
204
205   // buffer for output
206   fTOFTracks= new TClonesArray("AliTOFTrackV2");
207   // create a reference to fill the TClonesArray
208   TClonesArray &aTOFTracks = *fTOFTracks;
209
210   const Int_t maxIndex = 100000; // max number of primary tracks to be analysed 
211   // the index of the rtIndex array is the track label
212   // the content of the array is -1 for fake tracks
213   // and the track index for good tracks
214   Float_t dEdXarray[maxRecTracks];
215   Int_t rtIndex[maxIndex];
216   for(Int_t i = 0; i < maxIndex; i++) rtIndex[i] = -1;
217
218   AliKalmanTrack::SetConvConst(100/0.299792458/fField);
219
220   Int_t nRecTracks = (Int_t) fKalmanTree->GetEntries();
221   cout<<"Found "<<nRecTracks<<" entries in the track tree "<<endl;
222
223   // load the tracks into the array
224   for (Int_t i=0; i<nRecTracks; i++) {
225     AliTPCtrack *iotrack=new AliTPCtrack();
226     fBranchWithTracks->SetAddress(&iotrack);
227     fKalmanTree->GetEvent(i);
228     trackArray.AddLast(iotrack);
229     Int_t trackLabel = iotrack->GetLabel();
230     dEdXarray[i]=iotrack->GetdEdx(); // usefull for PID
231     // 
232     // start filling the TClonesArray of AliTOFTrackV2
233     Float_t trdXYZ[3]={0.,0.,0.};
234     Float_t trdPxPyPz[3]={0.,0.,0.};
235
236     // tpc outer wall positions
237     Double_t xk=iotrack->GetX();
238     // get the running coordinates in the lrf
239     // and alpha
240     Double_t x=xk;
241     Double_t y=iotrack->GetY();
242     Double_t z=iotrack->GetZ();
243     Double_t alpha=iotrack->GetAlpha();
244     GetGlobalXYZ(alpha, x, y, z);
245     Float_t tpcXYZ[3]={x,y,z};
246
247     // momentum at the end of TPC
248     Float_t lambda=TMath::ATan(iotrack->GetTgl());
249     Float_t invpt=TMath::Abs(iotrack->Get1Pt());
250     Float_t pt=-99.;
251     if (invpt) pt= 1./invpt; // pt
252     Float_t tpcmom=1./(invpt*TMath::Cos(lambda));
253     Float_t pz=tpcmom*TMath::Sin(lambda);
254     Float_t tpcPtPz[2]={pt,pz};
255
256     Int_t matchingStatus=-1;
257     if(trackLabel < 0) matchingStatus=0;
258     new(aTOFTracks[i]) AliTOFTrackV2(trackLabel,matchingStatus,tpcmom,dEdXarray[i],tpcXYZ,tpcPtPz,trdXYZ,trdPxPyPz);
259     //    printf("rt with %d clusters and label %d \n",
260     //     iotrack->GetNumberOfClusters(), trackLabel);
261
262     if(trackLabel < 0) continue;
263     if(trackLabel >= maxIndex) continue;
264     rtIndex[trackLabel] = i;
265     delete iotrack;
266   }
267
268   if(strstr(option,"MC")) Comparison(rtIndex);
269
270   // start loop on tracks
271   // and backpropagate them from TPC to TOF
272   // backpropagation is performed only
273   // for good tracks (fake tracks rejected)
274   AliTPCtrack *rt;
275
276   for (Int_t i=0; i<nRecTracks; i++) {
277
278     //******* tracking:  extract track coordinates, momentum, etc. 
279     rt = (AliTPCtrack*) trackArray.UncheckedAt(i);
280     // track length to be implemented
281     // Double_t tr_length = rt->GetLength();
282
283     Int_t tpcTrackLabel=rt->GetLabel();
284     // reject fake tracks
285     //if(tpcTrackLabel< 0) continue;
286
287     // starting backpropagation to TOF
288     // here we assume to have backpropagated tracks in TPC
289     // starting point xk=246.055
290     // starting back propagation
291     // outer wall of the TPC
292     Int_t iOuterTPCWall=rt->PropagateTo(261.53,40.,0.06124);
293     // frame with air just to the beginning of the TOF
294     Int_t iFrameWithAir=rt->PropagateTo(370.,36.66,1.2931e-3);
295     // trough the wall of the TOF plate
296     // thickness has changed according to the
297     // last value
298     Int_t iTOFWall=rt->PropagateTo(370.11,24.01,2.7);
299     
300     /*
301       // outer wall of the TPC
302       Int_t iOuterTPCWall=rt->PropagateTo(261.53,40.,0.06124);
303       // air in between TPC and TRD
304       Int_t iFrameWithAir=rt->PropagateTo(294.5,36.66,1.2931e-3);
305       // TRD itself
306       // mean density for the TRD calculated from
307       // TRD Technical Design Report
308       // page 11  -> global thickness
309       // page 23  -> different main layers thickness (Radiator Air/ Drift Chamber Gas /G10)
310       // page 139 -> material budget and radiation lengths
311       Int_t iTRD=rt->PropagateTo(369.1,171.7,0.33);
312       // air in between TRD and TOF
313       Int_t iFrameWithAirbis=rt->PropagateTo(370.,36.66,1.2931e-3);
314       // trough the wall of the TOF plate
315       Int_t iTOFWall=rt->PropagateTo(370.11,24.01,2.7);
316      */
317  
318     // select only cases when
319     // backpropagation succeded
320     // and particle is in the geometrical TOF acceptance along Z
321     // |z|<380 [cm]  
322     AliTOFTrackV2* oTOFtracks=(AliTOFTrackV2*)fTOFTracks->UncheckedAt(i);
323     Bool_t outOfZacceptance=(rt->GetZ()<=380.);
324     if(outOfZacceptance) oTOFtracks->SetMatchingStatus(-2);
325
326     if(iOuterTPCWall==1 && iFrameWithAir==1 && iTOFWall==1 && (!outOfZacceptance)){
327       Double_t cc[15];
328       // get sigmaY and sigmaZ
329       rt->GetExternalCovariance(cc);
330       Double_t sigmaY =TMath::Sqrt(cc[0]); // [cm]
331       Double_t sigmaZ =TMath::Sqrt(cc[2]); // [cm]
332
333       // arrays used by the DigitFinder
334       Int_t nSlot=1;
335       TArrayI *secArray= new TArrayI(nSlot);
336       TArrayI *plaArray= new TArrayI(nSlot);
337       TArrayI *strArray= new TArrayI(nSlot);
338       TArrayI *pdzArray= new TArrayI(nSlot);
339       TArrayI *pdxArray= new TArrayI(nSlot);
340
341       // make fNDummyTracks clones of the current backpropagated track
342       // make a copy of the current track
343       // smear according to the backpropagated area
344       for (Int_t j=0; j<fNDummyTracks; i++) {
345         AliTPCtrack *dummyrt=new AliTPCtrack(*rt);
346         // get ylrf and zlrf
347         Double_t ylrf= dummyrt->GetY();  // P0
348         Double_t zlrf= dummyrt->GetZ();  // P1
349
350         // smear according to sigmaY and sigmaZ
351         //Double_t ylrfNew=gRandom->Gaus(ylrf,fScaleSigmaFactor*sigmaY);
352         //Double_t zlrfNew=gRandom->Gaus(zlrf,fScaleSigmaFactor*sigmaZ);
353
354         // set Y and Z accordingly
355         // setter to be added in the class AliTPCtrack
356         // here I have to modify the AliTPCtrack class
357         // adding the setters for Y and Z
358         //dummyrt->SetY(ylrfNew);
359         //dummyrt->SetZ(zlrfNew);
360
361         // start fine-backpropagation inside the TOF
362         Bool_t padNotFound =kTRUE;
363         Bool_t isInStripsRegion=kTRUE;
364         Double_t xk=dummyrt->GetX();
365
366         while (padNotFound && isInStripsRegion){
367           xk+=fStep;
368           // here we assume a frame with air
369           dummyrt->PropagateTo(xk,36.66,1.2931e-3);
370           // get the running coordinates in the lrf
371           // and alpha
372           Double_t x=xk;
373           Double_t y=dummyrt->GetY();
374           Double_t z=dummyrt->GetZ();
375           Double_t alpha=dummyrt->GetAlpha();
376           GetGlobalXYZ(alpha, x, y, z);
377           
378           // check if the point falls into a pad
379           // using the G3 geometry
380           Int_t* volumeID = new Int_t[AliTOFConstants::fgkmaxtoftree];
381           // volumeID[0] -> TOF Sector range [1-18]
382           // volumeID[1] -> TOF Plate  range [1- 5]
383           // volumeID[2] -> TOF Strip  max range [1-20]
384           // volumeID[3] -> TOF Pad along Z range [1- 2]
385           // volumeID[4] -> TOF Pad along X range [1-48]
386           
387           Float_t zInPadFrame=0.;
388           Float_t xInPadFrame=0.;
389           IsInsideThePad((Float_t)x,(Float_t)y,(Float_t)z, volumeID, zInPadFrame, xInPadFrame);
390           // adding protection versus wrong volume numbering
391           // to be released in the next release after debugging
392           if(volumeID[4]){
393             padNotFound=kFALSE;
394             nSlot++;
395             secArray->Set(nSlot-1);
396             plaArray->Set(nSlot-1);
397             strArray->Set(nSlot-1);
398             pdzArray->Set(nSlot-1);
399             pdxArray->Set(nSlot-1);
400             
401             (*secArray)[nSlot-1]=volumeID[0];
402             (*plaArray)[nSlot-1]=volumeID[1];
403             (*strArray)[nSlot-1]=volumeID[2];
404             (*pdzArray)[nSlot-1]=volumeID[3];
405             (*pdxArray)[nSlot-1]=volumeID[4];
406             
407           } // track falls into a pad volume
408           
409           delete [] volumeID;
410           
411           // check on xk to stop the fine-propagation
412           if(xk>=(370.+stripRegionHeight)) isInStripsRegion=kFALSE;
413           
414         } // close the while for fine-propagation
415         
416         delete dummyrt;
417       } // end loop on test tracks
418       
419       // start TOF digit finder
420       Int_t assignedVol[5]={0,0,0,0,0};
421       Float_t tdc=-1.;
422       Int_t* digitTrackArray=0x0;
423       Bool_t assignedDigit=DigitFinder(secArray, plaArray, strArray, pdzArray, pdxArray, assignedVol, digitTrackArray, tdc);
424
425
426       if(assignedDigit){
427         // fill the tree for tracks with time of flight
428         // tof is given in tdc bin
429         // conversion to [ns]
430         Float_t binWidth=50.; // [ps]
431         Float_t timeOfFlight=tdc*binWidth/1000.;
432
433         // only the first track number contributing
434         // to the assigned digit
435         Int_t tofDigitTrackLabel=digitTrackArray[0];
436
437         // matching status for the current track
438         Int_t matching=4;
439         if(tpcTrackLabel==digitTrackArray[0] || tpcTrackLabel==digitTrackArray[1] || tpcTrackLabel==digitTrackArray[0]) matching=3;
440         oTOFtracks->UpdateTrack(tofDigitTrackLabel, matching, timeOfFlight);
441       } else {
442         // fill the TClonesArray for tracks with no time of flight
443         Int_t matching=1;
444         oTOFtracks->SetMatchingStatus(matching);
445       }
446       
447       // delete used memory for tmp arrays used by DigitFinder
448       delete secArray;
449       delete plaArray;
450       delete strArray;
451       delete pdzArray;
452       delete pdxArray;
453       
454     } // close the if for succeded backpropagation in TOF acceptance along z
455     
456   } // end loop on reconstructed tracks
457   
458   // free used memory for digitmap
459   delete fDigitsMap;
460
461   // save array with TOF tracks
462   Int_t output=SaveTracks();
463   if(output) cout << "Error writing TOF tracks " << endl;
464 }
465
466
467 //__________________________________________________________________
468 void AliTOFReconstructionerV2::Init(Option_t* opt)
469 {
470   //
471   // Initialize the AliTOFReconstructionerV2
472   //
473   //
474 }
475
476
477 //__________________________________________________________________
478 Int_t AliTOFReconstructionerV2::LoadTPCTracks()
479 {
480   //
481   // Connect the tree and the branch
482   // with reconstructed tracks
483   // backpropagated in the TPC
484
485   gBenchmark->Start("LoadTPCTracks");
486
487   TFile *kalFile    = TFile::Open(fTPCBackTracksFile.Data());
488   if (!kalFile->IsOpen()) {cerr<<"Can't open AliTPCBackTracks.root !\n"; return 3;}
489
490   // tracks from Kalman
491   Int_t event=0;
492   char treename[100]; sprintf(treename,"TreeT_TPCb_%d",event);
493   fKalmanTree=(TTree*)kalFile->Get(treename);
494   if (!fKalmanTree) {cerr<<"Can't get a tree with TPC back tracks !\n"; return 4;}
495
496   // otherwise you get always 0 for 1/pt
497   AliKalmanTrack::SetConvConst(100/0.299792458/fField);
498
499   fBranchWithTracks=fKalmanTree->GetBranch("tracks");
500   Int_t kalEntries =(Int_t)fKalmanTree->GetEntries();
501   cout<<"Number of loaded Tracks :"<< kalEntries <<endl;
502
503   gBenchmark->Stop("LoadTPCTracks");
504   gBenchmark->Show("LoadTPCTracks");   
505   return 0;
506 }
507
508 //__________________________________________________________________
509 Int_t AliTOFReconstructionerV2::LoadTRDTracks()
510 {
511   //
512   // Connect the tree and the branch
513   // with reconstructed tracks in TRD
514
515   gBenchmark->Start("LoadTRDTracks");
516
517   Int_t nEvent = 0;
518   const Int_t nPrimaries = 84210/16;
519   const Int_t maxIndex = nPrimaries;
520   Int_t rtIndex[maxIndex];
521
522   TFile *tf=TFile::Open("AliTRDtracks.root");
523
524   if (!tf->IsOpen()) {cerr<<"Can't open AliTRDtracks.root !\n"; return 3;}
525   TObjArray tarray(2000);
526   char   tname[100];
527   sprintf(tname,"TRDb_%d",nEvent);     
528   TTree *tracktree=(TTree*)tf->Get(tname);
529
530   TBranch *tbranch=tracktree->GetBranch("tracks");
531
532   Int_t nRecTracks = (Int_t) tracktree->GetEntries();
533   cerr<<"Found "<<nRecTracks<<" entries in the track tree"<<endl;
534
535   for (Int_t i=0; i<nRecTracks; i++) {
536     AliTRDtrack *iotrack=new AliTRDtrack();
537     tbranch->SetAddress(&iotrack);
538     tracktree->GetEvent(i);
539     tarray.AddLast(iotrack);
540     Int_t trackLabel = iotrack->GetLabel();
541
542     //    printf("rt with %d clusters and label %d \n",
543     //     iotrack->GetNumberOfClusters(), trackLabel);
544
545     if(trackLabel < 0) continue;
546     if(trackLabel >= maxIndex) continue;
547     rtIndex[trackLabel] = i;
548   }
549   tf->Close();                 
550   gBenchmark->Stop("LoadTRDTracks");
551   gBenchmark->Show("LoadTRDTracks");   
552   return 0;
553
554 }
555
556 //__________________________________________________________________
557 Int_t AliTOFReconstructionerV2::LoadTOFDigits()
558 {
559   //
560   // Connect the TClonesArray with TOF
561   // digits and fill the digit map
562   // used by the DigitFinder
563
564   Int_t rc=0;
565
566   gBenchmark->Start("LoadTOFDigits");
567   
568   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(fTOFDigitsFile.Data());
569   if(file){
570     cout<<"headerFile already open \n";
571   }
572   else {
573     if(!file)file=TFile::Open(fTOFDigitsFile.Data());
574   }
575   
576   // Get AliRun object from file
577   if (!gAlice) {
578     gAlice = (AliRun*)file->Get("gAlice");
579     if (gAlice) printf("AliRun object found on file\n");
580   }
581   
582   Int_t iEvNum = 0;
583   if (iEvNum == 0) iEvNum = (Int_t) gAlice->TreeE()->GetEntries();
584
585
586   AliTOFdigit *tofdigit;
587
588   AliTOF * tof = (AliTOF *) gAlice->GetDetector("TOF") ;
589
590   if (!tof) {
591     cout << "<LoadTOFDigits> No TOF detector found" << endl;
592     rc = 2;
593     return rc;
594   }
595
596   for (Int_t ievent = 0; ievent < iEvNum; ievent++) {
597
598     gAlice->GetEvent(ievent) ;
599     if(gAlice->TreeD()==0) {
600       cout << "<LoadTOFDigits> No  TreeD found" << endl;
601       rc = 4;
602       return rc;
603     }
604     
605
606     Int_t ndig;
607     gAlice->ResetDigits();
608     gAlice->TreeD()->GetEvent(ievent);
609     fTOFDigits   = tof->Digits();
610     
611     ndig=fTOFDigits->GetEntries();
612
613     // create the digit map
614     fDigitsMap = new AliTOFDigitMap(fTOFDigits);
615
616     
617     cout << "<LoadTOFDigits> found " << ndig
618          << " TOF digits for event " << ievent << endl;
619     
620     for (Int_t k=0; k<ndig; k++) {
621       tofdigit= (AliTOFdigit*) fTOFDigits->UncheckedAt(k);
622       Float_t tdc=tofdigit->GetTdc();
623       // adc value can be used for weighting
624       //Float_t adc=tofdigit->GetAdc();
625
626       // TOF digit volumes
627       Int_t    vol[5];       // location for a digit
628       Int_t sector    = tofdigit->GetSector(); // range [1-18]
629       Int_t plate     = tofdigit->GetPlate();  // range [1- 5]
630       Int_t strip     = tofdigit->GetStrip();  // range [1-20]
631       Int_t padx      = tofdigit->GetPadx();   // range [1-48]
632       Int_t padz      = tofdigit->GetPadz();   // range [1- 2]
633
634       vol[0] = sector;
635       vol[1] = plate;
636       vol[2] = strip;
637       vol[3] = padx;
638       vol[4] = padz;
639
640       // QA
641       Bool_t isDigitBad = (sector<1 || sector>18 || plate<1 || plate >5 || padz<1 || padz>2 || padx<1 || padx>48);
642
643       if (isDigitBad) {
644         cout << "<LoadTOFDigits>  strange digit found" << endl;
645         rc = 3;
646         return rc;
647       } // if (isDigitBad)
648
649       // Fill the digit map checking if the location is already used
650       // in this case we take the earliest signal
651       if (fDigitsMap->TestHit(vol) != kEmpty) {
652         // start comparison in between the 2 digit
653         AliTOFdigit *dig = static_cast<AliTOFdigit*>(fDigitsMap->GetHit(vol));
654         if(tdc < (dig->GetTdc())) fDigitsMap->SetHit(vol,k);
655         // we can add also the check on adc value
656         // by selecting the largest adc value
657       } else {
658         fDigitsMap->SetHit(vol,k);
659       }
660       // to be added protection versus 2-digit on the same pad
661       // we have to have memory also of the second digit
662       
663     } // for (k=0; k<ndig; k++)
664   
665   } // end loop on events
666
667   gBenchmark->Stop("LoadTOFDigits");
668   gBenchmark->Show("LoadTOFDigits");   
669   return rc;
670 }
671
672 //__________________________________________________________________
673 void AliTOFReconstructionerV2::IsInsideThePad(Float_t x, Float_t y, Float_t z, Int_t *nGeom, Float_t& zPad, Float_t& xPad) 
674 {
675   //   input: x,y,z - coordinates of a point in the mrf [cm]
676   //   output: array  nGeom[]
677   //          nGeom[0] - the TOF sector number, 1,2,...,18 along azimuthal direction starting from -90 deg.
678   //          nGeom[1] - the TOF module number, 1,2,3,4,5=C,B,A,B,C along z-direction
679   //          nGeom[2] - the TOF strip  number, 1,2,... along z-direction
680   //          nGeom[3] - the TOF padz  number,  1,2=NPZ across a strip
681   //          nGeom[4] - the TOF padx  number,  1,2,...,48=NPX along a strip
682   //          zPad, xPad - coordinates of the hit in the pad frame
683   //  numbering is adopted for the version 3.08 of AliRoot
684   //  example:
685   //   from Hits: sec,pla,str,padz,padx=4,2,14,2,35
686   //  Vol. n.0: ALIC, copy number 1
687   //  Vol. n.1: B077, copy number 1
688   //  Vol. n.2: B074, copy number 5
689   //  Vol. n.3: BTO2, copy number 1
690   //  Vol. n.4: FTOB, copy number 2
691   //  Vol. n.5: FLTB, copy number 0
692   //  Vol. n.6: FSTR, copy number 14
693   //  Vol. n.7: FSEN, copy number 0
694   //  Vol. n.8: FSEZ, copy number 2
695   //  Vol. n.9: FSEX, copy number 35
696   //  Vol. n.10: FPAD, copy number 0
697
698
699   Float_t xTOF[3];
700   Int_t sector=0,module=0,strip=0,padz=0,padx=0;
701   Int_t i,numed,nLevel,copyNumber;
702   Gcvolu_t* gcvolu;
703   char name[5];
704   name[4]=0;
705   
706   for (i=0; i<AliTOFConstants::fgkmaxtoftree; i++) nGeom[i]=0;
707   zPad=100.;
708   xPad=100.;
709   
710   xTOF[0]=x;
711   xTOF[1]=y;
712   xTOF[2]=z;
713   
714   TGeant3 * fG3Geom = (TGeant3*) gMC;
715
716   fG3Geom->Gmedia(xTOF, numed);
717   gcvolu=fG3Geom->Gcvolu();
718   nLevel=gcvolu->nlevel;
719   if(fdbg) {
720     for (Int_t i=0; i<nLevel; i++) {
721       strncpy(name,(char*) (&gcvolu->names[i]),4);
722       cout<<"Vol. n."<<i<<": "<<name<<", copy number "<<gcvolu->number[i]<<endl;
723     }
724   }
725   if(nLevel>=2) {
726     // sector type name: B071(1,2,...,10),B074(1,2,3,4,5-PHOS),B075(1,2,3-RICH)
727     strncpy(name,(char*) (&gcvolu->names[2]),4);
728     // volume copy: 1,2,...,10 for B071, 1,2,3,4,5 for B074, 1,2,3 for B075
729     copyNumber=gcvolu->number[2];
730    if(!strcmp(name,"B071")) {
731      if (copyNumber>=6 && copyNumber<=8) {
732        sector=copyNumber+10;
733      } else if (copyNumber>=1 && copyNumber<=5){
734        sector=copyNumber+7;
735      } else {
736        sector=copyNumber-8;
737      }
738    } else if(!strcmp(name,"B075")) {
739      sector=copyNumber+12;
740    } else if(!strcmp(name,"B074")) {
741      if (copyNumber>=1 && copyNumber<=3){
742        sector=copyNumber+4;
743      } else {
744        sector=copyNumber-1;
745      }
746    }
747   }
748   if(sector) {
749     nGeom[0]=sector;
750     if(nLevel>=4) {
751       // we'll use the module value in z-direction:
752       //                                    1    2    3    4    5
753       // the module order in z-direction: FTOC,FTOB,FTOA,FTOB,FTOC
754       // the module copy:                   2    2    0    1    1
755       // module type name: FTOA, FTOB, FTOC
756       strncpy(name,(char*) (&gcvolu->names[4]),4);
757       // module copy:  
758       copyNumber=gcvolu->number[4];
759       if(!strcmp(name,"FTOC")) {
760         if (copyNumber==2) {
761           module=1;
762         } else {
763           module=5;
764         }
765       } else if(!strcmp(name,"FTOB")) {
766         if (copyNumber==2) {
767           module=2;
768         } else {
769           module=4;
770         }
771       } else if(!strcmp(name,"FTOA")) {
772         module=3;
773       }
774     }
775   }
776   
777   if(module) {
778     nGeom[1]=module;
779     if(nLevel>=6) {
780       // strip type name: FSTR
781       strncpy(name,(char*) (&gcvolu->names[6]),4);
782       // strip copy:  
783       copyNumber=gcvolu->number[6];
784       if(!strcmp(name,"FSTR")) strip=copyNumber; 
785     }
786   }
787   
788   if(strip) {
789     nGeom[2]=strip;
790     if(nLevel>=8) {
791       // padz type name: FSEZ
792       strncpy(name,(char*) (&gcvolu->names[8]),4);
793       // padz copy:  
794       copyNumber=gcvolu->number[8];
795       if(!strcmp(name,"FSEZ")) padz=copyNumber; 
796     }
797   }
798   if(padz) {
799     nGeom[3]=padz;
800     if(nLevel>=9) {
801       // padx type name: FSEX
802       strncpy(name,(char*) (&gcvolu->names[9]),4);
803       // padx copy:  
804       copyNumber=gcvolu->number[9];
805       if(!strcmp(name,"FSEX")) padx=copyNumber; 
806     }
807   }
808   
809   if(padx) {
810     nGeom[4]=padx;
811     zPad=gcvolu->glx[2];
812     xPad=gcvolu->glx[0];
813   }
814   
815 }
816
817 //__________________________________________________________________
818 void AliTOFReconstructionerV2::GetGlobalXYZ(Double_t alpha, Double_t& x, Double_t& y, Double_t& z)
819 {
820   //
821   // return the current running coordinates of 
822   // the track in the global reference frame
823   // x, y and z have to initialized to the
824   // local frame coordinates by the caller
825   // alpha is the alpha coordinate in the TPC Kalman
826   // reference frame 
827
828   // it take into account differences in between
829   // TPC and TRD local coordinates frames
830   if      (alpha < -TMath::Pi()) alpha += 2*TMath::Pi();
831   else if (alpha >= TMath::Pi()) alpha -= 2*TMath::Pi();
832
833   // conversion
834   Double_t tmp=x*TMath::Cos(alpha) - y*TMath::Sin(alpha);
835   y=x*TMath::Sin(alpha) + y*TMath::Cos(alpha);
836   x=tmp;            
837 }
838
839 //__________________________________________________________________
840 Bool_t AliTOFReconstructionerV2::DigitFinder(TArrayI *secArray, TArrayI *plaArray, TArrayI *strArray, TArrayI *pdzArray, TArrayI *pdxArray, Int_t* assignedVol, Int_t* digitTrackArray, Float_t& tdc)
841 {
842   //
843   // Description
844   // input: arrays with sectors, plates, strips, padz, padx
845   // found during fine-propagation of probe tracks
846   //
847   // output kFALSE if signal is not found
848   //        kTRUE  if signal is found
849   // in this case the assignedVol array contains the digit volume numbers
850   // and digitTrackArray the track numbers (max 3) contributing to the
851   // digit
852
853   Bool_t dummy=kFALSE;
854   Int_t nFilledSlot=secArray->GetSize();
855
856   // start loop
857   Float_t maxWeight=-1.;
858   Int_t   indexOfMaxWeight=-1;
859   for (Int_t i = 0; i < nFilledSlot; i++) {
860     Int_t    vol[5];       // location for a digit
861     vol[0] = (*secArray)[i];
862     vol[1] = (*plaArray)[i];
863     vol[2] = (*strArray)[i];
864     vol[3] = (*pdxArray)[i];
865     vol[4] = (*pdzArray)[i];
866
867     // check for digit in the current location
868     if (fDigitsMap->TestHit(vol) != kEmpty) {
869
870       AliTOFdigit *dig = static_cast<AliTOFdigit*>(fDigitsMap->GetHit(vol));
871       Float_t adcWeight=dig->GetAdc();
872       if(adcWeight > maxWeight){
873         maxWeight=adcWeight;
874         indexOfMaxWeight=i;
875         tdc=dig->GetTdc();
876         digitTrackArray=dig->GetTracks();
877  
878       } // if(adcWeight > maxWeight)
879     } // close if (fDigitsMap->TestHit(vol) != kEmpty)
880
881   } // end loop
882
883   if(indexOfMaxWeight!=-1){
884     assignedVol[0]=(*secArray)[indexOfMaxWeight];
885     assignedVol[1]=(*plaArray)[indexOfMaxWeight];
886     assignedVol[2]=(*strArray)[indexOfMaxWeight];
887     assignedVol[3]=(*pdxArray)[indexOfMaxWeight];
888     assignedVol[4]=(*pdzArray)[indexOfMaxWeight];
889     dummy=kTRUE;
890   }
891
892   return dummy;
893 }
894
895 //__________________________________________________________________
896 Int_t AliTOFReconstructionerV2::SaveTracks(const Char_t *outname, const Int_t split)
897 {
898   //
899   // save reconstructed tracks into 
900   // outname file
901   //
902   TDirectory *savedir=gDirectory;
903   const Char_t *name="Writing Output";
904   cerr<<'\n'<<name<<"...\n";
905   gBenchmark->Start(name);
906   
907   TFile *out=TFile::Open(outname,"RECREATE");
908   if (!out->IsOpen()) {
909     cerr<<"AliTOFReconstructionerV2::SaveTracks(): ";
910     cerr<<"file for TOF tracks is not open !\n";
911      return 2;
912   }
913   
914   out->cd();
915   TTree T("T","tree with TOF tracks");
916   T.Branch("tracks",&fTOFTracks,256000,split);
917
918
919   T.Fill();
920   T.Write();
921   savedir->cd();
922   out->Close();
923   return 0;
924   gBenchmark->Stop(name);
925   gBenchmark->Show(name);
926 }
927
928 //__________________________________________________________________
929 void AliTOFReconstructionerV2::Comparison(Int_t* rtIndex)
930 {
931   //
932   // perform MC comparison
933   // used also for track length
934   // for the time being
935   // Connect the AliRoot file containing Geometry, Kine, Hits, and Digits
936
937
938   Int_t nEvent = 0;
939   Char_t *alifile = "galice.root";
940   
941   TFile *gafl = (TFile*) gROOT->GetListOfFiles()->FindObject(alifile);
942   if (!gafl) {
943     cout << "Open the ALIROOT-file " << alifile << endl;
944     gafl = new TFile(alifile);
945   }
946   else {
947     cout << alifile << " is already open" << endl;
948   }
949   
950   // Get AliRun object from file or create it if not on file
951   gAlice = (AliRun*) gafl->Get("gAlice");
952   if (gAlice)
953     cout << "AliRun object found on file" << endl;
954   else
955     gAlice = new AliRun("gAlice","Alice test program");
956   
957   
958   // Import the Trees for the event nEvent in the file
959   const Int_t nparticles = gAlice->GetEvent(nEvent);
960   if (nparticles <= 0) return;
961
962   // Get pointers to Alice detectors and Hits containers
963   AliDetector *TOF  = gAlice->GetDetector("TOF");
964   Int_t ntracks    = (Int_t) gAlice->TreeH()->GetEntries();     
965   
966   // Start loop on tracks in the hits containers
967   for (Int_t track=0; track < ntracks; track++) {
968     
969     if(TOF) {
970       for(AliTOFhitT0* tofHit = (AliTOFhitT0*)TOF->FirstHit(track); 
971           tofHit; 
972           tofHit=(AliTOFhitT0*)TOF->NextHit()) {
973         
974          Int_t ipart    = tofHit->GetTrack();
975          if(ipart >= 80000) continue;
976          if(rtIndex[ipart] < 0) continue; 
977
978          TParticle *part = gAlice->Particle(ipart);
979          
980          // get first the pdg code
981          Int_t pdgCode=part->GetPdgCode();
982          
983          // then track length
984          Float_t trackLength=tofHit->GetLen(); // [cm]
985
986          // update the tof TClonesArray with TOF tracks
987          AliTOFTrackV2* oTOFtracks=(AliTOFTrackV2*)fTOFTracks->UncheckedAt(rtIndex[ipart]);
988          oTOFtracks->UpdateTrack(pdgCode,trackLength);
989
990       } // loop on hits connected to the current track
991     } // if(TOF)
992   } // end loop on primary tracks
993 }
994
995 //__________________________________________________________________
996 Bool_t AliTOFReconstructionerV2::operator==(const AliTOFReconstructionerV2 & tofrecv2)const
997 {
998   // Equal operator.
999   // Reconstructioner are equal if their fField, fNDummyTracks, fScaleSigmaFactor and fStep are equal
1000  
1001   if( (fField==tofrecv2.fField)&&(fNDummyTracks==tofrecv2.fNDummyTracks)&&(fScaleSigmaFactor==tofrecv2.fScaleSigmaFactor)&&(fStep==tofrecv2.fStep))
1002     return kTRUE ;
1003   else
1004     return kFALSE ;
1005 }