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