]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
LHC10b,c pass3 support
authorantoniol <antoniol@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 5 Mar 2012 08:45:56 +0000 (08:45 +0000)
committerantoniol <antoniol@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 5 Mar 2012 08:45:56 +0000 (08:45 +0000)
ANALYSIS/TenderSupplies/AliTOFTenderSupply.cxx
ANALYSIS/TenderSupplies/AliTOFTenderSupply.h

index 276e667cd09d31c8d9f31714bd70ff96cf21f692..80fb2de0a1300aaa92d832a322e1233ef520cc07 100644 (file)
@@ -23,6 +23,9 @@
 // Contacts: Pietro.Antonioli@bo.infn.it                                     //
 //           Francesco.Noferini@bo.infn.it                                   //
 ///////////////////////////////////////////////////////////////////////////////
 // Contacts: Pietro.Antonioli@bo.infn.it                                     //
 //           Francesco.Noferini@bo.infn.it                                   //
 ///////////////////////////////////////////////////////////////////////////////
+#include <TFile.h>
+#include <TChain.h>
+
 #include <TMath.h>
 #include <TRandom.h>
 #include <AliLog.h>
 #include <TMath.h>
 #include <TRandom.h>
 #include <AliLog.h>
@@ -36,6 +39,7 @@
 #include <AliTOFcalib.h>
 #include <AliTOFT0maker.h>
 
 #include <AliTOFcalib.h>
 #include <AliTOFT0maker.h>
 
+#include <AliGeomManager.h>
 #include <AliCDBManager.h>
 #include <AliCDBEntry.h>
 #include <AliT0CalibSeasonTimeShift.h>
 #include <AliCDBManager.h>
 #include <AliCDBEntry.h>
 #include <AliT0CalibSeasonTimeShift.h>
@@ -53,14 +57,28 @@ AliTOFTenderSupply::AliTOFTenderSupply() :
   fIsMC(kFALSE),
   fTimeZeroType(AliESDpid::kBest_T0),
   fCorrectExpTimes(kTRUE),
   fIsMC(kFALSE),
   fTimeZeroType(AliESDpid::kBest_T0),
   fCorrectExpTimes(kTRUE),
+  fCorrectTRDBug(kFALSE),
   fLHC10dPatch(kFALSE),
   fT0DetectorAdjust(kFALSE),
   fDebugLevel(0),
   fAutomaticSettings(kTRUE),
   fLHC10dPatch(kFALSE),
   fT0DetectorAdjust(kFALSE),
   fDebugLevel(0),
   fAutomaticSettings(kTRUE),
+  fRecoPass(0),
+  fUserRecoPass(0),
   fTOFCalib(0x0),
   fTOFT0maker(0x0),
   fTOFres(100.),
   fTOFCalib(0x0),
   fTOFT0maker(0x0),
   fTOFres(100.),
