TOF/AliTOFReconstructionerV2.cxx

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