PWG2/RESONANCES/AliRsnCutESD2010.cxx

   1 //
   2 // Class AliRsnCutESD2010
   3 //
   4 // General implementation of a single cut strategy, which can be:
   5 // - a value contained in a given interval  [--> IsBetween()   ]
   6 // - a value equal to a given reference     [--> MatchesValue()]
   7 //
   8 // In all cases, the reference value(s) is (are) given as data members
   9 // and each kind of cut requires a given value type (Int, UInt, Double),
  10 // but the cut check procedure is then automatized and chosen thanks to
  11 // an enumeration of the implemented cut types.
  12 // At the end, the user (or any other point which uses this object) has
  13 // to use the method IsSelected() to check if this cut has been passed.
  14 //
  15 // authors: Martin Vala (martin.vala@cern.ch)
  16 //          Alberto Pulvirenti (alberto.pulvirenti@ct.infn.it)
  17 //
  18
  19 #include <Riostream.h>
  20
  21 #include "AliESDpid.h"
  22 #include "AliTOFT0maker.h"
  23 #include "AliTOFcalib.h"
  24 #include "AliCDBManager.h"
  25 #include "AliITSPIDResponse.h"
  26
  27 #include "AliRsnEvent.h"
  28 #include "AliRsnDaughter.h"
  29 #include "AliRsnCutESD2010.h"
  30
  31 ClassImp(AliRsnCutESD2010)
  32
  33 //_________________________________________________________________________________________________
  34 AliRsnCutESD2010::AliRsnCutESD2010
  35 (const char *name, Bool_t isMC) :
  36   AliRsnCut(name, AliRsnCut::kDaughter, 0.0, 0.0),
  37   fIsMC(isMC),
  38   fCheckITS(kTRUE),
  39   fCheckTPC(kTRUE),
  40   fCheckTOF(kTRUE),
  41   fUseGlobal(kTRUE),
  42   fUseITSSA(kTRUE),
  43   fPID(AliPID::kKaon),
  44   fMaxEta(1E6),
  45   fMaxITSband(3.0),
  46   fTPCpLimit(0.35),
  47   fMinTPCband(3.0),
  48   fMaxTPCband(5.0),
  49   fESDtrackCutsTPC(),
  50   fESDtrackCutsITS(),
  51   fESDpid(0x0),
  52   fTOFmaker(0x0),
  53   fTOFcalib(0x0),
  54   fTOFcalibrateESD(!isMC),
  55   fTOFcorrectTExp(kTRUE),
  56   fTOFuseT0(kTRUE),
  57   fTOFtuneMC(isMC),
  58   fTOFresolution(100.0),
  59   fMinTOF(-2.5),
  60   fMaxTOF( 3.5),
  61   fLastRun(-1)
  62 {
  63 //
  64 // Main constructor.
  65 //
  66
  67   SetMC(isMC);
  68
  69   // set default quality cuts for TPC+ITS tracks
  70   // TPC
  71   fESDtrackCutsTPC.SetRequireTPCStandAlone(kTRUE); // to get chi2 and ncls of kTPCin
  72   fESDtrackCutsTPC.SetMinNClustersTPC(70);
  73   fESDtrackCutsTPC.SetMaxChi2PerClusterTPC(4);
  74   fESDtrackCutsTPC.SetAcceptKinkDaughters(kFALSE);
  75   fESDtrackCutsTPC.SetRequireTPCRefit(kTRUE);
  76   // ITS
  77   fESDtrackCutsTPC.SetRequireITSRefit(kTRUE);
  78   fESDtrackCutsTPC.SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
  79   fESDtrackCutsTPC.SetMaxDCAToVertexXYPtDep("0.0350+0.0490/pt^1.0");
  80   fESDtrackCutsTPC.SetMaxDCAToVertexZ(1.e6);
  81   fESDtrackCutsTPC.SetDCAToVertex2D(kFALSE);
  82   fESDtrackCutsTPC.SetRequireSigmaToVertex(kFALSE);
  83   fESDtrackCutsTPC.SetEtaRange(-fMaxEta, fMaxEta);
  84
  85   // set default quality cuts for ITS standalone tracks
  86   fESDtrackCutsITS.SetRequireITSStandAlone(kTRUE);
  87   fESDtrackCutsITS.SetRequireITSPureStandAlone(kFALSE);
  88   fESDtrackCutsITS.SetRequireITSRefit(kTRUE);
  89   fESDtrackCutsITS.SetMinNClustersITS(4);
  90   fESDtrackCutsITS.SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny);
  91   fESDtrackCutsITS.SetMaxChi2PerClusterITS(2.5);
  92   fESDtrackCutsITS.SetMaxDCAToVertexXYPtDep("0.0595+0.0182/pt^1.55");
  93   fESDtrackCutsITS.SetEtaRange(-fMaxEta, fMaxEta);
  94 }
  95
  96 //_________________________________________________________________________________________________
  97 AliRsnCutESD2010::AliRsnCutESD2010
  98 (const AliRsnCutESD2010& copy) :
  99   AliRsnCut(copy),
 100   fIsMC(copy.fIsMC),
 101   fCheckITS(copy.fCheckITS),
 102   fCheckTPC(copy.fCheckTPC),
 103   fCheckTOF(copy.fCheckTOF),
 104   fUseGlobal(copy.fUseGlobal),
 105   fUseITSSA(copy.fUseITSSA),
 106   fPID(copy.fPID),
 107   fMaxEta(copy.fMaxEta),
 108   fMaxITSband(copy.fMaxITSband),
 109   fTPCpLimit(copy.fTPCpLimit),
 110   fMinTPCband(copy.fMinTPCband),
 111   fMaxTPCband(copy.fMaxTPCband),
 112   fESDtrackCutsTPC(copy.fESDtrackCutsTPC),
 113   fESDtrackCutsITS(copy.fESDtrackCutsITS),
 114   fESDpid(0x0),
 115   fTOFmaker(0x0),
 116   fTOFcalib(0x0),
 117   fTOFcalibrateESD(copy.fTOFcalibrateESD),
 118   fTOFcorrectTExp(copy.fTOFcorrectTExp),
 119   fTOFuseT0(copy.fTOFuseT0),
 120   fTOFtuneMC(copy.fTOFtuneMC),
 121   fTOFresolution(copy.fTOFresolution),
 122   fMinTOF(copy.fMinTOF),
 123   fMaxTOF(copy.fMaxTOF),
 124   fLastRun(-1)
 125 {
 126 //
 127 // Copy constructor.
 128 //
 129
 130   Int_t i = 0;
 131   for (i = 0; i < 5; i++) fTPCpar[i] = copy.fTPCpar[i];
 132 }
 133
 134 //_________________________________________________________________________________________________
 135 AliRsnCutESD2010& AliRsnCutESD2010::operator=(const AliRsnCutESD2010& copy)
 136 {
 137 //
 138 // Assignment operator
 139 //
 140
 141   AliRsnCut::operator=(copy);
 142
 143   fIsMC = copy.fIsMC;
 144   fCheckITS = copy.fCheckITS;
 145   fCheckTPC = copy.fCheckTPC;
 146   fCheckTOF = copy.fCheckTOF;
 147   fUseGlobal = copy.fUseGlobal;
 148   fUseITSSA = copy.fUseITSSA;
 149   fPID = copy.fPID;
 150   fMaxEta = copy.fMaxEta;
 151   fMaxITSband = copy.fMaxITSband;
 152   fTPCpLimit = copy.fTPCpLimit;
 153   fMinTPCband = copy.fMinTPCband;
 154   fMaxTPCband = copy.fMaxTPCband;
 155   fESDtrackCutsTPC = copy.fESDtrackCutsTPC;
 156   fESDtrackCutsITS = copy.fESDtrackCutsITS;
 157   fTOFcalibrateESD = copy.fTOFcalibrateESD;
 158   fTOFcorrectTExp = copy.fTOFcorrectTExp;
 159   fTOFuseT0 = copy.fTOFuseT0;
 160   fTOFtuneMC = copy.fTOFtuneMC;
 161   fTOFresolution = copy.fTOFresolution;
 162   fMinTOF = copy.fMinTOF;
 163   fMaxTOF = copy.fMaxTOF;
 164   fLastRun = copy.fLastRun;
 165
 166   Int_t i = 0;
 167   for (i = 0; i < 5; i++) fTPCpar[i] = copy.fTPCpar[i];
 168
 169   return (*this);
 170 }
 171
 172 //_________________________________________________________________________________________________
 173 void AliRsnCutESD2010::SetMC(Bool_t isMC)
 174 {
 175 //
 176 // Sets some aspects of cuts depending on the fact that runs on MC or not
 177 //
 178
 179   fIsMC = isMC;
 180
 181   if (isMC)
 182   {
 183     SetTPCpar(2.15898 / 50.0, 1.75295E1, 3.40030E-9, 1.96178, 3.91720);
 184     SetTOFcalibrateESD(kFALSE);
 185     SetTOFtuneMC(kTRUE);
 186   }
 187   else
 188   {
 189     SetTPCpar(1.41543 / 50.0, 2.63394E1, 5.0411E-11, 2.12543, 4.88663);
 190     SetTOFcalibrateESD(kTRUE);
 191     SetTOFtuneMC(kFALSE);
 192   }
 193 }
 194
 195 //_________________________________________________________________________________________________
 196 void AliRsnCutESD2010::SetEvent(AliRsnEvent *event)
 197 {
 198   // don't do anything if the event is the same as before
 199   if (fEvent != 0x0 && fEvent == event) return;
 200
 201   // retrieve the ESD event
 202   AliESDEvent *esd = event->GetRefESD();
 203   if (!esd)
 204   {
 205     fEvent = 0x0;
 206     return;
 207   }
 208   else
 209   {
 210     fEvent = event;
 211   }
 212
 213   // if absent, initialize ESD pid responst
 214   if (!fESDpid)
 215   {
 216     fESDpid = new AliESDpid;
 217     fESDpid->GetTPCResponse().SetBetheBlochParameters(fTPCpar[0],fTPCpar[1],fTPCpar[2],fTPCpar[3],fTPCpar[4]);
 218   }
 219
 220   // initialize DB to current run
 221   Int_t run = esd->GetRunNumber();
 222   if (run != fLastRun)
 223   {
 224     cout << "Run = " << run << " -- LAST = " << fLastRun << endl;
 225     fLastRun = run;
 226
 227     // setup TOF maker & calibration
 228     if (!fTOFcalib) fTOFcalib = new AliTOFcalib;
 229     fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
 230     if (!fTOFmaker) fTOFmaker = new AliTOFT0maker(fESDpid, fTOFcalib);
 231     fTOFmaker->SetTimeResolution(fTOFresolution);
 232
 233     AliCDBManager *cdb = AliCDBManager::Instance();
 234     cdb->ClearCache(); // suggestion by Annalisa
 235     cdb->Clear();      // suggestion by Annalisa
 236     cdb->SetDefaultStorage("raw://");
 237     cdb->SetRun(run);
 238
 239     fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
 240     fTOFcalib->Init();
 241   }
 242
 243   // if required, calibrate the TOF t0 maker with current event
 244   if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(esd);
 245   if (fTOFtuneMC) fTOFmaker->TuneForMC(esd);
 246   if (fTOFuseT0)
 247   {
 248     fTOFmaker->ComputeT0TOF(esd);
 249     fTOFmaker->ApplyT0TOF(esd);
 250     fESDpid->MakePID(esd, kFALSE, 0.);
 251   }
 252 }
 253
 254 //_________________________________________________________________________________________________
 255 Bool_t AliRsnCutESD2010::IsSelected(TObject *obj1, TObject* /*obj2*/)
 256 {
 257 //
 258 // Cut checker.
 259 //
 260
 261   // coherence check: require an ESD track
 262   AliRsnDaughter *daughter = dynamic_cast<AliRsnDaughter*>(obj1);
 263   if (!daughter) return kFALSE;
 264   AliESDtrack *track = daughter->GetRefESDtrack();
 265   if (!track) return kFALSE;
 266
 267   // if no reference event, skip
 268   if (!fEvent) return kFALSE;
 269
 270   // ITS: create the response function
 271   AliITSPIDResponse itsrsp(fIsMC);
 272
 273   // TOF: define fixed function for compatibility range
 274   //Double_t a1 = 0.01, a2 = -0.03;
 275   //Double_t b1 = 0.25, b2 =  0.25;
 276   //Double_t c1 = 0.05, c2 = -0.03;
 277   //Double_t ymax, ymin;
 278
 279   ULong_t  status;
 280   Int_t    k, nITS;
 281   Double_t times[10], tpcNSigma, tpcMaxNSigma, itsSignal, itsNSigma, mom, tofTime, tofSigma, tofRef, tofRel;
 282   Bool_t   okQuality, okTOF, okTPC, okITS, isTPC, isITSSA, isTOF;
 283   UChar_t  itsCluMap;
 284
 285   // get commonly used variables
 286   status  = (ULong_t)track->GetStatus();
 287   mom     = track->P();
 288   isTPC   = ((status & AliESDtrack::kTPCin) != 0);
 289   isITSSA = ((status & AliESDtrack::kTPCin)  == 0 && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
 290   isTOF   = (((status & AliESDtrack::kTOFout) != 0) && ((status & AliESDtrack::kTIME) != 0) /* && mom > TMath::Max(b1, b2)*/);
 291
 292   // check if the track type matches what is required
 293   if (!isTPC && !isITSSA)
 294   {
 295     AliDebug(AliLog::kDebug + 2, "Track is not either a TPC track or a ITS standalone. Rejected");
 296     return kFALSE;
 297   }
 298   else if (isTPC && !fUseGlobal)
 299   {
 300     AliDebug(AliLog::kDebug + 2, "Global tracks not used. Rejected");
 301     return kFALSE;
 302   }
 303   else if (isITSSA && !fUseITSSA)
 304   {
 305     AliDebug(AliLog::kDebug + 2, "ITS standalone not used. Rejected");
 306     return kFALSE;
 307   }
 308
 309   // does a preliminary check on TOF values, if necessary
 310   // then, if the reference time or TOF signal are meaningless
 311   // even if the 'isTOF' flag is true, switch it to false
 312   if (isTOF)
 313   {
 314     track->GetIntegratedTimes(times);
 315     tofTime  = (Double_t)track->GetTOFsignal();
 316     tofSigma = fTOFmaker->GetExpectedSigma(mom, times[(Int_t)fPID], AliPID::ParticleMass((Int_t)fPID));
 317     tofRef   = times[(Int_t)fPID];
 318     if (tofRef <= 0.0 && tofSigma <= 0.0) isTOF = kFALSE;
 319   }
 320
 321   // check quality (eta range must be adapted)
 322   AliESDtrackCuts *cuts = 0x0;
 323   if (isTPC)   cuts = &fESDtrackCutsTPC;
 324   if (isITSSA) cuts = &fESDtrackCutsITS;
 325   if (!cuts)   return kFALSE;
 326   cuts->SetEtaRange(-fMaxEta, fMaxEta);
 327   okQuality = cuts->IsSelected(track);
 328   AliDebug(AliLog::kDebug + 2, Form("Global quality cut = %s", (okQuality ? "GOOD" : "BAD")));
 329   if (!okQuality) return kFALSE;
 330
 331   if (isTPC) // this branch is entered by all global tracks
 332   {
 333     // check TPC dE/dx:
 334     if (fCheckTPC)
 335     {
 336       tpcNSigma = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, fPID));
 337       if (track->GetInnerParam()->P() > fTPCpLimit) tpcMaxNSigma = fMinTPCband; else tpcMaxNSigma = fMaxTPCband;
 338       okTPC = (tpcNSigma <= tpcMaxNSigma);
 339       AliDebug(AliLog::kDebug + 2, Form("TPC nsigma = %f -- max = %f -- cut %s", tpcNSigma, tpcMaxNSigma, (okTPC ? "passed" : "failed")));
 340     }
 341     else
 342     {
 343       // if TPC is not checked, it is as if all tracks do pass the cut
 344       okTPC = kTRUE;
 345       AliDebug(AliLog::kDebug + 2, "TPC not checked, track accepted");
 346     }
 347
 348     // check TOF (only if flags are OK)
 349     if (fCheckTOF)
 350     {
 351       if (isTOF)
 352       {
 353         // TOF can be checked only when track is matched there
 354         track->GetIntegratedTimes(times);
 355         tofTime  = (Double_t)track->GetTOFsignal();
 356         tofSigma = fTOFmaker->GetExpectedSigma(mom, times[(Int_t)fPID], AliPID::ParticleMass((Int_t)fPID));
 357         tofRef   = times[(Int_t)fPID];
 358         /*
 359         tofRel   = (tofTime - tofRef) / tofRef;
 360         ymax     = a1 / (mom - b1) + c1;
 361         ymin     = a2 / (mom - b2) + c2;
 362         okTOF    = (tofRel >= ymin && tofRel <= ymax);
 363         */
 364         tofRel   = (tofTime - tofRef) / tofSigma;
 365         okTOF    = (tofRel >= fMinTOF && tofRel <= fMaxTOF);
 366         AliDebug(AliLog::kDebug + 2, Form("TOF nsigma = %f -- range = %f %f -- cut %s", tofRel, fMinTOF, fMaxTOF, (okTOF ? "passed" : "failed")));
 367       }
 368       else
 369       {
 370         // if TOF is not matched, the answer depends on TPC:
 371         // - if TPC is required, track is checked only there and TOF simply ignored
 372         // - if TPC is not required, track is rejected when TOF does not match it, if TOF check is required
 373         if (fCheckTPC) okTOF = kTRUE; else okTOF = kFALSE;
 374       }
 375     }
 376     else
 377     {
 378       okTOF = kTRUE;
 379     }
 380
 381     // properly combine the outcome of TPC and TOF cuts
 382     return okTPC && okTOF;
 383   }
 384   else if (isITSSA) // this branch is entered by all ITS standalone tracks
 385   {
 386     // check dE/dx only if this is required, otherwise ITS standalone are just added but not checked for PID
 387     if (fCheckITS)
 388     {
 389       itsSignal = track->GetITSsignal();
 390       itsCluMap = track->GetITSClusterMap();
 391       nITS      = 0;
 392       for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
 393       if (nITS < 3) return kFALSE;
 394       itsNSigma = itsrsp.GetNumberOfSigmas(mom, itsSignal, fPID, nITS, kTRUE);
 395       okITS = (TMath::Abs(itsNSigma) <= fMaxITSband);
 396       AliDebug(AliLog::kDebug + 2, Form("ITS nsigma = %f -- max = %f -- cut %s", itsNSigma, fMaxITSband, (okITS ? "passed" : "failed")));
 397     }
 398     else
 399     {
 400       okITS = kTRUE;
 401     }
 402
 403     return okITS;
 404   }
 405   else
 406   {
 407     // if we are here, the track is surely bad
 408     return kFALSE;
 409   }
 410 }