-  fT0IntercalibrationShift(0)
+  fT0IntercalibrationShift(0),
+  fGeomSet(kFALSE),
+  fIsEnteringInTRD(kFALSE),
+  fInTRD(kFALSE),
+  fIsComingOutTRD(kFALSE),
+  fOutTRD(kFALSE),
+  fRhoTRDin(288.38), // cm
+  fRhoTRDout(366.38), // cm
+  fStep(0.5),
+  fMagField(0.),
+  fCDBkey(0)
+
 
 
 {
 
 
 {
@@ -80,14 +98,27 @@ AliTOFTenderSupply::AliTOFTenderSupply(const char *name, const AliTender *tender
   fIsMC(kFALSE),
   fTimeZeroType(AliESDpid::kBest_T0),
   fCorrectExpTimes(kTRUE),
   fIsMC(kFALSE),
   fTimeZeroType(AliESDpid::kBest_T0),
   fCorrectExpTimes(kTRUE),
+  fCorrectTRDBug(kFALSE),
   fLHC10dPatch(kFALSE),
   fT0DetectorAdjust(kFALSE),
   fDebugLevel(0),
   fAutomaticSettings(kTRUE),
   fLHC10dPatch(kFALSE),
   fT0DetectorAdjust(kFALSE),
   fDebugLevel(0),
   fAutomaticSettings(kTRUE),
+  fRecoPass(0),
+  fUserRecoPass(0),
   fTOFCalib(0x0),
   fTOFT0maker(0x0),
   fTOFres(100.),
   fTOFCalib(0x0),
   fTOFT0maker(0x0),
   fTOFres(100.),
-  fT0IntercalibrationShift(0)
+  fT0IntercalibrationShift(0),
+  fGeomSet(kFALSE),
+  fIsEnteringInTRD(kFALSE),
+  fInTRD(kFALSE),
+  fIsComingOutTRD(kFALSE),
+  fOutTRD(kFALSE),
+  fRhoTRDin(288.38), // cm
+  fRhoTRDout(366.38), // cm
+  fStep(0.5),
+  fMagField(0.),
+  fCDBkey(0)
  
 {
   //
  
 {
   //
@@ -112,12 +143,17 @@ void AliTOFTenderSupply::Init()
   Int_t run = fTender->GetRun();
   if (run == 0) return;                // to skip first init, when we don't have yet a run number
 
   Int_t run = fTender->GetRun();
   if (run == 0) return;                // to skip first init, when we don't have yet a run number
 
+  fGeomSet=kFALSE;                         
+
   if (fAutomaticSettings) {
   if (fAutomaticSettings) {
+    if (fUserRecoPass == 0) DetectRecoPass();
+    else fRecoPass = fUserRecoPass;
     if (run<114737) {
       tenderUnsupported = kTRUE;
     }
     else if (run>=114737&&run<=117223) {      //period="LHC10B";
     if (run<114737) {
       tenderUnsupported = kTRUE;
     }
     else if (run>=114737&&run<=117223) {      //period="LHC10B";
-      fCorrectExpTimes=kTRUE;
+      if (fRecoPass == 2) {fCorrectExpTimes=kTRUE; fCorrectTRDBug=kFALSE;}
+      else if (fRecoPass == 3) {fCorrectExpTimes=kFALSE; fCorrectTRDBug=kTRUE;}
       fLHC10dPatch=kFALSE;
       fTOFres=100.;
       fTimeZeroType=AliESDpid::kTOF_T0;
       fLHC10dPatch=kFALSE;
       fTOFres=100.;
       fTimeZeroType=AliESDpid::kTOF_T0;
@@ -125,7 +161,8 @@ void AliTOFTenderSupply::Init()
       fT0DetectorAdjust=kTRUE;
     }
     else if (run>=118503&&run<=121040) { //period="LHC10C";
       fT0DetectorAdjust=kTRUE;
     }
     else if (run>=118503&&run<=121040) { //period="LHC10C";
-      fCorrectExpTimes=kTRUE;
+      if (fRecoPass == 2) {fCorrectExpTimes=kTRUE; fCorrectTRDBug=kFALSE;}
+      else if (fRecoPass == 3) {fCorrectExpTimes=kFALSE; fCorrectTRDBug=kTRUE;}
       fLHC10dPatch=kFALSE;
       fTOFres=100.;
       fTimeZeroType=AliESDpid::kTOF_T0;
       fLHC10dPatch=kFALSE;
       fTOFres=100.;
       fTimeZeroType=AliESDpid::kTOF_T0;
@@ -201,7 +238,7 @@ void AliTOFTenderSupply::Init()
     AliFatal(" ------- TOF tender not to be used in this run, issuing FATAL error -------- ");
   }
 
     AliFatal(" ------- TOF tender not to be used in this run, issuing FATAL error -------- ");
   }
 
-  // Check if another detector already created the esd pid object
+  // Check if another tender wagon already created the esd pid object
   // if not we create it and set it to the ESD input handler
   fESDpid=fTender->GetESDhandler()->GetESDpid();
   if (!fESDpid) {
   // if not we create it and set it to the ESD input handler
   fESDpid=fTender->GetESDhandler()->GetESDpid();
   if (!fESDpid) {
@@ -221,18 +258,19 @@ void AliTOFTenderSupply::Init()
   fTOFCalib->SetRemoveMeanT0(!(fIsMC));        // must be kFALSE on MC (default is kTRUE)
   fTOFCalib->SetCalibrateTOFsignal(!(fIsMC));  // must be kFALSE on MC (no new calibration) (default is kTRUE)
   fTOFCalib->SetCorrectTExp(fCorrectExpTimes); // apply a fine tuning on the expected times at low momenta
   fTOFCalib->SetRemoveMeanT0(!(fIsMC));        // must be kFALSE on MC (default is kTRUE)
   fTOFCalib->SetCalibrateTOFsignal(!(fIsMC));  // must be kFALSE on MC (no new calibration) (default is kTRUE)
   fTOFCalib->SetCorrectTExp(fCorrectExpTimes); // apply a fine tuning on the expected times at low momenta
-  
+                                               // (this is done for LHC10b/LHC10c pass2)
 
   // Configure TOFT0 maker class
   //  if (!fTOFT0maker) fTOFT0maker = new AliTOFT0maker(fESDpid,fTOFCalib); // create if needed
   if (!fTOFT0maker) fTOFT0maker = new AliTOFT0maker(fESDpid); // without passing AliTOFCalib it uses the diamond
   fTOFT0maker->SetTimeResolution(fTOFres);     // set TOF resolution for the PID
 
   // Configure TOFT0 maker class
   //  if (!fTOFT0maker) fTOFT0maker = new AliTOFT0maker(fESDpid,fTOFCalib); // create if needed
   if (!fTOFT0maker) fTOFT0maker = new AliTOFT0maker(fESDpid); // without passing AliTOFCalib it uses the diamond
   fTOFT0maker->SetTimeResolution(fTOFres);     // set TOF resolution for the PID
-  
+  fTOFT0maker->SetTOFT0algorithm(2);
 
   AliInfo("|******************************************************|");
   AliInfo(Form("|    Alice TOF Tender Initialisation (Run %d)  |",fTender->GetRun()));
   AliInfo("|    Settings:                                         |");
   AliInfo(Form("|    Correct Exp Times              :  %d               |",fCorrectExpTimes));
 
   AliInfo("|******************************************************|");
   AliInfo(Form("|    Alice TOF Tender Initialisation (Run %d)  |",fTender->GetRun()));
   AliInfo("|    Settings:                                         |");
   AliInfo(Form("|    Correct Exp Times              :  %d               |",fCorrectExpTimes));
+  AliInfo(Form("|    Correct TRD Bug                :  %d               |",fCorrectTRDBug));
   AliInfo(Form("|    LHC10d patch                   :  %d               |",fLHC10dPatch));
   AliInfo(Form("|    TOF resolution for TOFT0 maker :  %5.2f (ps)     |",fTOFres));
   AliInfo(Form("|    timeZero selection             :  %d               |",fTimeZeroType));
   AliInfo(Form("|    LHC10d patch                   :  %d               |",fLHC10dPatch));
   AliInfo(Form("|    TOF resolution for TOFT0 maker :  %5.2f (ps)     |",fTOFres));
   AliInfo(Form("|    timeZero selection             :  %d               |",fTimeZeroType));
@@ -291,13 +329,18 @@ void AliTOFTenderSupply::ProcessEvent()
 
 
   fTOFCalib->CalibrateESD(event);   //recalculate TOF signal (no harm for MC, see settings inside init)
 
 
   fTOFCalib->CalibrateESD(event);   //recalculate TOF signal (no harm for MC, see settings inside init)
-  if (fLHC10dPatch && !(fIsMC)) RecomputeTExp(event);
+
+
+  // patches for various reconstruction bugs
+  if (fLHC10dPatch && !(fIsMC)) RecomputeTExp(event);   // LHC10d pass2: fake full TRD geometry
+  if (fCorrectTRDBug && !(fIsMC)) FixTRDBug(event);     // LHC10b,c pass3: wrong TRD dE/dx 
 
   Double_t startTime = 0.;
 
   Double_t startTime = 0.;
-  if(fIsMC) startTime = fTOFCalib->TuneForMC(event,fTOFres);
+  if(fIsMC) startTime = fTOFCalib->TuneForMC(event,fTOFres);   // this is for old MC when we didn't jitter startTime in MC
 
   if (fDebugLevel > 1) Printf(" TofTender: startTime %f",startTime);
   if (fDebugLevel > 1) Printf(" TofTender: T0 time (orig) %f %f %f",event->GetT0TOF(0),event->GetT0TOF(1),event->GetT0TOF(2));
 
   if (fDebugLevel > 1) Printf(" TofTender: startTime %f",startTime);
   if (fDebugLevel > 1) Printf(" TofTender: T0 time (orig) %f %f %f",event->GetT0TOF(0),event->GetT0TOF(1),event->GetT0TOF(2));
+
   // event by event TO detector treatment  
   if(event->GetT0TOF()){   // protection: we adjust T0 only if it is there....
 
   // event by event TO detector treatment  
   if(event->GetT0TOF()){   // protection: we adjust T0 only if it is there....
 
@@ -345,12 +388,18 @@ void AliTOFTenderSupply::ProcessEvent()
   fTOFT0maker->ComputeT0TOF(event);
   fTOFT0maker->WriteInESD(event);
 
   fTOFT0maker->ComputeT0TOF(event);
   fTOFT0maker->WriteInESD(event);
 
-  // recalculate PID probabilities
+  //set preferred startTime
   fESDpid->SetTOFResponse(event, (AliESDpid::EStartTimeType_t)fTimeZeroType);
 
   fESDpid->SetTOFResponse(event, (AliESDpid::EStartTimeType_t)fTimeZeroType);
 
+  // recalculate PID probabilities
   // this is for safety, especially if the user doesn't attach a PID tender after TOF tender  
   Int_t ntracks=event->GetNumberOfTracks();
   // this is for safety, especially if the user doesn't attach a PID tender after TOF tender  
   Int_t ntracks=event->GetNumberOfTracks();
+  //  AliESDtrack *track = NULL;
+  //  Float_t tzeroTrack = 0;
   for(Int_t itrack = 0; itrack < ntracks; itrack++){
   for(Int_t itrack = 0; itrack < ntracks; itrack++){
+    //    track = event->GetTrack(itrack);
+    //    tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
+    //    Printf("================> Track # %d mom: %f tzeroTrack %f",itrack,track->P(),tzeroTrack);
     fESDpid->MakeTOFPID(event->GetTrack(itrack),0);   
   }
   
     fESDpid->MakeTOFPID(event->GetTrack(itrack),0);   
   }
   
@@ -409,3 +458,718 @@ void AliTOFTenderSupply::RecomputeTExp(AliESDtrack *track) const
 }
 
 
 }
 
 
+//______________________________________________________________________________
+void AliTOFTenderSupply::DetectRecoPass()
+{
+  //
+  // Detect reconstruction information
+  //
+  
+  //reset information
+  fRecoPass=0;
+  
+  //Get the current file to check the reconstruction pass (UGLY, but not stored in ESD... )
+  AliAnalysisManager *mgr=AliAnalysisManager::GetAnalysisManager();
+  AliVEventHandler *inputHandler=mgr->GetInputEventHandler();
+  if (!inputHandler) return;
+  
+  TTree *tree= (TTree*)inputHandler->GetTree();
+  TFile *file= (TFile*)tree->GetCurrentFile();
+  
+  if (!file) {
+    AliFatal("Current file not found");
+  }
+  
+  //find pass from file name (UGLY, but not stored in ESD... )
+  TString fileName(file->GetName());
+  if (fileName.Contains("pass1") ) {
+    fRecoPass=1;
+  } else if (fileName.Contains("pass2") ) {
+    fRecoPass=2;
+  } else if (fileName.Contains("pass3") ) {
+    fRecoPass=3;
+  }
+  if (fRecoPass == 0) {
+    AliInfo(Form("From file name %s reco pass cannot be detected",fileName.Data()));
+    AliInfo("Change file name or use SetUserRecoPass method");
+    AliFatal("------------- TOF tender cannot run with reco pass unspecified, issuing FATAL error ---------- ");
+  }
+}
+
+
+//______________________________________________________________________________
+void AliTOFTenderSupply::InitGeom()
+{
+
+  if (fGeomSet == kTRUE) return;
+
+  //  Printf("\n \n ----- calling InitGeom to fix TRD Bug ----- \n \n");
+  AliCDBManager * man = AliCDBManager::Instance();
+  man->SetDefaultStorage("raw://");
+  fCDBkey = man->SetLock(kFALSE, fCDBkey);
+  Int_t run = fTender->GetRun();
+  //  Printf(" ---------> run is %d",run);
+  man->SetRun(run);
+  AliCDBEntry *entry = (AliCDBEntry*)man->Get("GRP/Geometry/Data");
+  if (entry) {
+    AliGeomManager::LoadGeometry();
+    AliGeomManager::ApplyAlignObjsFromCDB("ITS TPC TRD TOF");
+    //    fCDBkey = man->SetLock(kTRUE, fCDBkey);
+    //    Printf("\n \n ----- Geometry loaded ------ \n \n");
+  }
+  fGeomSet=kTRUE;
+
+}
+
+
+//______________________________________________________________________________
+void AliTOFTenderSupply::FixTRDBug(AliESDEvent* event)
+//
+// recompute texp fixing wrong dE/dx from TRD (LHC10b,c pass3)
+//
+{
+
+  if (fGeomSet == kFALSE) InitGeom();
+
+  //  Printf("Running FixTRD bug ");
+  /* loop over tracks */
+  AliESDtrack *track = NULL;
+  for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) {
+    track = event->GetTrack(itrk);
+    FixTRDBug(track);
+  }
+}
+
+
+//_____________________________________________________
+void AliTOFTenderSupply::FixTRDBug(AliESDtrack *track)
+{
+  // 
+  //
+  //
+
+
+    ULong_t status=track->GetStatus();
+    if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) &&
+          ( (status & AliVTrack::kTPCrefit)==AliVTrack::kTPCrefit ) &&
+          ( (status & AliVTrack::kTPCout)==AliVTrack::kTPCout ) &&
+          ( (status & AliVTrack::kTOFout)==AliVTrack::kTOFout ) &&
+          ( (status & AliVTrack::kTIME)==AliVTrack::kTIME ) ) ) return;
+
+    fIsEnteringInTRD=kFALSE;
+    fInTRD=kFALSE;
+    fIsComingOutTRD=kFALSE;
+    fOutTRD=kFALSE;
+
+    //    Printf("Track reached TOF %f",track->P());
+    Double_t correctionTimes[5] = {0.,0.,0.,0.,0.}; // to be added to the expected times
+    FindTRDFix(track, correctionTimes);
+    Double_t expectedTimes[5] = {0.,0.,0.,0.,0.}; track->GetIntegratedTimes(expectedTimes);
+    //    Printf("Exp. times: %f %f %f %f %f",
+    //    expectedTimes[0],expectedTimes[1],expectedTimes[2],expectedTimes[3],expectedTimes[4]);
+    //    Printf("Corr. times: %f %f %f %f %f",
+    //    correctionTimes[0],correctionTimes[1],correctionTimes[2],correctionTimes[3],correctionTimes[4]);
+
+    for (Int_t jj=0; jj<5; jj++) expectedTimes[jj]+=correctionTimes[jj];
+    track->SetIntegratedTimes(expectedTimes);
+  
+}
+
+
+//________________________________________________________________________
+void AliTOFTenderSupply::FindTRDFix(AliESDtrack *track,Double_t *corrections)
+{
+
+  Double_t pT = track->Pt();
+  ULong_t status=track->GetStatus();
+  Bool_t isTRDout = (status & AliVTrack::kTRDout)==AliVTrack::kTRDout;
+
+  Double_t length = 0.;
+
+  Double_t xyzIN[3]={0.,0.,0.};
+  fIsEnteringInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyzIN);
+
+  Double_t xyzOUT[3]={0.,0.,0.};
+  fIsComingOutTRD = track->GetXYZAt(fRhoTRDout,fMagField,xyzOUT);
+
+  if (fIsEnteringInTRD && fIsComingOutTRD) {
+
+
+    Double_t phiIN = TMath::Pi()+TMath::ATan2(-xyzIN[1],-xyzIN[0]);
+    phiIN *= TMath::RadToDeg();
+    fInTRD = ( (phiIN>=  0. && phiIN<= 40.) ||
+              (phiIN>=140. && phiIN<=220.) ||
+              (phiIN>=340. && phiIN<=360.) ); // TRD SMs installed @ 2010
+
+    Double_t phiOUT = TMath::Pi()+TMath::ATan2(-xyzOUT[1],-xyzOUT[0]);
+    phiOUT *= TMath::RadToDeg();
+    fOutTRD = ( (phiOUT>=  0. && phiOUT<= 40.) ||
+               (phiOUT>=140. && phiOUT<=220.) ||
+               (phiOUT>=340. && phiOUT<=360.) ); // TRD SMs installed @ 2010
+
+    length = 0.;
+
+    if (fInTRD || fOutTRD) {
+
+      if ( ( fInTRD && fOutTRD ) || ( fInTRD && !fOutTRD ) ) {
+       length = EstimateLengthInTRD1(track);
+      } else if ( !fInTRD && fOutTRD ) {
+       length = EstimateLengthInTRD2(track);
+      }
+
+    } else { // ( !fInTRD && !fOutTRD )
+
+      length = EstimateLengthOutTRD(track);
+
+    }
+
+  }
+  //  Printf("estimated length in TRD %f [isTRDout %d]",length,isTRDout);
+  CorrectDeltaTimes(pT,length,isTRDout,corrections);
+
+}
+
+//________________________________________________________________________
+void AliTOFTenderSupply::CorrectDeltaTimes(Double_t pT,
+                                               Double_t length,
+                                               Bool_t flagTRDout,
+                                               Double_t *corrections)
+{
+
+  corrections[2] = CorrectExpectedPionTime(pT,length,flagTRDout);
+  corrections[0] = corrections[2]; // x electrons used pion corrections
+  corrections[1] = corrections[2]; // x muons used pion corrections
+  corrections[3] = CorrectExpectedKaonTime(pT,length,flagTRDout);
+  corrections[4] = CorrectExpectedProtonTime(pT,length,flagTRDout);
+
+}
+
+//________________________________________________________________________
+Double_t AliTOFTenderSupply::CorrectExpectedPionTime(Double_t pT,
+                                                         Double_t length,
+                                                         Bool_t isTRDout)
+{
+  // correction for expected time for pions
+
+  Double_t delta=0.;
+
+  //  Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
+  if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
+
+    Float_t p[2]={0.,0.};
+
+    if (isTRDout) {
+
+      if (pT<0.30) {
+       p[0] = 180.; p[1] = 0.;
+      } else if (pT>=0.30 && pT<0.35) {
+       p[0] = 740.; p[1] = -1800.;
+      } else if (pT>=0.35 && pT<0.40) {
+       p[0] = 488.; p[1] =-1080.;
+      } else if (pT>=0.40 && pT<0.45) {
+       p[0] = 179.; p[1] = -307.;
+      } else if (pT>=0.45 && pT<0.50) {
+       p[0] =  97.; p[1] = -123.;
+      } else { //if (pT>=0.50)
+       p[0] = 120.; p[1] = -172.;
+      }
+
+    } else {
+
+      if (pT<0.30) {
+       p[0] =  70.; p[1] =    0.;
+      } else if (pT>=0.30 && pT<0.35) {
+       p[0] = 339.; p[1] = -927.;
+      } else if (pT>=0.35 && pT<0.40) {
+       p[0] =  59.; p[1] = -121.;
+      } else if (pT>=0.40 && pT<0.50) {
+       p[0] =  21.; p[1] =  -24.;
+      } else { //if (pT>=0.50)
+       p[0] =  42.; p[1] =  -67.;
+      }
+
+    }
+
+    delta = p[0]+p[1]*pT;
+    //    Printf("Pion time: %f %f %f %f",p[0],p[1],length,delta);
+
+  } else {
+
+    Float_t p[2] = {0.,0.};
+
+    if ( fInTRD &&  fOutTRD) { // zone 1
+
+      if (isTRDout) {
+
+       if (length<130.) {
+         p[0] = 0.; p[1] = 0.;
+       } else if (length>=130. && length<170.) {
+         p[0] = -20.5; p[1] = 0.25;
+       } else {//if (length>=170.)
+         p[0] = 22.; p[1] = 0.;
+       }
+
+      } else { // !isTRDout
+
+       p[0] = 20.; p[1] = 0.;
+
+      }
+
+    } else if (!fInTRD && !fOutTRD) { // zone 2
+
+      p[0] = 0.; p[1] = 0.;
+
+    } else if ( fInTRD &&  !fOutTRD) { // zone 3
+
+      if (isTRDout) {
+
+       if (length< 75.) {
+         p[0] = 17.; p[1] =  0.;
+       } else if (length>= 75. && length< 95.) {
+         p[0] = 81.; p[1] = -0.85;
+       } else if (length>= 95. && length<155.) {
+         p[0] =  0.; p[1] =  0.;
+       } else {//if (length>=155.)
+         p[0] = 10.; p[1] =  0.;
+       }
+
+      } else { // !isTRDout
+
+       p[0] = 0.; p[1] = 0.;
+
+      }
+
+    } else if (!fInTRD &&  fOutTRD) { // zone 4
+
+      if (isTRDout) {
+
+       if (length<80.) {
+         p[0] =  0.; p[1] =  0.;
+       } else {//if (length>=80.)
+         p[0] = 10.; p[1] =  0.;
+       }
+
+      } else { // !isTRDout
+
+       if (length<30.) {
+         p[0] =  0.; p[1] =  0.;
+       } else {//if (length>=30.)
+         p[0] =  6.; p[1] =  0.;
+       }
+
+      }
+
+    }
+
+    delta = p[0]+p[1]*length;
+    //    Printf("Pion time: %f %f %f %f",p[0],p[1],length,delta);
+
+  }
+
+  return delta;
+
+}
+
+//________________________________________________________________________
+Double_t AliTOFTenderSupply::CorrectExpectedKaonTime(Double_t pT,
+                                                         Double_t length,
+                                                         Bool_t isTRDout)
+{
+  // correction for expected time for kaons
+
+  Double_t delta=0.;
+  //  Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
+
+  if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
+
+    Float_t p[2]={0.,0.};
+
+    if (isTRDout) {
+
+      if (pT<0.4) {
+       p[0] =  900.; p[1] = 0.;
+      } else if (pT>=0.40 && pT<0.45) {
+       p[0] = 3100.; p[1] = -6000.;
+      } else if (pT>=0.45 && pT<0.50) {
+       p[0] = 1660.; p[1] = -2800.;
+      } else if (pT>=0.50 && pT<0.55) {
+       p[0] =  860.; p[1] = -1200.;
+      } else { //if (pT>=0.55)
+       p[0] =  200.; p[1] = 0.;
+      }
+
+    } else {
+
+      if (pT<0.30) {
+       p[0] =   0.; p[1] =    0.;
+      } else if (pT>=0.30 && pT<0.32) {
+       p[0] = 570.; p[1] =    0.;
+      } else if (pT>=0.32 && pT<0.35) {
+       p[0] = 3171.; p[1] = -8133.;
+      } else if (pT>=0.35 && pT<0.40) {
+       p[0] = 1815.; p[1] = -4260.;
+      } else if (pT>=0.40 && pT<0.45) {
+       p[0] =  715.; p[1] =  -1471.;
+      } else if (pT>=0.45 && pT<0.50) {
+       p[0] =  233.; p[1] =  -407.;
+      } else if (pT>=0.50 && pT<0.55) {
+       p[0] =  408.; p[1] =  -752.;
+      } else { //if (pT>=0.55)
+       p[0] =  408.-752.*0.55; p[1] = 0.;
+      }
+
+    }
+
+    delta = p[0]+p[1]*pT;
+    //    Printf("Kaon time: %f %f %f %f",p[0],p[1],length,delta);
+
+  } else {
+
+    Float_t p[2] = {0.,0.};
+
+    if ( fInTRD &&  fOutTRD) { // zone 1
+
+      if (isTRDout) {
+
+       if (length<95.) {
+         p[0] = 20.; p[1] = 0.;
+       } else if (length>=95. && length<195.) {
+         p[0] = -24.0; p[1] = 0.10+0.0041*length;
+       } else {//if (length>=195.)
+         p[0] =  150.; p[1] = 0.;
+       }
+
+      } else { // !isTRDout
+
+       p[0] = 40.; p[1] = 0.;
+
+      }
+
+    } else if (!fInTRD && !fOutTRD) { // zone 2
+
+      p[0] = 0.; p[1] = 0.;
+
+    } else if ( fInTRD &&  !fOutTRD) { // zone 3
+
+      if (isTRDout) {
+
+       if (length< 15.) {
+         p[0] = 180.; p[1] =  0.;
+       } else if (length>= 15. && length< 55.) {
+         p[0] = 215.; p[1] = -2.5;
+       } else {//if (length>=55.)
+         p[0] = 78.; p[1] =  0.;
+       }
+
+      } else { // !isTRDout
+
+       p[0] = 0.; p[1] = 0.;
+
+      }
+
+    } else if (!fInTRD &&  fOutTRD) { // zone 4
+
+      if (isTRDout) {
+
+       if (length< 55.) {
+         p[0] =  0.; p[1] =  0.;
+       } else if (length>= 55. && length<115.) {
+         p[0] = -85.; p[1] = 1.9;
+       } else {//if (length>=115.)
+         p[0] = 100.; p[1] =  0.;
+       }
+
+      } else { // !isTRDout
+
+       p[0] =  0.; p[1] =  0.;
+
+      }
+
+    }
+
+    delta = p[0]+p[1]*length;
+    //    Printf("Kaon time: %f %f %f %f",p[0],p[1],length,delta);
+
+  }
+
+  return delta;
+
+}
+
+//________________________________________________________________________
+Double_t AliTOFTenderSupply::CorrectExpectedProtonTime(Double_t pT,
+                                                           Double_t length,
+                                                           Bool_t isTRDout)
+{
+  // correction for expected time for protons
+
+  Double_t delta=0.;
+  //  Printf("Flags Ent In ComingOut Out %d %d %d %d",fIsEnteringInTRD,fInTRD,fIsComingOutTRD,fOutTRD);
+
+  if (!fIsEnteringInTRD || !fIsComingOutTRD) { // zone 5
+    Float_t p[2]={0.,0.};
+
+
+    if (isTRDout) {
+
+      if (pT<0.375) {
+       p[0] = 1000.; p[1] = 0.;
+      }        else if (pT>=0.375 && pT<0.45) {
+       p[0] = 1500.; p[1] = 0.;
+      } else if (pT>=0.45 && pT<0.50) {
+       p[0] = 4650.; p[1] = -7000.;
+      } else if (pT>=0.50 && pT<0.55) {
+       p[0] = 3150.; p[1] = -4000.;
+      } else { //if (pT>=0.55)
+       p[0] = 3150. -4000.*0.55; p[1] = 0.;
+      }
+
+    } else {
+
+      if (pT<0.32) {
+       p[0] = 2963.-5670.*0.032; p[1] = 0.;
+      } else if (pT>=0.32 && pT<0.35) {
+       p[0] = 2963.; p[1] =  -5670.;
+      } else if (pT>=0.35 && pT<0.40) {
+       p[0] = 4270.; p[1] =  -9400.;
+      } else if (pT>=0.40 && pT<0.45) {
+       p[0] = 1550.; p[1] =  -2600.;
+      } else if (pT>=0.45 && pT<0.50) {
+       p[0] = 1946.; p[1] =  -3480.;
+      } else if (pT>=0.50 && pT<0.55) {
+       p[0] = 1193.; p[1] =  -1974.;
+      } else { //if (pT>=0.55)
+       p[0] = 1193.-1974.*0.55; p[1] = 0.;
+      }
+
+    }
+
+    delta = p[0]+p[1]*pT;
+    //    Printf("Proton time: %f %f %f %f",p[0],p[1],length,delta);
+
+  } else {
+
+    Float_t p[2] = {0.,0.};
+
+    if ( fInTRD &&  fOutTRD) { // zone 1
+
+      if (isTRDout) {
+
+       if (length<90.) {
+         p[0] = 0.; p[1] = 0.;
+       } else if (length>=90. && length<200.) {
+         p[0] = 1063.; p[1] = -32.+0.30*length-0.00072*length*length;
+       } else {//if (length>=200.)
+         p[0] =  900.; p[1] = 0.;
+       }
+
+      } else { // !isTRDout
+
+       p[0] = 80.; p[1] = 0.;
+
+      }
+
+    } else if (!fInTRD && !fOutTRD) { // zone 2
+
+      if (isTRDout) {
+       p[0] = 0.; p[1] = 0.;
+      } else {
+       if (length<125.) {
+         p[0] = 0.; p[1] = 0.;
+       } else if (length>=125. && length<180.) {
+         p[0] = -132.; p[1] = 1.3;
+       } else {
+         p[0] = 100.; p[1] = 0.;
+       }
+
+      }
+
+    } else if ( fInTRD &&  !fOutTRD) { // zone 3
+
+      if (isTRDout) {
+
+       if (length< 30.) {
+         p[0] = 670.; p[1] =  0.;
+       } else if (length>= 30. && length<155.) {
+         p[0] = 944.; p[1] = -11.+0.064*length;
+       } else {//if (length>=155.)
+         p[0] = 780.; p[1] =  0.;
+       }
+
+      } else { // !isTRDout
+
+       if (length< 30.) {
+         p[0] = 140.; p[1] = -4.5;
+       } else {
+         p[0] = 0.; p[1] = 0.;
+       }
+      
+      }
+
+    } else if (!fInTRD &&  fOutTRD) { // zone 4
+
+      if (isTRDout) {
+
+       if (length< 45.) {
+         p[0] = 130.; p[1] =  0.;
+       } else if (length>= 45. && length<120.) {
+         p[0] = -190.; p[1] = 6.5;
+       } else {//if (length>=120.)
+         p[0] = 750.; p[1] =  0.;
+       }
+
+      } else { // !isTRDout
+
+       if (length<75.5) {
+         p[0] =    0.; p[1] =  0.;
+       } else if (length>= 75.5 && length<90.) {
+         p[0] = -830.; p[1] = 11.;
+       } else {
+         p[0] =  160.; p[1] =  0.;
+       }
+
+      }
+
+    }
+
+    delta = p[0]+p[1]*length;
+    //    Printf("Proton time: %f %f %f %f",p[0],p[1],length,delta);
+
+  }
+
+  return delta;
+
+}
+
+//________________________________________________________________________
+Double_t AliTOFTenderSupply::EstimateLengthInTRD1(AliESDtrack *track)
+{
+
+  Double_t xyz0[3]={0.,0.,0.};
+  Bool_t stayInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyz0);
+
+  Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
+  phi0 *= TMath::RadToDeg();
+  stayInTRD = stayInTRD && ( (phi0>=  0. && phi0<= 40.) ||
+                            (phi0>=140. && phi0<=220.) ||
+                            (phi0>=340. && phi0<=360.) );
+
+  Double_t trackLengthInTRD = 0.;
+  Int_t iStep=0;
+
+  Double_t b[3];track->GetBxByBz(b);
+
+  Double_t xyz1[3]={0.,0.,0.};
+  Double_t rho = fRhoTRDin;
+  while (stayInTRD && rho<=fRhoTRDout) {
+    iStep++;
+    rho += fStep;
+
+    for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
+    stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
+    Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
+    phi1 *= TMath::RadToDeg();
+    stayInTRD = stayInTRD && ( (phi1>=  0. && phi1<= 40.) ||
+                              (phi1>=140. && phi1<=220.) ||
+                              (phi1>=340. && phi1<=360.) );
+
+    Double_t l2  = TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0]) +
+                              (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]) +
+                              (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
+    trackLengthInTRD += l2;
+
+    for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
+  }
+
+  return trackLengthInTRD;
+
+}
+
+//________________________________________________________________________
+Double_t AliTOFTenderSupply::EstimateLengthInTRD2(AliESDtrack *track)
+{
+
+  Double_t xyz0[3]={0.,0.,0.};
+  Bool_t stayInTRD = track->GetXYZAt(fRhoTRDout,fMagField,xyz0);
+
+  Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
+  phi0 *= TMath::RadToDeg();
+  stayInTRD = stayInTRD && ( (phi0>=  0. && phi0<= 40.) ||
+                            (phi0>=140. && phi0<=220.) ||
+                            (phi0>=340. && phi0<=360.) );
+
+  Double_t trackLengthInTRD = 0.;
+  Int_t iStep=0;
+
+  Double_t b[3];track->GetBxByBz(b);
+
+  Double_t xyz1[3]={0.,0.,0.};
+  Double_t rho = fRhoTRDout;
+  while (stayInTRD && rho>=fRhoTRDin) {
+    iStep++;
+    rho -= fStep;
+
+    for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
+    stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
+    Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
+    phi1 *= TMath::RadToDeg();
+    stayInTRD = stayInTRD && ( (phi1>=  0. && phi1<= 40.) ||
+                              (phi1>=140. && phi1<=220.) ||
+                              (phi1>=340. && phi1<=360.) );
+
+    Double_t l2  = TMath::Sqrt((xyz0[0]-xyz1[0])*(xyz0[0]-xyz1[0]) +
+                              (xyz0[1]-xyz1[1])*(xyz0[1]-xyz1[1]) +
+                              (xyz0[2]-xyz1[2])*(xyz0[2]-xyz1[2]));
+    trackLengthInTRD += l2;
+
+    for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
+  }
+
+  return trackLengthInTRD;
+
+}
+
+//________________________________________________________________________
+Double_t AliTOFTenderSupply::EstimateLengthOutTRD(AliESDtrack *track)
+{
+
+  Double_t xyz0[3]={0.,0.,0.};
+  Bool_t stayInTRD = track->GetXYZAt(fRhoTRDin,fMagField,xyz0);
+
+  Double_t phi0 = TMath::Pi()+TMath::ATan2(-xyz0[1],-xyz0[0]);
+  phi0 *= TMath::RadToDeg();
+  stayInTRD = stayInTRD && !( (phi0>=  0. && phi0<= 40.) ||
+                             (phi0>=140. && phi0<=220.) ||
+                             (phi0>=340. && phi0<=360.) );
+
+  Double_t trackLengthInTRD = 0.;
+  Int_t iStep=0;
+
+  Double_t b[3];track->GetBxByBz(b);
+
+  Double_t xyz1[3]={0.,0.,0.};
+  Double_t rho = fRhoTRDin;
+  while (stayInTRD && rho<=fRhoTRDout) {
+    iStep++;
+    rho += fStep;
+
+    for (Int_t ii=0; ii<3; ii++) xyz1[ii]=0.;
+    stayInTRD = track->GetXYZAt(rho,fMagField,xyz1);
+    Double_t phi1 = TMath::Pi()+TMath::ATan2(-xyz1[1],-xyz1[0]);
+    phi1 *= TMath::RadToDeg();
+    stayInTRD = stayInTRD && !( (phi1>=  0. && phi1<= 40.) ||
+                               (phi1>=140. && phi1<=220.) ||
+                               (phi1>=340. && phi1<=360.) );
+
+    Double_t l2  = TMath::Sqrt((xyz1[0]-xyz0[0])*(xyz1[0]-xyz0[0]) +
+                              (xyz1[1]-xyz0[1])*(xyz1[1]-xyz0[1]) +
+                              (xyz1[2]-xyz0[2])*(xyz1[2]-xyz0[2]));
+    trackLengthInTRD += l2;
+
+    for (Int_t ii=0; ii<3; ii++) xyz0[ii]=xyz1[ii];
+  }
+
+  return trackLengthInTRD;
+
+}
+
index b6e1f48eaad1dd29d3c8033b820817a08900f64c..4d6a286cec007ccc68dfb8ac6c854fee4450aace 100644 (file)
@@ -45,6 +45,9 @@ public:
     return;
   }
   void SetAutomaticSettings(Bool_t flag=kTRUE){fAutomaticSettings=flag;}
     return;
   }
   void SetAutomaticSettings(Bool_t flag=kTRUE){fAutomaticSettings=flag;}
