]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TOF/AliTOFT0maker.cxx
Separate trees for ITSsa and global tracks
[u/mrichter/AliRoot.git] / TOF / AliTOFT0maker.cxx
index 3ee891b925bf8bb71b84c2907bd479b710df17dc..88d542f768e2e8c38344fe12de8917505bd5a143 100644 (file)
@@ -46,6 +46,7 @@
 #include "AliTOFT0v1.h"
 #include "AliTOFT0maker.h"
 #include "AliPID.h"
+#include "AliLog.h"
 #include "AliESDpid.h"
 #include "AliESDEvent.h"
 #include "TFile.h"
@@ -53,6 +54,7 @@
 #include "AliTOFcalib.h"
 #include "AliTOFRunParams.h"
 #include "TRandom.h"
+#include "AliTOFHeader.h"
 
 ClassImp(AliTOFT0maker)
            
@@ -71,7 +73,8 @@ AliTOFT0maker::AliTOFT0maker():
   fKmask(0),
   fT0width(150.),
   fT0spreadExt(-1.),
-  fT0fillExt(0)
+  fT0fillExt(0),
+  fTOFT0algorithm(1)
 {
   // ctr
   fCalculated[0] = 0;
@@ -79,6 +82,9 @@ AliTOFT0maker::AliTOFT0maker():
   fCalculated[2] = 0;
   fCalculated[3] = 0;
 
+  fT0cur[0]=0.;
+  fT0cur[1]=0.;
+
   if(AliPID::ParticleMass(0) == 0) new AliPID();
 
   fPIDesd = new AliESDpid();
@@ -104,7 +110,8 @@ AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
     fKmask(0),
     fT0width(150.),
     fT0spreadExt(-1.),
-    fT0fillExt(0)
+    fT0fillExt(0),
+    fTOFT0algorithm(1)
 {
   // ctr
   fCalculated[0] = 0;
@@ -112,6 +119,9 @@ AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
   fCalculated[2] = 0;
   fCalculated[3] = 0;
 
+  fT0cur[0]=0.;
+  fT0cur[1]=0.;
+
   if(AliPID::ParticleMass(0) == 0) new AliPID();
 
   if(!fPIDesd){
@@ -130,7 +140,6 @@ AliTOFT0maker::AliTOFT0maker(AliESDpid *externalPID, AliTOFcalib *tofCalib):
 AliTOFT0maker::~AliTOFT0maker()
 {
   // dtor
-
   delete fT0TOF;
   if (!fExternalPIDFlag) delete fPIDesd;
 }
@@ -139,7 +148,6 @@ Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t
   //
   // Remake TOF PID probabilities
   //
-
   Double_t t0tof[6];
 
   if(fKmask) ApplyMask(esd);
@@ -164,10 +172,23 @@ Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t
       t0fill = fT0fillExt;
     }
   }
+  else if(esd){
+      Float_t t0spread = esd->GetSigma2DiamondZ(); // vertex pread ^2
+      if(t0spread > 0) t0spread = TMath::Sqrt(t0spread)/0.0299792458;
 
-  fT0TOF->Init(esd);
-  AliTOFT0v1* t0maker= fT0TOF;
+      if(fT0spreadExt > 0) SetT0FillWidth(fT0spreadExt);
+      else{
+         SetT0FillWidth(t0spread);
+         t0fill = fT0fillExt;
+      }
+  }
+
+  Float_t thrGood = TMath::Max(Float_t(500.),fT0width*3);
 
+
+  fT0TOF->Init(esd);
+  AliTOFT0v1* t0maker = fT0TOF;
+  if (fTOFT0algorithm==2) t0maker->SetOptimization(kTRUE);
   t0maker->DefineT0("all",1.5,3.0);
   t0tof[0] = t0maker->GetResult(0);
   t0tof[1] = t0maker->GetResult(1);
@@ -196,11 +217,13 @@ Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t
   fCalculated[6]=sigmaFill; // sigma t0 fill
   fCalculated[7] = t0tof[3];  // n TOF tracks used for T0
 
+  if(fCalculated[7] > 30) thrGood = 10000000;
+
   //statistics
   fCalculated[8] = t0tof[4]; // real time in s
   fCalculated[9] = t0tof[5]; // cpu time in s
 
-  if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < 500 && fCalculated[1] < fTimeResolution*1.2){
+  if(fCalculated[1] < sigmaFill && TMath::Abs(fCalculated[0] - t0fill) < thrGood && fCalculated[1] < fTimeResolution*1.2){
     fT0sigma=fCalculated[1];
     lT0Current=fCalculated[0];
   }
@@ -231,7 +254,7 @@ Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t
     }
   }
 
