TOF/tzero15tracks.C

   1 void tzero15tracks(const char* datafile, Float_t globalTimeRes=1.2e-10)
   2 {
   3
   4   //
   5   // Calculate the t0 distribution
   6   // selecting 15 tracks
   7   // (primary particles reaching TOF)
   8   // NB: see the description of the analogous AliTOFT0 class
   9   // NB II: apply this macro to a realistic Pb-Pb event (e.g. Hijing event)
  10   // Use case
  11   // - start root
  12   // root [0] .L tzero15tracksopt.C
  13   // // exec the macro with the default time resolution 120 ps
  14   // root [1] tzero15tracksopt("hij25evCentralTOFv4T0Set10.04T.root")
  15
  16   // // exec the macro with 150 ps global time resolution
  17   // // NB: time resolution has to be given in [s]
  18   // root [1] tzero15tracksopt("hij25evCentralTOFv4T0Set10.04T.root", 1.5e-10))
  19
  20
  21   // Dynamically link some shared libs
  22   if (gClassTable->GetID("AliRun") < 0) {
  23     gROOT->LoadMacro("loadlibs.C");
  24     loadlibs();
  25   }
  26
  27   // Connect the Root Galice file containing Geometry, Kine and Hits
  28   Int_t ngood;
  29   Int_t nmisidentified=0;
  30   Int_t nmisidentified0=0;
  31   Int_t nmisidentified1=0;
  32   Int_t nmisidentified2=0;
  33   Int_t nmisidentified3=0;
  34   Int_t nmisidentified4=0;
  35   Int_t nmisidentified5=0;
  36   Int_t nmisidentified6=0;
  37   Int_t nmisidentified7=0;
  38   Int_t nmisidentified8=0;
  39   Int_t nmisidentified9=0;
  40   Int_t nmisidentified10=0;
  41   Int_t nmisidentified11=0;
  42   Int_t nmisidentified12=0;
  43   Int_t nmisidentified13=0;
  44   Int_t nmisidentified14=0;
  45   Int_t nbytes = 0;
  46   Int_t ipartold = -1;
  47   Int_t dblehit = 0;
  48   Int_t j,hit,ipart;
  49   Int_t nhits;
  50   Int_t selected=0;
  51   Int_t istop=0;
  52   const Float_t timeresolution=globalTimeRes; // in [s]
  53   Float_t timeresolutioninns=timeresolution*(1.e+9); // convert in [ns]
  54   cout << "Global Time Resolution " << timeresolutioninns << " ns" << endl; //
  55
  56   // output (histos!) filename
  57   char outFileName[100];
  58   strcpy(outFileName,"ht015tr");
  59   strcat(outFileName,datafile);
  60
  61   gBenchmark->Start("TOFT0");
  62
  63   TH1F* htzerobest= new TH1F("htzerobest","T0 for best assignment",200,-1.,1.);
  64   TH1F* hchibest  = new TH1F("hchibest","ChiSquare Min Distribution",80,0.,40.);
  65   TH1F* hchibestconflevel  = new TH1F("hchibestconflevel","ChiSquare Min Confidence Level",10,0.,1.);
  66
  67   Float_t t0best=999.;
  68   Int_t assparticle[15]={3,3,3,3,3,3,3,3,3,3,3,3,3,3,3};
  69   Int_t truparticle[15]={3,3,3,3,3,3,3,3,3,3,3,3,3,3,3};
  70   Float_t timeofflight[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  71   Float_t momentum[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  72   Float_t timezero[15];
  73   Float_t weightedtimezero[15];
  74   Float_t beta[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  75   Float_t sqMomError[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  76   Float_t sqTrackError[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  77   Float_t massarray[3]={0.13957,0.493677,0.9382723};
  78   Float_t chisquare=999.;
  79   Float_t tracktoflen[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
  80   TRandom* rnd = new TRandom();
  81   TParticle *particle;
  82   AliTOFhitT0  *tofHit;
  83
  84   TFile *file = TFile::Open(datafile,"old");
  85
  86   // Get AliRun object from file or create it if not on file
  87   gAlice = (AliRun*)file->Get("gAlice");
  88
  89   // Import the Kine and Hits Trees for the event evNumber in the file
  90
  91   Int_t nallevents = (Int_t) gAlice->TreeE()->GetEntries();
  92
  93   for (Int_t nev=0; nev<nallevents;nev++) {
  94
  95     Int_t nparticles = gAlice->GetEvent(nev);
  96     if (nparticles <= 0) return;
  97     //     ngood=0;
  98     // Get pointers to Alice detectors and Hits containers
  99
 100     AliDetector *TOF  = gAlice->GetDetector("TOF");
 101     TObjArray *Particles = gAlice->Particles();
 102
 103     if (TOF) TClonesArray *TOFhits   = TOF->Hits();
 104
 105     TTree *TH = gAlice->TreeH();
 106     Int_t ntracks = TH->GetEntries();
 107
 108
 109     // Start loop on primary tracks in the hits containers
 110     Int_t lasttrack=-1;
 111     Int_t nset=0;
 112     for (Int_t track=lasttrack+1; track<ntracks;track++) {
 113
 114       if(nset>=3) break; // check on the number of set analyzed
 115
 116       gAlice->ResetHits();
 117       nbytes += TH->GetEvent(track);
 118
 119       // =======>Histogram TOF
 120       if (TOF) {
 121         nhits = TOFhits->GetEntriesFast();
 122
 123         for (hit=0;hit<nhits;hit++) {
 124
 125           tofHit   = (AliTOFhitT0*)TOFhits->UncheckedAt(hit);
 126           ipart    = tofHit->GetTrack();
 127
 128           // check to discard the case when the same particle is selected more than one
 129           // time
 130
 131           if (ipart != ipartold){
 132
 133             particle = (TParticle*)gAlice->Particle(ipart);
 134
 135             Float_t idealtime=tofHit->GetTof();
 136             //     Float_t time=idealtime;
 137             Float_t time   = rnd->Gaus(idealtime, timeresolution);
 138             Float_t toflen=tofHit->GetLen();
 139             toflen=toflen/100.; // toflen given in m because GEANT gives cm!
 140             Int_t pdg   = particle->GetPdgCode();
 141             Int_t abspdg   =TMath::Abs(pdg);
 142             Float_t idealmom  = particle->P();
 143             Float_t momres=idealmom*0.025; // 2.5% res token into account for all momenta
 144             //Float_t mom =idealmom+2.*(rnd->Rndm()-0.5)*momres;
 145             Float_t mom =rnd->Gaus(idealmom,momres);
 146             Bool_t isgoodpart=(abspdg==211 || abspdg==2212 || abspdg==321);
 147
 148             time*=1.e+9; // tof given in nanoseconds
 149             if (particle->GetFirstMother() < 0 && isgoodpart && mom<=1.75 && mom>=1.25){
 150
 151               selected+=1;
 152               istop=selected;
 153               if(istop>15) break;
 154               Int_t index=selected-1;
 155               timeofflight[index]=time;
 156               tracktoflen[index]=toflen;
 157               momentum[index]=mom;
 158               switch (abspdg) {
 159               case 211:
 160                 truparticle[index]=0;
 161                 break ;
 162               case 321:
 163                 truparticle[index]=1;
 164                 break ;
 165               case 2212:
 166                 truparticle[index]=2;
 167                 break ;
 168               }
 169
 170             }
 171             ipartold = ipart;
 172
 173             if(istop==15){ // start analysis on current set
 174               nset+=1;
 175               lasttrack=track;
 176               istop=0;
 177               selected=0;
 178               cout << "starting t0 calculation for current set" << endl;
 179               for (Int_t i1=0; i1<3;i1++) {
 180                 beta[0]=momentum[0]/TMath::Sqrt(massarray[i1]*massarray[i1]+momentum[0]*momentum[0]);
 181                 for (Int_t i2=0; i2<3;i2++) {
 182                   beta[1]=momentum[1]/TMath::Sqrt(massarray[i2]*massarray[i2]+momentum[1]*momentum[1]);
 183                   for (Int_t i3=0; i3<3;i3++) {
 184                     beta[2]=momentum[2]/TMath::Sqrt(massarray[i3]*massarray[i3]+momentum[2]*momentum[2]);
 185                     for (Int_t i4=0; i4<3;i4++) {
 186                       beta[3]=momentum[3]/TMath::Sqrt(massarray[i4]*massarray[i4]+momentum[3]*momentum[3]);
 187                       for (Int_t i5=0; i5<3;i5++) {
 188                         beta[4]=momentum[4]/TMath::Sqrt(massarray[i5]*massarray[i5]+momentum[4]*momentum[4]);
 189                         for (Int_t i6=0; i6<3;i6++) {
 190                           beta[5]=momentum[5]/TMath::Sqrt(massarray[i6]*massarray[i6]+momentum[5]*momentum[5]);
 191                           for (Int_t i7=0; i7<3;i7++) {
 192                             beta[6]=momentum[6]/TMath::Sqrt(massarray[i7]*massarray[i7]+momentum[6]*momentum[6]);
 193                             for (Int_t i8=0; i8<3;i8++) {
 194                               beta[7]=momentum[7]/TMath::Sqrt(massarray[i8]*massarray[i8]+momentum[7]*momentum[7]);
 195                               for (Int_t i9=0; i9<3;i9++) {
 196                                 beta[8]=momentum[8]/TMath::Sqrt(massarray[i9]*massarray[i9]+momentum[8]*momentum[8]);
 197                                 for (Int_t i10=0; i10<3;i10++) {
 198                                   beta[9]=momentum[9]/TMath::Sqrt(massarray[i10]*massarray[i10]+momentum[9]*momentum[9]);
 199                                   for (Int_t i11=0; i11<3;i11++) {
 200                                     beta[10]=momentum[10]/TMath::Sqrt(massarray[i11]*massarray[i11]+momentum[10]*momentum[10]);
 201                                     for (Int_t i12=0; i12<3;i12++) {
 202                                       beta[11]=momentum[11]/TMath::Sqrt(massarray[i12]*massarray[i12]+momentum[11]*momentum[11]);
 203                                       for (Int_t i13=0; i13<3;i13++) {
 204                                         beta[12]=momentum[12]/TMath::Sqrt(massarray[i13]*massarray[i13]+momentum[12]*momentum[12]);
 205                                         for (Int_t i14=0; i14<3;i14++) {
 206                                           beta[13]=momentum[13]/TMath::Sqrt(massarray[i14]*massarray[i14]+momentum[13]*momentum[13]);
 207                                           for (Int_t i15=0; i15<3;i15++) {
 208                                             beta[14]=momentum[14]/TMath::Sqrt(massarray[i15]*massarray[i15]+momentum[14]*momentum[14]);
 209
 210                                             Float_t meantzero=0.;
 211                                             Float_t sumAllweights=0.;
 212                                             for (Int_t itz=0; itz<15;itz++) {
 213                                               sqMomError[itz]=((1.-beta[itz]*beta[itz])*0.025)*((1.-beta[itz]*beta[itz])*0.025)*(tracktoflen[itz]/(0.299792*beta[itz]))*(tracktoflen[itz]/(0.299792*beta[itz])); // this gives the square of the momentum error in nanoseconds
 214                                               sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); // total error for the current track
 215                                               sumAllweights+=1./sqTrackError[itz];
 216                                               timezero[itz]=(tracktoflen[itz]/(beta[itz]*0.299792))-timeofflight[itz];
 217                                               weightedtimezero[itz]=((tracktoflen[itz]/(beta[itz]*0.299792))-timeofflight[itz])/sqTrackError[itz];// weighted time zero for current track
 218                                               meantzero+=weightedtimezero[itz];
 219                                             }
 220
 221                                             meantzero=meantzero/sumAllweights; // it is given in [ns]
 222
 223                                             Float_t dummychisquare=0.;
 224                                             // calculate the chisquare for the current assignment
 225                                             for (Int_t itz=0; itz<15;itz++) {
 226                                               dummychisquare+=(timezero[itz]-meantzero)*(timezero[itz]-meantzero)/sqTrackError[itz];
 227                                             }
 228
 229
 230                                             if(dummychisquare<=chisquare){
 231                                               assparticle[0]=i1;
 232                                               assparticle[1]=i2;
 233                                               assparticle[2]=i3;
 234                                               assparticle[3]=i4;
 235                                               assparticle[4]=i5;
 236                                               assparticle[5]=i6;
 237                                               assparticle[6]=i7;
 238                                               assparticle[7]=i8;
 239                                               assparticle[8]=i9;
 240                                               assparticle[9]=i10;
 241                                               assparticle[10]=i11;
 242                                               assparticle[11]=i12;
 243                                               assparticle[12]=i13;
 244                                               assparticle[13]=i14;
 245                                               assparticle[14]=i15;
 246                                               chisquare=dummychisquare;
 247                                               t0best=meantzero;
 248                                             }
 249                                           }
 250                                         }
 251                                       }
 252                                     }
 253                                   }
 254                                 }
 255                               }
 256                             }
 257                           }
 258                         }
 259
 260                       }
 261                     }
 262                   }
 263                 }
 264               }
 265
 266               cout << "true" << truparticle[0] << truparticle[1] << truparticle[2] << truparticle[3] << truparticle[4] << truparticle[5] << truparticle[6] << truparticle[7] << truparticle[8] << truparticle[9] <<truparticle[10] << truparticle[11] << truparticle[1
 267 2] << truparticle[13] << truparticle[14] <<endl;
 268               cout << "best" << assparticle[0] << assparticle[1] << assparticle[2] << assparticle[3] << assparticle[4] << assparticle[5] << assparticle[6] << assparticle[7] << assparticle[8] << assparticle[9] <<assparticle[10] << assparticle[11] << assparticle[1
 269 2] << assparticle[13] << assparticle[14] << endl;
 270               if(truparticle[0]==assparticle[0] && truparticle[1]==assparticle[1] && truparticle[2]==assparticle[2]  && truparticle[3]==assparticle[3] && truparticle[4]==assparticle[4]&& truparticle[5]==assparticle[5] && truparticle[6]==assparticle[6] && trupart
 271 icle[7]==assparticle[7]  && truparticle[8]==assparticle[8] && truparticle[9]==assparticle[9] &&truparticle[10]==assparticle[10] && truparticle[11]==assparticle[11] && truparticle[12]==assparticle[12]  && truparticle[13]==assparticle[13] && truparticle[14]
 272 ==assparticle[14]) ngood+=1;
 273               if(truparticle[0]!=assparticle[0]) nmisidentified0+=1;
 274               if(truparticle[1]!=assparticle[1]) nmisidentified1+=1;
 275               if(truparticle[2]!=assparticle[2]) nmisidentified2+=1;
 276               if(truparticle[3]!=assparticle[3]) nmisidentified3+=1;
 277               if(truparticle[4]!=assparticle[4]) nmisidentified4+=1;
 278               if(truparticle[5]!=assparticle[5]) nmisidentified5+=1;
 279               if(truparticle[6]!=assparticle[6]) nmisidentified6+=1;
 280               if(truparticle[7]!=assparticle[7]) nmisidentified7+=1;
 281               if(truparticle[8]!=assparticle[8]) nmisidentified8+=1;
 282               if(truparticle[9]!=assparticle[9]) nmisidentified9+=1;
 283               if(truparticle[10]!=assparticle[10]) nmisidentified10+=1;
 284               if(truparticle[11]!=assparticle[11]) nmisidentified11+=1;
 285               if(truparticle[12]!=assparticle[12]) nmisidentified12+=1;
 286               if(truparticle[13]!=assparticle[13]) nmisidentified13+=1;
 287               if(truparticle[14]!=assparticle[14]) nmisidentified14+=1;
 288               cout << "chisquare for current set" << chisquare << endl;
 289               htzerobest->Fill(t0best);
 290               hchibest->Fill(chisquare);
 291               Double_t dblechisquare=(Double_t)chisquare;
 292               Float_t confLevel=(Float_t)TMath::Prob(dblechisquare,14);
 293               cout << "conf level " << confLevel << endl;
 294               hchibestconflevel->Fill(confLevel);
 295               chisquare=999.;
 296               t0best=999.;
 297
 298             } // end for the current set. close if(istop==15)
 299
 300
 301           } // end condition on ipartold
 302
 303
 304         }// end loop on hits for the current track
 305       }// end TOF
 306       if(istop>=15) break;
 307     }// end loop on Tracks
 308     cout << "ngood for current event " << ngood << endl;
 309
 310
 311   }//end loop on events
 312
 313   nmisidentified=(nmisidentified0+nmisidentified1+nmisidentified2+nmisidentified3+nmisidentified4+nmisidentified5+nmisidentified6+nmisidentified7+nmisidentified8+nmisidentified9+nmisidentified10+nmisidentified11+nmisidentified12+nmisidentified13+nmisident
 314 ified14);
 315   cout << "total misidentified " << nmisidentified << endl;
 316   TFile *houtfile = new TFile(outFileName,"recreate");
 317   houtfile->cd();
 318
 319   htzerobest->Write();
 320   hchibest->Write();
 321   hchibestconflevel->Write();
 322   houtfile->Close();
 323
 324   gBenchmark->Stop("TOFT0");
 325   cout << "15 tracks:" << endl ;
 326   cout << "   took " << gBenchmark->GetCpuTime("TOFT0") << " seconds in order to calculate T0 "
 327        <<  gBenchmark->GetCpuTime("TOFT0")/nallevents << " seconds per event " << endl ;
 328   cout << endl ;
 329
 330 }
 331
 332