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