+  void SetUserRecoPass(Int_t flag=0){fUserRecoPass=flag;}
+  Int_t GetRecoPass(void){return fRecoPass;}
+  void DetectRecoPass();
   virtual void SetTimeZeroType(AliESDpid::EStartTimeType_t tofTimeZeroType) {fTimeZeroType = tofTimeZeroType;}
 
   /* theoretical expected time: related stuff for LHC10d patch */
   virtual void SetTimeZeroType(AliESDpid::EStartTimeType_t tofTimeZeroType) {fTimeZeroType = tofTimeZeroType;}
 
   /* theoretical expected time: related stuff for LHC10d patch */
@@ -52,18 +55,33 @@ public:
   static Float_t GetExpTimeTh(Float_t m, Float_t p, Float_t L) {return L / 2.99792457999999984e-02 / GetBetaTh(m, p);}; // get exp time th
   void RecomputeTExp(AliESDEvent *event) const;
   void RecomputeTExp(AliESDtrack *track) const;
   static Float_t GetExpTimeTh(Float_t m, Float_t p, Float_t L) {return L / 2.99792457999999984e-02 / GetBetaTh(m, p);}; // get exp time th
   void RecomputeTExp(AliESDEvent *event) const;
   void RecomputeTExp(AliESDtrack *track) const;
