From bcb56cb9e81cdb2fb3883fa6854f5f0ab86b444b Mon Sep 17 00:00:00 2001 From: vicinanz Date: Tue, 12 Feb 2002 14:00:29 +0000 Subject: [PATCH] New t_zero macro added to repository --- TOF/tzero15tracks.C | 321 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 321 insertions(+) create mode 100644 TOF/tzero15tracks.C diff --git a/TOF/tzero15tracks.C b/TOF/tzero15tracks.C new file mode 100644 index 00000000000..a11668b4fb2 --- /dev/null +++ b/TOF/tzero15tracks.C @@ -0,0 +1,321 @@ +void tzero15tracks(const char* datafile, Float_t globalTimeRes=1.2e-10) +{ + + // + // Calculate the t0 distribution + // + // (primary particles reaching TOF) + // + + + // Dynamically link some shared libs + if (gClassTable->GetID("AliRun") < 0) { + gROOT->LoadMacro("loadlibs.C"); + loadlibs(); + } + + // Connect the Root Galice file containing Geometry, Kine and Hits + Int_t ngood; + Int_t nmisidentified=0; + Int_t nmisidentified0=0; + Int_t nmisidentified1=0; + Int_t nmisidentified2=0; + Int_t nmisidentified3=0; + Int_t nmisidentified4=0; + Int_t nmisidentified5=0; + Int_t nmisidentified6=0; + Int_t nmisidentified7=0; + Int_t nmisidentified8=0; + Int_t nmisidentified9=0; + Int_t nmisidentified10=0; + Int_t nmisidentified11=0; + Int_t nmisidentified12=0; + Int_t nmisidentified13=0; + Int_t nmisidentified14=0; + Int_t nbytes = 0; + Int_t ipartold = -1; + Int_t dblehit = 0; + Int_t j,hit,ipart; + Int_t nhits; + Int_t selected=0; + Int_t istop=0; + const Float_t timeresolution=globalTimeRes; // in [s] + Float_t timeresolutioninns=timeresolution*(1.e+9); // convert in [ns] + cout << "Global Time Resolution " << timeresolutioninns << " ns" << endl; // + + // output (histos!) filename + char outFileName[100]; + strcpy(outFileName,"ht015tr"); + strcat(outFileName,datafile); + + gBenchmark->Start("TOFT0"); + + TH1F* htzerobest= new TH1F("htzerobest","T0 for best assignment",200,-1.,1.); + TH1F* hchibest = new TH1F("hchibest","ChiSquare Min Distribution",80,0.,40.); + TH1F* hchibestconflevel = new TH1F("hchibestconflevel","ChiSquare Min Confidence Level",10,0.,1.); + + Float_t t0best=999.; + Int_t assparticle[15]={3,3,3,3,3,3,3,3,3,3,3,3,3,3,3}; + Int_t truparticle[15]={3,3,3,3,3,3,3,3,3,3,3,3,3,3,3}; + Float_t timeofflight[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; + Float_t momentum[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; + Float_t timezero[15]; + Float_t weightedtimezero[15]; + Float_t beta[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; + Float_t sqMomError[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; + Float_t sqTrackError[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; + Float_t massarray[3]={0.13957,0.493677,0.9382723}; + Float_t chisquare=999.; + Float_t tracktoflen[15]={0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}; + TRandom* rnd = new TRandom(); + TParticle *particle; + AliTOFhitT0 *tofHit; + + TFile *file = TFile::Open(datafile,"old"); + + // Get AliRun object from file or create it if not on file + gAlice = (AliRun*)file->Get("gAlice"); + + // Import the Kine and Hits Trees for the event evNumber in the file + + Int_t nallevents = (Int_t) gAlice->TreeE()->GetEntries(); + + for (Int_t nev=0; nevGetEvent(nev); + if (nparticles <= 0) return; + // ngood=0; + // Get pointers to Alice detectors and Hits containers + + AliDetector *TOF = gAlice->GetDetector("TOF"); + TObjArray *Particles = gAlice->Particles(); + + if (TOF) TClonesArray *TOFhits = TOF->Hits(); + + TTree *TH = gAlice->TreeH(); + Int_t ntracks = TH->GetEntries(); + + + // Start loop on primary tracks in the hits containers + Int_t lasttrack=-1; + Int_t nset=0; + for (Int_t track=lasttrack+1; track=3) break; // check on the number of set analyzed + + gAlice->ResetHits(); + nbytes += TH->GetEvent(track); + + // =======>Histogram TOF + if (TOF) { + nhits = TOFhits->GetEntriesFast(); + + for (hit=0;hitUncheckedAt(hit); + ipart = tofHit->GetTrack(); + + // check to discard the case when the same particle is selected more than one + // time + + if (ipart != ipartold){ + + particle = (TParticle*)gAlice->Particle(ipart); + + Float_t idealtime=tofHit->GetTof(); + // Float_t time=idealtime; + Float_t time = rnd->Gaus(idealtime, timeresolution); + Float_t toflen=tofHit->GetLen(); + toflen=toflen/100.; // toflen given in m because GEANT gives cm! + Int_t pdg = particle->GetPdgCode(); + Int_t abspdg =TMath::Abs(pdg); + Float_t idealmom = particle->P(); + Float_t momres=idealmom*0.025; // 2.5% res token into account for all momenta + //Float_t mom =idealmom+2.*(rnd->Rndm()-0.5)*momres; + Float_t mom =rnd->Gaus(idealmom,momres); + Bool_t isgoodpart=(abspdg==211 || abspdg==2212 || abspdg==321); + + time*=1.e+9; // tof given in nanoseconds + if (particle->GetFirstMother() < 0 && isgoodpart && mom<=1.7 && mom>=1.25){ + + selected+=1; + istop=selected; + if(istop>15) break; + Int_t index=selected-1; + timeofflight[index]=time; + tracktoflen[index]=toflen; + momentum[index]=mom; + switch (abspdg) { + case 211: + truparticle[index]=0; + break ; + case 321: + truparticle[index]=1; + break ; + case 2212: + truparticle[index]=2; + break ; + } + + } + ipartold = ipart; + + if(istop==15){ // start analysis on current set + nset+=1; + lasttrack=track; + istop=0; + selected=0; + cout << "starting t0 calculation for current set" << endl; + for (Int_t i1=0; i1<3;i1++) { + for (Int_t i2=0; i2<3;i2++) { + for (Int_t i3=0; i3<3;i3++) { + for (Int_t i4=0; i4<3;i4++) { + for (Int_t i5=0; i5<3;i5++) { + for (Int_t i6=0; i6<3;i6++) { + for (Int_t i7=0; i7<3;i7++) { + for (Int_t i8=0; i8<3;i8++) { + for (Int_t i9=0; i9<3;i9++) { + for (Int_t i10=0; i10<3;i10++) { + for (Int_t i11=0; i11<3;i11++) { + for (Int_t i12=0; i12<3;i12++) { + for (Int_t i13=0; i13<3;i13++) { + for (Int_t i14=0; i14<3;i14++) { + for (Int_t i15=0; i15<3;i15++) { + beta[0]=momentum[0]/sqrt(massarray[i1]*massarray[i1]+momentum[0]*momentum[0]); + beta[1]=momentum[1]/sqrt(massarray[i2]*massarray[i2]+momentum[1]*momentum[1]); + beta[2]=momentum[2]/sqrt(massarray[i3]*massarray[i3]+momentum[2]*momentum[2]); + beta[3]=momentum[3]/sqrt(massarray[i4]*massarray[i4]+momentum[3]*momentum[3]); + beta[4]=momentum[4]/sqrt(massarray[i5]*massarray[i5]+momentum[4]*momentum[4]); + beta[5]=momentum[5]/sqrt(massarray[i6]*massarray[i6]+momentum[5]*momentum[5]); + beta[6]=momentum[6]/sqrt(massarray[i7]*massarray[i7]+momentum[6]*momentum[6]); + beta[7]=momentum[7]/sqrt(massarray[i8]*massarray[i8]+momentum[7]*momentum[7]); + beta[8]=momentum[8]/sqrt(massarray[i9]*massarray[i9]+momentum[8]*momentum[8]); + beta[9]=momentum[9]/sqrt(massarray[i10]*massarray[i10]+momentum[9]*momentum[9]); + beta[10]=momentum[10]/sqrt(massarray[i11]*massarray[i11]+momentum[10]*momentum[10]); + beta[11]=momentum[11]/sqrt(massarray[i12]*massarray[i12]+momentum[11]*momentum[11]); + beta[12]=momentum[12]/sqrt(massarray[i13]*massarray[i13]+momentum[12]*momentum[12]); + beta[13]=momentum[13]/sqrt(massarray[i14]*massarray[i14]+momentum[13]*momentum[13]); + beta[14]=momentum[14]/sqrt(massarray[i15]*massarray[i15]+momentum[14]*momentum[14]); + + + Float_t meantzero=0.; + Float_t sumAllweights=0.; + for (Int_t itz=0; itz<15;itz++) { + 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 + sqTrackError[itz]=(timeresolutioninns*timeresolutioninns+sqMomError[itz]); // total error for the current track + sumAllweights+=1./sqTrackError[itz]; + // redefining beta, it is useful in order to calculate t zero + beta[itz]*=0.299792; + timezero[itz]=(tracktoflen[itz]/beta[itz])-timeofflight[itz]; + weightedtimezero[itz]=((tracktoflen[itz]/beta[itz])-timeofflight[itz])/sqTrackError[itz];// weighted time zero for current track + meantzero+=weightedtimezero[itz]; + } + meantzero=meantzero/sumAllweights; // it is given in [ns] + + Float_t dummychisquare=0.; + // calculate the chisquare for the current assignment + for (Int_t itz=0; itz<15;itz++) { + dummychisquare+=(timezero[itz]-meantzero)*(timezero[itz]-meantzero)/sqTrackError[itz]; + } + + + if(dummychisquare<=chisquare){ + assparticle[0]=i1; + assparticle[1]=i2; + assparticle[2]=i3; + assparticle[3]=i4; + assparticle[4]=i5; + assparticle[5]=i6; + assparticle[6]=i7; + assparticle[7]=i8; + assparticle[8]=i9; + assparticle[9]=i10; + assparticle[10]=i11; + assparticle[11]=i12; + assparticle[12]=i13; + assparticle[13]=i14; + assparticle[14]=i15; + chisquare=dummychisquare; + t0best=meantzero; + } + } + } + } + } + } + } + } + } + } + } + + } + } + } + } + } + + cout << "true" << truparticle[0] << truparticle[1] << truparticle[2] << truparticle[3] << truparticle[4] << truparticle[5] << truparticle[6] << truparticle[7] << truparticle[8] << truparticle[9] <Fill(t0best,weight); + htzerobest->Fill(t0best); + hchibest->Fill(chisquare); + Double_t dblechisquare=(Double_t)chisquare; + Float_t confLevel=(Float_t)TMath::Prob(dblechisquare,14); + cout << "conf level " << confLevel << endl; + hchibestconflevel->Fill(confLevel); + chisquare=999.; + t0best=999.; + + } // end for the current set. close if(istop==5) + + + } // end condition on ipartold + + + }// end loop on hits for the current track + }// end TOF + if(istop>=15) break; + }// end loop on Tracks + cout << "ngood for current event " << ngood << endl; + + + }//end loop on events + + nmisidentified=(nmisidentified0+nmisidentified1+nmisidentified2+nmisidentified3+nmisidentified4+nmisidentified5+nmisidentified6+nmisidentified7+nmisidentified8+nmisidentified9+nmisidentified10+nmisidentified11+nmisidentified12+nmisidentified13+nmisidentified14); + cout << "total misidentified " << nmisidentified << endl; + TFile *houtfile = new TFile(outFileName,"recreate"); + houtfile->cd(); + + htzerobest->Write(); + hchibest->Write(); + hchibestconflevel->Write(); + houtfile->Close(); + + gBenchmark->Stop("TOFT0"); + cout << "15 tracks:" << endl ; + cout << " took " << gBenchmark->GetCpuTime("TOFT0") << " seconds in order to calculate T0 " + << gBenchmark->GetCpuTime("TOFT0")/nallevents << " seconds per event " << endl ; + cout << endl ; + +} + + -- 2.43.0