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