+  void FixTRDBug(AliESDEvent *event);
+  void FixTRDBug(AliESDtrack *track);
+  void InitGeom();
+  void FindTRDFix(AliESDtrack *track,Double_t *corr);
+  Double_t EstimateLengthInTRD1(AliESDtrack *track);
+  Double_t EstimateLengthInTRD2(AliESDtrack *track);
+  Double_t EstimateLengthOutTRD(AliESDtrack *track);
+  void CorrectDeltaTimes(Double_t pT, Double_t length, Bool_t isTRDout, Double_t *corrections);
+  Double_t CorrectExpectedProtonTime(Double_t pT,Double_t length, Bool_t isTRDout);
+  Double_t CorrectExpectedKaonTime(Double_t pT,Double_t length, Bool_t isTRDout);
+  Double_t CorrectExpectedPionTime(Double_t pT,Double_t length, Bool_t isTRDout);
+
 
 private:
   AliESDpid          *fESDpid;         //! ESD pid object
 
 
 private:
   AliESDpid          *fESDpid;         //! ESD pid object
 
+  
   Bool_t fIsMC;              // flag for MC data
   Int_t  fTimeZeroType;      // flag to select timeZero type 
   Bool_t fCorrectExpTimes;   // flag to apply Expected Time correction 
   Bool_t fIsMC;              // flag for MC data
   Int_t  fTimeZeroType;      // flag to select timeZero type 
   Bool_t fCorrectExpTimes;   // flag to apply Expected Time correction 