-  if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < 500){
+  if(fT0sigma < sigmaFill && TMath::Abs(lT0Current - t0fill) < thrGood){
     fCalculated[1]=fT0sigma;
     fCalculated[0]=lT0Current;
   }
@@ -245,6 +268,8 @@ Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t
   }
 
   // T0 pt bin
+  Float_t *t0values = new Float_t[fNmomBins];
+  Float_t *t0resolution = new Float_t[fNmomBins];
   if(fCalculated[7] < 100){
     for(Int_t i=0;i<fNmomBins;i++){
       t0maker->DefineT0("all",fPIDesd->GetTOFResponse().GetMinMom(i),fPIDesd->GetTOFResponse().GetMaxMom(i));
@@ -252,12 +277,11 @@ Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t
       t0tof[1] = t0maker->GetResult(1);
       t0tof[2] = t0maker->GetResult(2);
       t0tof[3] = t0maker->GetResult(3);
 
       Float_t t0bin =-1000*t0tof[0]; // best t0
       Float_t t0binRes =1000*t0tof[1]; // sigma best t0
       
-      if(t0binRes < sigmaFill  && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < 500){
+      if(t0binRes < sigmaFill  && t0binRes < fTimeResolution * 1.2 && TMath::Abs(t0bin - t0fill) < thrGood){
        // Ok T0
        if(t0sigma < 1000){
          Double_t w1 = 1./t0sigma/t0sigma;
@@ -277,16 +301,23 @@ Double_t* AliTOFT0maker::ComputeT0TOF(AliESDEvent *esd,Double_t t0time,Double_t
           t0binRes= t0sigma;     
        }
       }
-      fPIDesd->GetTOFResponse().SetT0bin(i,t0bin);
-      fPIDesd->GetTOFResponse().SetT0binRes(i,t0binRes);
+      t0values[i] = t0bin;
+      t0resolution[i] = t0binRes;
     }
   }
   else{
     for(Int_t i=0;i<fNmomBins;i++){
-      fPIDesd->GetTOFResponse().SetT0bin(i,lT0Current);
-      fPIDesd->GetTOFResponse().SetT0binRes(i,fT0sigma);
+      t0values[i] = lT0Current;
+      t0resolution[i] = fT0sigma;
     }
   }
+  for(Int_t i=0;i<fNmomBins;i++){
+    fPIDesd->GetTOFResponse().SetT0bin(i,t0values[i]);
+    fPIDesd->GetTOFResponse().SetT0binRes(i,t0resolution[i]);
+  }
+  
+  delete[] t0values;
+  delete[] t0resolution;
 
   return fCalculated;
 }
@@ -338,7 +369,7 @@ void AliTOFT0maker::ApplyT0TOF(AliESDEvent *esd){
   //
 }
 //____________________________________________________________________________ 
-void  AliTOFT0maker::LoadChannelMap(char *filename){
+void  AliTOFT0maker::LoadChannelMap(const char *filename){
   // Load the histo with the channel off map
   TFile *f= new TFile(filename);
   if(!f){
@@ -402,12 +433,14 @@ AliTOFT0maker::TuneForMC(AliESDEvent *esd){ // return true T0 event
     if ((t->GetStatus()&AliESDtrack::kTOFout)==0) continue;
     
     /* check if channel is enabled */
-    if (fTOFcalib && !fTOFcalib->IsChannelEnabled(t->GetTOFCalChannel())) {
-      /* reset TOF status */
-      t->ResetStatus(AliESDtrack::kTOFin);
-      t->ResetStatus(AliESDtrack::kTOFout);
-      t->ResetStatus(AliESDtrack::kTOFmismatch);
-      t->ResetStatus(AliESDtrack::kTOFpid);
+    if (fTOFcalib){
+       if(!fTOFcalib->IsChannelEnabled(t->GetTOFCalChannel())) {
+           /* reset TOF status */
+           t->ResetStatus(AliESDtrack::kTOFin);
+           t->ResetStatus(AliESDtrack::kTOFout);
+           t->ResetStatus(AliESDtrack::kTOFmismatch);
+           t->ResetStatus(AliESDtrack::kTOFpid);
+       }
     }
 
     Double_t time=t->GetTOFsignal();
@@ -424,3 +457,91 @@ AliTOFT0maker::TuneForMC(AliESDEvent *esd){ // return true T0 event
   //
   return t0;
 }
+//_________________________________________________________________________
+void     AliTOFT0maker::WriteInESD(AliESDEvent *esd){
+  //
+  // Write t0Gen, t0ResGen, nt0;
+  //       t0resESD[0:nt0], it0ESD[0:nt0]
+  // in the AliESDEvent
+  //
+    Int_t nMomBins = fPIDesd->GetTOFResponse().GetNmomBins();
+
+    Int_t nt0=0;
+    Float_t *t0 = new Float_t[nMomBins];
+    Float_t *t0res = new Float_t[nMomBins];
+    Int_t *multT0 = new Int_t[nMomBins];
+
+    for(Int_t i=0;i<nMomBins;i++){
+//     printf("START %i) %f %f\n",i,fT0event[i],fT0resolution[i]);
+       multT0[i]=0;
+       Bool_t kNewT0=kTRUE;
+       for(Int_t j=0;j < nt0;j++){
+           if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0[j])<0.1){
+               kNewT0=kFALSE;
+               multT0[j]++;
+               j=nMomBins*10;
+           }
+       }
+       if(kNewT0){
+           t0[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
+           t0res[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
+           nt0++;
+       }
+    }
+
+    Int_t iMultT0=0,nmult=0;
+    for(Int_t j=0;j < nt0;j++){ // find the most frequent
+       if(multT0[j] > nmult){
+           iMultT0 = j;
+           nmult = multT0[j];
+       }
+    }
+
+    Float_t *t0ESD = new Float_t[nMomBins];
+    Float_t *t0resESD = new Float_t[nMomBins];
+    Int_t *it0ESD = new Int_t[nMomBins];
+    
+    Float_t t0Gen,t0ResGen;
+    t0Gen = t0[iMultT0];
+    t0ResGen = t0res[iMultT0];
+    nt0=0;
+    //  printf("T0 to ESD\n%f %f\n",t0Gen,t0ResGen);
+    for(Int_t i=0;i<nMomBins;i++){
+       if(TMath::Abs(fPIDesd->GetTOFResponse().GetT0bin(i) - t0Gen)>0.1){
+           t0ESD[nt0]=fPIDesd->GetTOFResponse().GetT0bin(i);
+           t0resESD[nt0]=fPIDesd->GetTOFResponse().GetT0binRes(i);
+           it0ESD[nt0]=i;
+//         printf("%i) %f %f %i\n",nt0,t0ESD[nt0],t0resESD[nt0],it0ESD[nt0]);
+           nt0++;
+       }
+    }
+
+    // Write  t0Gen,t0ResGen; nt0; t0resESD[0:nt0],it0ESD[0:nt0] in the AliESDEvent
+
+   AliTOFHeader *tofHeader =
+       new AliTOFHeader(t0Gen,t0ResGen,nt0,
+                        t0ESD,t0resESD,it0ESD,fTimeResolution,fT0width);
+  
+    esd->SetTOFHeader(tofHeader);
+
+    delete tofHeader;
+
+    AliDebug(1,Form("resTOF=%f T0spread=%f t0Gen=%f t0resGen=%f",fTimeResolution,fT0width,t0Gen,t0ResGen));
+    AliDebug(1,Form("%d ",nt0));
+    for (Int_t ii=0; ii<nt0; ii++)
+      AliDebug(1,Form("pBin=%d t0val=%f t0res=%f",it0ESD[ii],t0ESD[ii],t0resESD[ii]));
+    
+    delete[] t0;
+    t0 = NULL;
+    delete[] t0res;
+    t0res = NULL;
+    delete[] multT0;
+    multT0 = NULL;
+    delete[] t0ESD;
+    t0ESD = NULL;
+    delete[] t0resESD;
+    t0resESD = NULL;
+    delete[] it0ESD;
+    it0ESD = NULL;
+
+}