PWG2/RESONANCES/extra/AliRsnAnalysisPhi7TeV.cxx

   1 //
   2 // Implementation file for implementation of data analysis aft 900 GeV
   3 //
   4 // Author: A. Pulvirenti
   5 //
   6
   7 #include "Riostream.h"
   8 #include <iomanip>
   9
  10 #include "TH1.h"
  11 #include "TH3.h"
  12 #include "TFile.h"
  13 #include "TTree.h"
  14 #include "TParticle.h"
  15 #include "TRandom.h"
  16 #include "TLorentzVector.h"
  17
  18 #include "AliLog.h"
  19 #include "AliESDpid.h"
  20 #include "AliESDEvent.h"
  21 #include "AliESDVertex.h"
  22 #include "AliESDtrack.h"
  23 #include "AliESDtrackCuts.h"
  24 #include "AliStack.h"
  25 #include "AliMCEvent.h"
  26 #include "AliTOFT0maker.h"
  27 #include "AliTOFcalib.h"
  28 #include "AliCDBManager.h"
  29 #include "AliITSPIDResponse.h"
  30
  31 #include "AliRsnEvent.h"
  32 #include "AliRsnTarget.h"
  33 #include "AliRsnDaughter.h"
  34 #include "AliRsnCutESD2010.h"
  35
  36 #include "AliRsnAnalysisPhi7TeV.h"
  37
  38 //__________________________________________________________________________________________________
  39 AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const char *name, Bool_t isMC) :
  40   AliAnalysisTaskSE(name),
  41   fUseMC(kFALSE),
  42   fCheckITS(kTRUE),
  43   fCheckTPC(kTRUE),
  44   fCheckTOF(kTRUE),
  45   fAddITSSA(kFALSE),
  46   fMaxVz(1E6),
  47   fMaxITSband(1E6),
  48   fMaxITSmom(0.0),
  49   fTPCpLimit(0.35),
  50   fMinTPCband(-1E6),
  51   fMaxTPCband( 1E6),
  52   fMinTOF(-3.0),
  53   fMaxTOF( 3.0),
  54   fOutList(0x0),
  55   fUnlike(0x0),
  56   fLikePP(0x0),
  57   fLikeMM(0x0),
  58   fTrues(0x0),
  59   fHEvents(0x0),
  60   fESDtrackCutsTPC(),
  61   fESDtrackCutsITS(),
  62   fESDpid(0x0),
  63   fTOFmaker(0x0),
  64   fTOFcalib(0x0),
  65   fTOFcalibrateESD(!isMC),
  66   fTOFcorrectTExp(kTRUE),
  67   fTOFuseT0(kTRUE),
  68   fTOFtuneMC(isMC),
  69   fTOFresolution(100.0),
  70   fDaughter(),
  71   fRsnCuts()
  72 {
  73 //
  74 // Constructor
  75 //
  76
  77   SetUseMC(isMC);
  78   DefineOutput(1, TList::Class());
  79 }
  80
  81 //__________________________________________________________________________________________________
  82 AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const AliRsnAnalysisPhi7TeV& copy) :
  83   AliAnalysisTaskSE(copy),
  84   fUseMC(copy.fUseMC),
  85   fCheckITS(copy.fCheckITS),
  86   fCheckTPC(copy.fCheckTPC),
  87   fCheckTOF(copy.fCheckTOF),
  88   fAddITSSA(copy.fAddITSSA),
  89   fMaxVz(copy.fMaxVz),
  90   fMaxITSband(copy.fMaxITSband),
  91   fMaxITSmom(copy.fMaxITSmom),
  92   fTPCpLimit(copy.fTPCpLimit),
  93   fMinTPCband(copy.fMinTPCband),
  94   fMaxTPCband(copy.fMaxTPCband),
  95   fMinTOF(copy.fMinTOF),
  96   fMaxTOF(copy.fMaxTOF),
  97   fOutList(0x0),
  98   fUnlike(0x0),
  99   fLikePP(0x0),
 100   fLikeMM(0x0),
 101   fTrues(0x0),
 102   fHEvents(0x0),
 103   fESDtrackCutsTPC(copy.fESDtrackCutsTPC),
 104   fESDtrackCutsITS(copy.fESDtrackCutsITS),
 105   fESDpid(0x0),
 106   fTOFmaker(0x0),
 107   fTOFcalib(0x0),
 108   fTOFcalibrateESD(kFALSE),
 109   fTOFcorrectTExp(kFALSE),
 110   fTOFuseT0(kFALSE),
 111   fTOFtuneMC(kFALSE),
 112   fTOFresolution(0.0),
 113   fDaughter(),
 114   fRsnCuts()
 115 {
 116 //
 117 // Copy constructor
 118 //
 119
 120   SetUseMC(copy.fUseMC);
 121 }
 122
 123 //__________________________________________________________________________________________________
 124 AliRsnAnalysisPhi7TeV& AliRsnAnalysisPhi7TeV::operator=(const AliRsnAnalysisPhi7TeV& copy)
 125 {
 126 //
 127 // Assignment operator
 128 //
 129
 130   fUseMC = copy.fUseMC;
 131   fCheckITS = copy.fCheckITS;
 132   fCheckTPC = copy.fCheckTPC;
 133   fCheckTOF = copy.fCheckTOF;
 134   fAddITSSA = copy.fAddITSSA;
 135
 136   fMaxVz   = copy.fMaxVz;
 137   fMaxITSband = copy.fMaxITSband;
 138   fMaxITSmom  = copy.fMaxITSmom;
 139
 140   fTPCpLimit  = copy.fTPCpLimit;
 141   fMinTPCband = copy.fMinTPCband;
 142   fMaxTPCband = copy.fMaxTPCband;
 143
 144   fESDtrackCutsTPC = copy.fESDtrackCutsTPC;
 145   fESDtrackCutsITS = copy.fESDtrackCutsITS;
 146
 147   fTOFcalibrateESD = copy.fTOFcalibrateESD;
 148   fTOFcorrectTExp = copy.fTOFcorrectTExp;
 149   fTOFuseT0 = copy.fTOFuseT0;
 150   fTOFtuneMC = copy.fTOFtuneMC;
 151   fTOFresolution = copy.fTOFresolution;
 152
 153   SetUseMC(copy.fUseMC);
 154
 155   return (*this);
 156 }
 157
 158 //__________________________________________________________________________________________________
 159 AliRsnAnalysisPhi7TeV::~AliRsnAnalysisPhi7TeV()
 160 {
 161 //
 162 // Destructor
 163 //
 164 }
 165
 166 //_________________________________________________________________________________________________
 167 void AliRsnAnalysisPhi7TeV::SetUseMC(Bool_t isMC)
 168 {
 169 //
 170 // Sets some aspects of cuts depending on the fact that runs on MC or not
 171 //
 172
 173   fUseMC = isMC;
 174
 175   if (isMC)
 176   {
 177     SetTPCpar(2.15898 / 50.0, 1.75295E1, 3.40030E-9, 1.96178, 3.91720);
 178     SetTOFcalibrateESD(kFALSE);
 179     SetTOFtuneMC(kTRUE);
 180   }
 181   else
 182   {
 183     SetTPCpar(1.41543 / 50.0, 2.63394E1, 5.0411E-11, 2.12543, 4.88663);
 184     SetTOFcalibrateESD(kTRUE);
 185     SetTOFtuneMC(kFALSE);
 186   }
 187 }
 188
 189 //__________________________________________________________________________________________________
 190 void AliRsnAnalysisPhi7TeV::UserCreateOutputObjects()
 191 {
 192 //
 193 // Create the output data container
 194 //
 195
 196   // initialize random
 197   //gRandom->SetSeed(0);
 198
 199   // create output trees
 200   OpenFile(1);
 201   fOutList    = new TList;
 202   fHEvents    = new TH1I("hEvents", "Event details", 6, 0, 6);
 203   fVertexX[0] = new TH1F("hVertexTracksX", "X position of primary vertex (tracks)", 200,  -2,  2);
 204   fVertexY[0] = new TH1F("hVertexTracksY", "Y position of primary vertex (tracks)", 200,  -2,  2);
 205   fVertexZ[0] = new TH1F("hVertexTracksZ", "Z position of primary vertex (tracks)", 400, -40, 40);
 206   fVertexX[1] = new TH1F("hVertexSPDX", "X position of primary vertex (SPD)", 200,  -2,  2);
 207   fVertexY[1] = new TH1F("hVertexSPDY", "Y position of primary vertex (SPD)", 200,  -2,  2);
 208   fVertexZ[1] = new TH1F("hVertexSPDZ", "Z position of primary vertex (SPD)", 400, -40, 40);
 209   //fUnlike     = new TH3F("hPM", "+-   pairs", 500, 0.9, 1.4, 100, 0.0, 10.0, 24, -1.2, 1.2);
 210   //fLikePP     = new TH3F("hPP", "++   pairs", 500, 0.9, 1.4, 100, 0.0, 10.0, 24, -1.2, 1.2);
 211   //fLikeMM     = new TH3F("hMM", "--   pairs", 500, 0.9, 1.4, 100, 0.0, 10.0, 24, -1.2, 1.2);
 212   //fTrues      = new TH3F("hTR", "True pairs", 500, 0.9, 1.4, 100, 0.0, 10.0, 24, -1.2, 1.2);
 213   TFile *ffile = TFile::Open("template.root");
 214   if (ffile)
 215   {
 216     TH3F *tmp = (TH3F*)ffile->Get("template");
 217     fUnlike = (TH3F*)tmp->Clone("hPM");
 218     fLikePP = (TH3F*)tmp->Clone("hPP");
 219     fLikeMM = (TH3F*)tmp->Clone("hMM");
 220     fTrues  = (TH3F*)tmp->Clone("hTR");
 221   }
 222
 223   fHEvents->GetXaxis()->SetBinLabel(1, "Good vertex with tracks");
 224   fHEvents->GetXaxis()->SetBinLabel(2, "Good vertex with SPD");
 225   fHEvents->GetXaxis()->SetBinLabel(3, "Far vertex with tracks");
 226   fHEvents->GetXaxis()->SetBinLabel(4, "Far vertex with SPD");
 227   fHEvents->GetXaxis()->SetBinLabel(5, "No good vertex");
 228   fHEvents->GetXaxis()->SetBinLabel(6, "Empty event");
 229
 230   fOutList->Add(fHEvents);
 231   fOutList->Add(fVertexX[0]);
 232   fOutList->Add(fVertexY[0]);
 233   fOutList->Add(fVertexZ[0]);
 234   fOutList->Add(fVertexX[1]);
 235   fOutList->Add(fVertexY[1]);
 236   fOutList->Add(fVertexZ[1]);
 237   fOutList->Add(fUnlike);
 238   fOutList->Add(fLikePP);
 239   fOutList->Add(fLikeMM);
 240   fOutList->Add(fTrues);
 241
 242   // setup RSN-related objects
 243
 244   fRsnCuts.SetMC       (fUseMC);
 245   fRsnCuts.SetCheckITS (fCheckITS);
 246   fRsnCuts.SetCheckTPC (fCheckTPC);
 247   fRsnCuts.SetCheckTOF (fCheckTOF);
 248   fRsnCuts.SetUseITSTPC(kTRUE);
 249   fRsnCuts.SetUseITSSA (fAddITSSA);
 250   fRsnCuts.SetPID      (AliPID::kKaon);
 251   fRsnCuts.SetMaxITSPIDmom(fMaxITSmom);
 252   fRsnCuts.SetITSband(fMaxITSband);
 253   fRsnCuts.SetTPCpLimit(fTPCpLimit);
 254   fRsnCuts.SetTPCrange(fMinTPCband, fMaxTPCband);
 255   fRsnCuts.SetTPCpar(fTPCpar[0],fTPCpar[1],fTPCpar[2],fTPCpar[3],fTPCpar[4]);
 256   //fRsnCuts.SetTOFcalibrateESD(fTOFcalibrateESD);
 257   //fRsnCuts.SetTOFcorrectTExp (fTOFcorrectTExp);
 258   //fRsnCuts.SetTOFuseT0       (fTOFuseT0);
 259   //fRsnCuts.SetTOFtuneMC      (fTOFtuneMC);
 260   //fRsnCuts.SetTOFresolution  (fTOFresolution);
 261   fRsnCuts.SetTOFrange       (fMinTOF, fMaxTOF);
 262   fRsnCuts.GetCutsTPC()->Copy(fESDtrackCutsTPC);
 263   fRsnCuts.GetCutsITS()->Copy(fESDtrackCutsITS);
 264 }
 265
 266 //__________________________________________________________________________________________________
 267 void AliRsnAnalysisPhi7TeV::UserExec(Option_t *)
 268 {
 269 //
 270 // Main execution function.
 271 // Fills the fHEvents data member with the following legenda:
 272 // 0 -- event OK, prim vertex with tracks
 273 // 1 -- event OK, prim vertex with SPD
 274 // 2 -- event OK but vz large
 275 // 3 -- event bad
 276 //
 277
 278   // retrieve ESD event and related stack (if available)
 279   AliESDEvent *esd   = dynamic_cast<AliESDEvent*>(fInputEvent);
 280   AliStack    *stack = (fMCEvent ? fMCEvent->Stack() : 0x0);
 281
 282   // check the event
 283   Int_t eval = EventEval(esd);
 284   fHEvents->Fill(eval);
 285
 286   // if the event is good for analysis, process it
 287   if (eval == kGoodTracksPrimaryVertex || eval == kGoodSPDPrimaryVertex)
 288   {
 289     ProcessESD(esd, stack);
 290   }
 291
 292   // update histogram container
 293   PostData(1, fOutList);
 294 }
 295
 296 //__________________________________________________________________________________________________
 297 void AliRsnAnalysisPhi7TeV::Terminate(Option_t *)
 298 {
 299 //
 300 // Terminate
 301 //
 302 }
 303
 304 //__________________________________________________________________________________________________
 305 Int_t AliRsnAnalysisPhi7TeV::EventEval(AliESDEvent *esd)
 306 {
 307 //
 308 // Checks if the event is good for analysis.
 309 // Returns one of the flag values defined in the header
 310 //
 311
 312   static Int_t evNum = 0;
 313   evNum++;
 314
 315   // reject empty events
 316   Int_t ntracks = esd->GetNumberOfTracks();
 317   if (!ntracks) return kEmptyEvent;
 318
 319   // get the best primary vertex:
 320   // first try the one with tracks
 321   const AliESDVertex *vTrk  = esd->GetPrimaryVertexTracks();
 322   const AliESDVertex *vSPD  = esd->GetPrimaryVertexSPD();
 323   Double_t            vzTrk = 1000000.0;
 324   Double_t            vzSPD = 1000000.0;
 325   Int_t               ncTrk = -1;
 326   Int_t               ncSPD = -1;
 327   if (vTrk) ncTrk = (Int_t)vTrk->GetNContributors();
 328   if (vSPD) ncSPD = (Int_t)vSPD->GetNContributors();
 329   if (vTrk) vzTrk = TMath::Abs(vTrk->GetZv());
 330   if (vSPD) vzSPD = TMath::Abs(vSPD->GetZv());
 331   if(vTrk && ncTrk > 0)
 332   {
 333     // fill the histograms
 334     fVertexX[0]->Fill(vTrk->GetXv());
 335     fVertexY[0]->Fill(vTrk->GetYv());
 336     fVertexZ[0]->Fill(vTrk->GetZv());
 337
 338     // check VZ position
 339     if (vzTrk <= fMaxVz)
 340       return kGoodTracksPrimaryVertex;
 341     else
 342       return kFarTracksPrimaryVertex;
 343   }
 344   else if (vSPD && ncSPD > 0)
 345   {
 346     // fill the histograms
 347     fVertexX[1]->Fill(vSPD->GetXv());
 348     fVertexY[1]->Fill(vSPD->GetYv());
 349     fVertexZ[1]->Fill(vSPD->GetZv());
 350
 351     // check VZ position
 352     if (vzSPD <= fMaxVz)
 353       return kGoodSPDPrimaryVertex;
 354     else
 355       return kFarSPDPrimaryVertex;
 356   }
 357   else
 358     return kNoGoodPrimaryVertex;
 359 }
 360
 361 //__________________________________________________________________________________________________
 362 void AliRsnAnalysisPhi7TeV::ProcessESD
 363 (AliESDEvent *esd, AliStack *stack)
 364 {
 365 //
 366 // This function works with the ESD object
 367 //
 368
 369   static Int_t lastRun = -1;
 370   static Int_t evnum = 0;
 371   evnum++;
 372
 373   // get current run
 374   Int_t run = esd->GetRunNumber();
 375
 376   // if absent, initialize ESD pid response
 377   if (!fESDpid)
 378   {
 379     AliITSPIDResponse itsresponse(fUseMC);
 380
 381     fESDpid = new AliESDpid;
 382     fESDpid->GetTPCResponse().SetBetheBlochParameters(fTPCpar[0],fTPCpar[1],fTPCpar[2],fTPCpar[3],fTPCpar[4]);
 383     fESDpid->GetITSResponse() = itsresponse;
 384   }
 385
 386   // if the run number has changed, update it now and give a message
 387   if (run != lastRun)
 388   {
 389     AliInfo("============================================================================================");
 390     AliInfo(Form("*** CHANGING RUN NUMBER: PREVIOUS = %d --> CURRENT = %d ***", lastRun, run));
 391     AliInfo("============================================================================================");
 392     lastRun = run;
 393
 394     AliCDBManager::Instance()->SetDefaultStorage("raw://");
 395     AliCDBManager::Instance()->SetRun(lastRun);
 396
 397     if (fTOFmaker) delete fTOFmaker;
 398     if (fTOFcalib) delete fTOFcalib;
 399
 400     fTOFcalib = new AliTOFcalib();
 401     if (fTOFcorrectTExp) fTOFcalib->SetCorrectTExp(kTRUE);
 402     fTOFcalib->Init();
 403
 404     fTOFmaker = new AliTOFT0maker(fESDpid, fTOFcalib);
 405     fTOFmaker->SetTimeResolution(fTOFresolution);
 406   }
 407
 408   if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(esd);
 409   if (fTOFtuneMC) fTOFmaker->TuneForMC(esd);
 410   if (fTOFuseT0)
 411   {
 412     fTOFmaker->ComputeT0TOF(esd);
 413     fTOFmaker->ApplyT0TOF(esd);
 414     fESDpid->MakePID(esd, kFALSE, 0.);
 415   }
 416
 417   // RSN event
 418   AliRsnEvent event;
 419   event.SetRef(esd);
 420   event.SetRefMC(fMCEvent);
 421   //fRsnCuts.ProcessEvent(esd);
 422
 423   // loop on all tracks
 424   Int_t           i1, i2, label, pdg, ntracks = esd->GetNumberOfTracks();
 425   Bool_t          okTrack;
 426   AliPID          pid;
 427   Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
 428   Double_t        phimass = 1.019455;
 429   Double_t        angle, cosangle;
 430   AliMCParticle  *p1 = 0x0, *p2 = 0x0;
 431   AliESDtrack    *t1 = 0x0, *t2 = 0x0;
 432   TLorentzVector  v1, v2, vsum, vref;
 433
 434   // external loop (T1)
 435   for (i1 = 0; i1 < ntracks; i1++)
 436   {
 437     t1 = esd->GetTrack(i1);
 438     if (!t1) continue;
 439
 440     // setup RSN
 441     fDaughter.SetRef(t1);
 442     if (stack)
 443     {
 444       p1 = (AliMCParticle*)fMCEvent->GetTrack(TMath::Abs(t1->GetLabel()));
 445       if (p1) fDaughter.SetRefMC(p1);
 446     }
 447     fDaughter.SetMass(kmass);
 448     v1.SetXYZM(t1->Px(), t1->Py(), t1->Pz(), kmass);
 449
 450     // check track
 451     okTrack = OkTrack(t1, AliPID::kKaon);
 452
 453     // internal loop (T2)
 454     for (i2 = i1+1; i2 < ntracks; i2++)
 455     {
 456       t2 = esd->GetTrack(i2);
 457       if (!t2) continue;
 458
 459       // check track
 460       if (!OkTrack(t2, AliPID::kKaon)) continue;
 461
 462       // if unlike pair, check if it is true
 463       pdg = 0;
 464       if ((t1->Charge() == t2->Charge()) && p1 && p2)
 465       {
 466         TParticle *part1 = p1->Particle();
 467         TParticle *part2 = p2->Particle();
 468         if (TMath::Abs(part1->GetPdgCode()) == 321 && TMath::Abs(part2->GetPdgCode()) == 321 && part1->GetFirstMother() == part2->GetFirstMother())
 469         {
 470           label = part1->GetFirstMother();
 471           if (label >= 0 && label < stack->GetNtrack())
 472           {
 473             TParticle *mum = stack->Particle(label);
 474             pdg = mum->GetPdgCode();
 475           }
 476         }
 477       }
 478       pdg = TMath::Abs(pdg);
 479
 480       // combine momenta
 481       v1.SetXYZM(t1->Px(), t1->Py(), t1->Pz(), kmass);
 482       v2.SetXYZM(t2->Px(), t2->Py(), t2->Pz(), kmass);
 483       angle = v1.Angle(v2.Vect());
 484       cosangle = TMath::Abs(TMath::Cos(angle));
 485       if (cosangle <= 0.02) continue;
 486       vsum = v1 + v2;
 487       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
 488
 489       // fill appropriate histogram
 490       if (t1->Charge() != t2->Charge())
 491       {
 492         fUnlike->Fill(vsum.M(), vsum.Perp(), vref.Rapidity());
 493         if (pdg == 333) fTrues->Fill(vsum.M(), vsum.Perp(), vref.Rapidity());
 494       }
 495       else if (t1->Charge() > 0)
 496       {
 497         fLikePP->Fill(vsum.M(), vsum.Perp(), vref.Rapidity());
 498       }
 499       else
 500       {
 501         fLikeMM->Fill(vsum.M(), vsum.Perp(), vref.Rapidity());
 502       }
 503     }
 504   }
 505
 506   PostData(1, fOutList);
 507 }
 508
 509 //______________________________________________________________________________
 510 Bool_t AliRsnAnalysisPhi7TeV::OkQuality(AliESDtrack *track)
 511 {
 512 //
 513 // Check track quality parameters.
 514 // Rejects all tracks which are not either TPC+ITS nor ITS standalone.
 515 // If tracks of any type are not flagged to be used, they are rejected anyway.
 516 //
 517
 518   if (IsITSTPC(track))
 519     return fESDtrackCutsTPC.IsSelected(track);
 520   else if (IsITSSA (track))
 521   {
 522     if (fAddITSSA)
 523       return fESDtrackCutsITS.IsSelected(track);
 524     else
 525       return kFALSE;
 526   }
 527   else
 528     return kFALSE;
 529 }
 530
 531 //______________________________________________________________________________
 532 Bool_t AliRsnAnalysisPhi7TeV::OkITSPID (AliESDtrack *track, AliPID::EParticleType pid)
 533 {
 534 //
 535 // Check ITS particle identification with 3sigma cut
 536 //
 537
 538   // reject not ITS standalone tracks
 539   if (!IsITSSA(track)) return kFALSE;
 540
 541   // count PID layers and reject if they are too few
 542   Int_t   k, nITSpidLayers = 0;
 543   UChar_t itsCluMap = track->GetITSClusterMap();
 544   for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITSpidLayers;
 545   if (nITSpidLayers < 3)
 546   {
 547     AliDebug(AliLog::kDebug+2, "Rejecting track with too few ITS pid layers");
 548     return kFALSE;
 549   }
 550
 551   // check the track type (ITS+TPC or ITS standalone)
 552   // and reject it if it is of none of the allowed types
 553   Bool_t isSA = kFALSE;
 554   if (IsITSTPC(track)) isSA = kFALSE;
 555   else if (IsITSSA(track)) isSA = kTRUE;
 556   else
 557   {
 558     AliWarning("Track is neither ITS+TPC nor ITS standalone");
 559     return kFALSE;
 560   }
 561
 562   // create the PID response object and compute nsigma
 563   AliITSPIDResponse &itsrsp = fESDpid->GetITSResponse();
 564   Double_t mom    = track->P();
 565   Double_t nSigma = itsrsp.GetNumberOfSigmas(mom, track->GetITSsignal(), pid, nITSpidLayers, isSA);
 566
 567   // evaluate the cut
 568   Bool_t ok = (TMath::Abs(nSigma) <= fMaxITSband);
 569
 570   // debug message
 571   AliDebug(AliLog::kDebug + 2, Form("ITS nsigma = %f -- max = %f -- cut %s", nSigma, fMaxITSband, (ok ? "passed" : "failed")));
 572
 573   // outcome
 574   return ok;
 575 }
 576
 577 //______________________________________________________________________________
 578 Bool_t AliRsnAnalysisPhi7TeV::OkTPCPID (AliESDtrack *track, AliPID::EParticleType pid)
 579 {
 580 //
 581 // Check TPC particle identification with {3|5}sigmacut,
 582 // depending on the track total momentum.
 583 //
 584
 585   // reject not TPC tracks
 586   if (!IsITSTPC(track)) return kFALSE;
 587
 588   // setup TPC PID response
 589   AliTPCPIDResponse &tpcrsp = fESDpid->GetTPCResponse();
 590   tpcrsp.SetBetheBlochParameters(fTPCpar[0],fTPCpar[1],fTPCpar[2],fTPCpar[3],fTPCpar[4]);
 591
 592   // get momentum and number of sigmas and choose the reference band
 593   Double_t mom       = track->GetInnerParam()->P();
 594   Double_t nSigma    = tpcrsp.GetNumberOfSigmas(mom, track->GetTPCsignal(), track->GetTPCsignalN(), pid);
 595   Double_t maxNSigma = fMaxTPCband;
 596   if (mom < fTPCpLimit) maxNSigma = fMinTPCband;
 597
 598   // evaluate the cut
 599   Bool_t ok = (TMath::Abs(nSigma) <= maxNSigma);
 600
 601   // debug message
 602   AliDebug(AliLog::kDebug + 2, Form("TPC nsigma = %f -- max = %f -- cut %s", nSigma, maxNSigma, (ok ? "passed" : "failed")));
 603
 604   // outcome
 605   return ok;
 606 }
 607
 608 //______________________________________________________________________________
 609 Bool_t AliRsnAnalysisPhi7TeV::OkTOFPID (AliESDtrack *track, AliPID::EParticleType pid)
 610 {
 611 //
 612 // Check TOF particle identification if matched there.
 613 //
 614
 615   // reject not TOF-matched tracks
 616   if (!MatchTOF(track)) return kFALSE;
 617
 618   // setup TOF PID response
 619   AliTOFPIDResponse &tofrsp = fESDpid->GetTOFResponse();
 620
 621   // get info for computation
 622   Double_t momentum = track->P();
 623   Double_t time     = track->GetTOFsignal();
 624   Double_t timeint[AliPID::kSPECIES];
 625   tofrsp.GetStartTime(momentum);
 626   track->GetIntegratedTimes(timeint);
 627
 628   // check the cut
 629   Double_t timeDiff = time - timeint[(Int_t)pid];
 630   Double_t sigmaRef = tofrsp.GetExpectedSigma(momentum, timeint[(Int_t)pid], AliPID::ParticleMass(pid));
 631   Double_t nSigma   = timeDiff / sigmaRef;
 632
 633   // evaluate the cut
 634   Bool_t ok = (nSigma >= fMinTOF && nSigma <= fMaxTOF);
 635
 636   // debug message
 637   AliDebug(AliLog::kDebug + 2, Form("TOF nsigma = %f -- range = %f - %f -- cut %s", nSigma, fMinTOF, fMaxTOF, (ok ? "passed" : "failed")));
 638   //if (print) cout << Form("**PHI** TOF nsigma = %f -- range = %f - %f -- cut %s", nSigma, fMinTOF, fMaxTOF, (ok ? "passed" : "failed")) << endl;
 639
 640   // outcome
 641   return ok;
 642 }
 643
 644 //______________________________________________________________________________
 645 Bool_t AliRsnAnalysisPhi7TeV::OkTrack(AliESDtrack *track, AliPID::EParticleType pid)
 646 {
 647 //
 648 // Global track cut check
 649 //
 650
 651   Bool_t okITS, okTPC, okTOF;
 652
 653   // check quality and track type and reject tracks not passing this step
 654   if (!OkQuality(track))
 655   {
 656     AliDebug(AliLog::kDebug+2, "Failed quality cut");
 657     //printf("[PHI] Track %4d: REJECTED\n", i);
 658     return kFALSE;
 659   }
 660
 661   // ITS PID can be checked always
 662   // if PID is not required, the flag is sed as
 663   // if the cut was alsways passed
 664   okITS = OkITSPID(track, AliPID::kKaon);
 665   if (!fCheckITS) okITS = kTRUE;
 666
 667   // TPC PID can be checked only for TPC+ITS tracks
 668   // if PID is not required, the flag is sed as
 669   // if the cut was alsways passed
 670   okTPC = kFALSE;
 671   if (IsITSTPC(track)) okTPC = OkTPCPID(track, AliPID::kKaon);
 672   if (!fCheckTPC) okTPC = kTRUE;
 673
 674   // TOF PID can be checked only if TOF is matched
 675   // if PID is not required, the flag is sed as
 676   // if the cut was alsways passed
 677   okTOF = kFALSE;
 678   if (IsITSTPC(track) && MatchTOF(track)) okTOF = OkTOFPID(track, AliPID::kKaon);
 679   if (!fCheckTOF) okTOF = kTRUE;
 680
 681   // now combine all outcomes according to the different possibilities:
 682   // -- ITS standalone:
 683   //    --> only ITS PID, always
 684   // -- ITS + TPC:
 685   //    --> ITS PID, only for momenta lower than 'fMaxITSPIDmom' and when the ITSpid flag is active
 686   //    --> TPC PID, always --> MASTER (first to be checked, if fails, track is rejected)
 687   //    --> TOF PID, only if matched
 688   if (IsITSSA(track))
 689   {
 690     if (!okITS)
 691     {
 692       AliDebug(AliLog::kDebug+2, "ITS standalone track --> ITS PID failed");
 693       //printf("[PHI] Track %4d: REJECTED\n", i);
 694       return kFALSE;
 695     }
 696   }
 697   else // checking IsITSTPC() is redundant due to OkQuality() cut check
 698   {
 699     if (!okTPC)
 700     {
 701       AliDebug(AliLog::kDebug+2, "ITS+TPC track --> TPC PID failed");
 702       //printf("[PHI] Track %4d: REJECTED\n", i);
 703       return kFALSE;
 704     }
 705     else if (MatchTOF(track) && !okTOF)
 706     {
 707       AliDebug(AliLog::kDebug+2, "ITS+TPC track --> TOF matched but TOF PID failed");
 708       //printf("[PHI] Track %4d: REJECTED\n", i);
 709       return kFALSE;
 710     }
 711     else if (track->IsOn(AliESDtrack::kITSpid) && track->P() <= fMaxITSmom && !okITS)
 712     {
 713       AliDebug(AliLog::kDebug+2, Form("ITS+TPC track --> Momentum lower than limit (%.2f) and ITS PID failed", fMaxITSmom));
 714       //printf("[PHI] Track %4d: REJECTED\n", i);
 715       return kFALSE;
 716     }
 717   }
 718
 719   // this point is reached only if function didn't exit before due to somecheck not passed
 720   return kTRUE;
 721 }