+  Bool_t fCorrectTRDBug;     // flag to fix wrong dE/dx inside TRD
   Bool_t fLHC10dPatch;       // flag to apply special patch for LHC10d (reconstructed with wrong geometry)
   Bool_t fT0DetectorAdjust;  // flag to apply offsets to T0 data (works only on some periods)
   Int_t  fDebugLevel;        // debug purposes 0= no output, 1 Info, 2 lot of info....
   Bool_t fAutomaticSettings; // enable/disable automatic (per run) settings
   Bool_t fLHC10dPatch;       // flag to apply special patch for LHC10d (reconstructed with wrong geometry)
   Bool_t fT0DetectorAdjust;  // flag to apply offsets to T0 data (works only on some periods)
   Int_t  fDebugLevel;        // debug purposes 0= no output, 1 Info, 2 lot of info....
   Bool_t fAutomaticSettings; // enable/disable automatic (per run) settings
-
+  Int_t  fRecoPass;          // reconstruction pass: the tender applies different recipes depending on the pass
+  Int_t  fUserRecoPass;      // when reco pass is selected by user
 
   // variables for TOF calibrations and timeZero setup
   AliTOFcalib     *fTOFCalib;       // recalibrate TOF signal with OCDB
 
   // variables for TOF calibrations and timeZero setup
   AliTOFcalib     *fTOFCalib;       // recalibrate TOF signal with OCDB
