remove option C for Clear for trigger array for the moment, causes malloc error
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / 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 }