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