@@ -76,11 +94,22 @@ private:
   static Float_t fgT0Aresolution;   // T0 resolution A-Side (MC)
   static Float_t fgT0Cresolution;   // T0 resolution C-Side (MC)
 
   static Float_t fgT0Aresolution;   // T0 resolution A-Side (MC)
   static Float_t fgT0Cresolution;   // T0 resolution C-Side (MC)
 
+  // variables to steer TRD bug fix
+  Bool_t fGeomSet;                 // steer loading GRP entry
+  Bool_t fIsEnteringInTRD;
+  Bool_t fInTRD;
+  Bool_t fIsComingOutTRD;
+  Bool_t fOutTRD;
+  Float_t fRhoTRDin;                // cm
+  Float_t fRhoTRDout;               // cm
+  Float_t fStep;                    // cm
+  Double_t fMagField;               // magnetic field value [kGauss]
+  ULong_t fCDBkey;
 
   AliTOFTenderSupply(const AliTOFTenderSupply&c);
   AliTOFTenderSupply& operator= (const AliTOFTenderSupply&c);
 
 
   AliTOFTenderSupply(const AliTOFTenderSupply&c);
   AliTOFTenderSupply& operator= (const AliTOFTenderSupply&c);
 
-  ClassDef(AliTOFTenderSupply, 7);
+  ClassDef(AliTOFTenderSupply, 8);
 };
 
 
 };