Splitted the PID cuts into three classes for ITS, TPC and TOF, and added the possibil...
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / extra / AliRsnAnalysisPhi7TeVNoPID.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 "TTree.h"
12 #include "TParticle.h"
13 #include "TRandom.h"
14 #include "TLorentzVector.h"
15
16 #include "AliLog.h"
17 #include "AliESDpid.h"
18 #include "AliESDEvent.h"
19 #include "AliESDVertex.h"
20 #include "AliESDtrack.h"
21 #include "AliStack.h"
22 #include "AliMCEvent.h"
23 #include "AliTOFT0maker.h"
24 #include "AliTOFcalib.h"
25 #include "AliCDBManager.h"
26 #include "AliITSPIDResponse.h"
27
28 #include "AliRsnAnalysisPhi7TeVNoPID.h"
29
30 //__________________________________________________________________________________________________
31 AliRsnAnalysisPhi7TeVNoPID::AliRsnAnalysisPhi7TeVNoPID(const char *name) :
32   AliAnalysisTaskSE(name),
33   fUseMC(kFALSE),
34   fPDG(0),
35   fCh(0),
36   fIM(0.0),
37   fPt(0.0),
38   fY(0.0),
39   fEta(0.0),
40   fMaxVz(1E6),
41   fMaxITSband(1E6),
42   fTPCpLimit(0.35),
43   fMinTPCband(-1E6),
44   fMaxTPCband( 1E6),
45   fRsnTreeComp(0x0),
46   fRsnTreeTrue(0x0),
47   fOutList(0x0),
48   fHEvents(0x0),
49   fESDtrackCutsTPC(),
50   fESDtrackCutsITS(),
51   fESDpid(0x0),
52   fTOFmaker(0x0),
53   fTOFcalib(0x0),
54   fTOFcalibrateESD(kFALSE),
55   fTOFcorrectTExp(kFALSE),
56   fTOFuseT0(kFALSE),
57   fTOFtuneMC(kFALSE),
58   fTOFresolution(0.0)
59   
60 {
61 //
62 // Constructor
63 //
64
65   DefineOutput(1, TTree::Class());
66   DefineOutput(2, TTree::Class());
67   DefineOutput(3, TList::Class());
68 }
69
70 //__________________________________________________________________________________________________
71 AliRsnAnalysisPhi7TeVNoPID::AliRsnAnalysisPhi7TeVNoPID(const AliRsnAnalysisPhi7TeVNoPID& copy) :
72   AliAnalysisTaskSE(copy),
73   fUseMC(copy.fUseMC),
74   fPDG(0),
75   fCh(0),
76   fIM(0.0),
77   fPt(0.0),
78   fY(0.0),
79   fEta(0.0),
80   fMaxVz(copy.fMaxVz),
81   fMaxITSband(copy.fMaxITSband),
82   fTPCpLimit(copy.fTPCpLimit),
83   fMinTPCband(copy.fMinTPCband),
84   fMaxTPCband(copy.fMaxTPCband),
85   fRsnTreeComp(0x0),
86   fRsnTreeTrue(0x0),
87   fOutList(0x0),
88   fHEvents(0x0),
89   fESDtrackCutsTPC(copy.fESDtrackCutsTPC),
90   fESDtrackCutsITS(copy.fESDtrackCutsITS),
91   fESDpid(0x0),
92   fTOFmaker(0x0),
93   fTOFcalib(0x0),
94   fTOFcalibrateESD(kFALSE),
95   fTOFcorrectTExp(kFALSE),
96   fTOFuseT0(kFALSE),
97   fTOFtuneMC(kFALSE),
98   fTOFresolution(0.0)
99 {
100 //
101 // Copy constructor
102 //
103 }
104
105 //__________________________________________________________________________________________________
106 AliRsnAnalysisPhi7TeVNoPID& AliRsnAnalysisPhi7TeVNoPID::operator=(const AliRsnAnalysisPhi7TeVNoPID& copy)
107 {
108 //
109 // Assignment operator
110 //
111
112   fUseMC = copy.fUseMC;
113
114   fMaxVz   = copy.fMaxVz;
115   fMaxITSband = copy.fMaxITSband;
116   
117   fTPCpLimit  = copy.fTPCpLimit;
118   fMinTPCband = copy.fMinTPCband;
119   fMaxTPCband = copy.fMaxTPCband;
120   
121   fESDtrackCutsTPC = copy.fESDtrackCutsTPC;
122   fESDtrackCutsITS = copy.fESDtrackCutsITS;
123   
124   fTOFcalibrateESD = copy.fTOFcalibrateESD;
125   fTOFcorrectTExp = copy.fTOFcorrectTExp;
126   fTOFuseT0 = copy.fTOFuseT0;
127   fTOFtuneMC = copy.fTOFtuneMC;
128   fTOFresolution = copy.fTOFresolution;
129
130   return (*this);
131 }
132
133 //__________________________________________________________________________________________________
134 AliRsnAnalysisPhi7TeVNoPID::~AliRsnAnalysisPhi7TeVNoPID()
135 {
136 //
137 // Destructor
138 //
139
140   if (fRsnTreeComp) delete fRsnTreeComp;
141   if (fRsnTreeTrue) delete fRsnTreeTrue;
142   if (fHEvents)     delete fHEvents;
143   if (fESDpid)      delete fESDpid;
144 }
145
146 //__________________________________________________________________________________________________
147 void AliRsnAnalysisPhi7TeVNoPID::UserCreateOutputObjects()
148 {
149 //
150 // Create the output data container
151 //
152
153   // setup TPC response
154   fESDpid = new AliESDpid;
155   fESDpid->GetTPCResponse().SetBetheBlochParameters(fTPCpar[0],fTPCpar[1],fTPCpar[2],fTPCpar[3],fTPCpar[4]);
156
157   // setup TOF maker & calibration
158   fTOFcalib = new AliTOFcalib;
159   fTOFmaker = new AliTOFT0maker(fESDpid, fTOFcalib);
160   fTOFmaker->SetTimeResolution(fTOFresolution);
161   
162   // initialize random
163   gRandom->SetSeed(0);
164
165   // create output trees
166   OpenFile(1);
167   fRsnTreeComp = new TTree("rsnTree", "Pairs");
168
169   fRsnTreeComp->Branch("pdg", &fPDG, "pdg/S"   );
170   fRsnTreeComp->Branch("ch" , &fCh , "ch/S"    );
171   fRsnTreeComp->Branch("im" , &fIM , "im/F"    );
172   fRsnTreeComp->Branch("y"  , &fY  , "y/F"     );
173   fRsnTreeComp->Branch("pt" , &fPt , "pt/F"    );
174   fRsnTreeComp->Branch("eta", &fEta, "eta/F"   );
175   fRsnTreeComp->Branch("its", &fITS, "its[2]/S");
176   
177   fRsnTreeComp->Branch("p"     , &fP        , "p[2]/F");
178   fRsnTreeComp->Branch("ptpc"  , &fPTPC     , "ptpc[2]/F");
179   fRsnTreeComp->Branch("tpcpid", &fTPCnsigma, "tpcpid[2]/F");
180   fRsnTreeComp->Branch("itspid", &fITSnsigma, "itspid[2]/F");
181   fRsnTreeComp->Branch("tofpid", &fTOFdiff  , "tofpid[2]/F");
182
183   OpenFile(2);
184   fRsnTreeTrue = new TTree("rsnTrue", "True pairs");
185
186   fRsnTreeTrue->Branch("im" , &fIM , "im/F" );
187   fRsnTreeTrue->Branch("y"  , &fY  , "y/F"  );
188   fRsnTreeTrue->Branch("pt" , &fPt , "pt/F" );
189   fRsnTreeTrue->Branch("eta", &fEta, "eta/F");
190
191   OpenFile(3);
192   fOutList    = new TList;
193   fHEvents    = new TH1I("hEvents", "Event details", 5, 0, 5);
194   fVertexX[0] = new TH1F("hVertexTracksX", "X position of primary vertex (tracks)", 200,  -2,  2);
195   fVertexY[0] = new TH1F("hVertexTracksY", "Y position of primary vertex (tracks)", 200,  -2,  2);
196   fVertexZ[0] = new TH1F("hVertexTracksZ", "Z position of primary vertex (tracks)", 400, -40, 40);
197   fVertexX[1] = new TH1F("hVertexSPDX", "X position of primary vertex (SPD)", 200,  -2,  2);
198   fVertexY[1] = new TH1F("hVertexSPDY", "Y position of primary vertex (SPD)", 200,  -2,  2);
199   fVertexZ[1] = new TH1F("hVertexSPDZ", "Z position of primary vertex (SPD)", 400, -40, 40);
200   
201   fHEvents->GetXaxis()->SetBinLabel(1, "Good vertex with tracks");
202   fHEvents->GetXaxis()->SetBinLabel(2, "Good vertex with SPD");
203   fHEvents->GetXaxis()->SetBinLabel(3, "Far vertex with tracks");
204   fHEvents->GetXaxis()->SetBinLabel(4, "Far vertex with SPD");
205   fHEvents->GetXaxis()->SetBinLabel(5, "No good vertex");
206
207   fOutList->Add(fHEvents);
208   fOutList->Add(fVertexX[0]);
209   fOutList->Add(fVertexY[0]);
210   fOutList->Add(fVertexZ[0]);
211   fOutList->Add(fVertexX[1]);
212   fOutList->Add(fVertexY[1]);
213   fOutList->Add(fVertexZ[1]);
214 }
215
216 //__________________________________________________________________________________________________
217 void AliRsnAnalysisPhi7TeVNoPID::UserExec(Option_t *)
218 {
219 //
220 // Main execution function.
221 // Fills the fHEvents data member with the following legenda:
222 // 0 -- event OK, prim vertex with tracks
223 // 1 -- event OK, prim vertex with SPD
224 // 2 -- event OK but vz large
225 // 3 -- event bad
226 //
227
228   static Int_t evNum = 0;
229   evNum++;
230
231   // retrieve ESD event and related stack (if available)
232   AliESDEvent *esd   = dynamic_cast<AliESDEvent*>(fInputEvent);
233   AliStack    *stack = (fMCEvent ? fMCEvent->Stack() : 0x0);
234   
235   // check the event
236   Int_t eval = EventEval(esd);
237   fHEvents->Fill(eval);
238   
239   // if the event is good for analysis, process it
240   if (eval == kGoodTracksPrimaryVertex || eval == kGoodSPDPrimaryVertex)
241   {
242     ProcessESD(esd, stack);
243     ProcessMC(stack);
244   }
245   
246   // update histogram container
247   PostData(3, fOutList);
248 }
249
250 //__________________________________________________________________________________________________
251 void AliRsnAnalysisPhi7TeVNoPID::Terminate(Option_t *)
252 {
253 //
254 // Terminate
255 //
256 }
257
258 //__________________________________________________________________________________________________
259 Int_t AliRsnAnalysisPhi7TeVNoPID::EventEval(AliESDEvent *esd)
260 {
261 //
262 // Checks if the event is good for analysis.
263 // Returns one of the flag values defined in the header
264 //
265
266   static Int_t evNum = 0;
267   evNum++;
268
269   // debug message
270   AliDebug(AliLog::kDebug + 1, Form("Event %d -- number of tracks = %d", evNum, esd->GetNumberOfTracks()));
271   
272   // get the best primary vertex:
273   // first try the one with tracks
274   const AliESDVertex *vTrk  = esd->GetPrimaryVertexTracks();
275   const AliESDVertex *vSPD  = esd->GetPrimaryVertexSPD();
276   Double_t            vzTrk = 1000.0;
277   Double_t            vzSPD = 1000.0;
278   if (vTrk) vzTrk = TMath::Abs(vTrk->GetZv());
279   if (vSPD) vzSPD = TMath::Abs(vSPD->GetZv());
280   AliDebug(AliLog::kDebug + 1, Form("Event %d -- vertex with tracks: contributors = %d, abs(vz) = %f", evNum, vTrk->GetNContributors(), vzTrk));
281   AliDebug(AliLog::kDebug + 1, Form("Event %d -- vertex with SPD,    contributors = %d, abs(vz) = %f", evNum, vSPD->GetNContributors(), vzSPD));
282   if(vTrk->GetNContributors() > 0)
283   {
284     // fill the histograms
285     fVertexX[0]->Fill(vTrk->GetXv());
286     fVertexY[0]->Fill(vTrk->GetYv());
287     fVertexZ[0]->Fill(vTrk->GetZv());
288     
289     // check VZ position
290     if (vzTrk <= fMaxVz)
291       return kGoodTracksPrimaryVertex;
292     else
293       return kFarTracksPrimaryVertex;
294   }
295   else if (vSPD->GetNContributors() > 0)
296   {
297     // fill the histograms
298     fVertexX[1]->Fill(vSPD->GetXv());
299     fVertexY[1]->Fill(vSPD->GetYv());
300     fVertexZ[1]->Fill(vSPD->GetZv());
301     
302     // check VZ position
303     if (vzSPD <= fMaxVz)
304       return kGoodSPDPrimaryVertex;
305     else
306       return kFarSPDPrimaryVertex;
307   }
308   else
309     return kNoGoodPrimaryVertex;
310 }
311
312 //__________________________________________________________________________________________________
313 void AliRsnAnalysisPhi7TeVNoPID::ProcessESD
314 (AliESDEvent *esd, AliStack *stack)
315 {
316 //
317 // This function works with the ESD object
318 //
319
320   // ITS stuff #1 create the response function
321   Bool_t isMC = (stack != 0x0);
322   AliITSPIDResponse itsrsp(isMC);
323
324   // TOF stuff #1: init OCDB
325   Int_t run = esd->GetRunNumber();
326   AliCDBManager *cdb = AliCDBManager::Instance();
327   cdb->SetDefaultStorage("raw://");
328   cdb->SetRun(run);
329   // TOF stuff #2: init calibration
330   fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
331   fTOFcalib->Init();
332   // TOF stuff #3: calibrate
333   if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(esd);
334   if (fTOFtuneMC) fTOFmaker->TuneForMC(esd);
335   if (fTOFuseT0) 
336   {
337     fTOFmaker->ComputeT0TOF(esd);
338     fTOFmaker->ApplyT0TOF(esd);
339     fESDpid->MakePID(esd, kFALSE, 0.);
340   }
341
342   // prepare to look on all tracks to select the ones
343   // which pass all the cuts
344   Int_t   ntracks = esd->GetNumberOfTracks();
345   TArrayI pos(ntracks);
346   TArrayI neg(ntracks);
347   TArrayI itspos(ntracks);
348   TArrayI itsneg(ntracks);
349   
350   // loop on all tracks
351   ULong_t  status;
352   Int_t    i, k, charge, npos = 0, nneg = 0, nITS;
353   Double_t times[10], itsSignal, mom, tofTime, tofRef;
354   Bool_t   isTPC, isITSSA;
355   UChar_t  itsCluMap;
356   for (i = 0; i < ntracks; i++)
357   {
358     AliESDtrack *track = esd->GetTrack(i);
359     if (!track) continue;
360     
361     // get commonly used variables
362     status  = (ULong_t)track->GetStatus();
363     mom     = track->P();
364     isTPC   = ((status & AliESDtrack::kTPCin)  != 0);
365     isITSSA = ((status & AliESDtrack::kTPCin)  == 0 && (status & AliESDtrack::kITSrefit) != 0 && (status & AliESDtrack::kITSpureSA) == 0 && (status & AliESDtrack::kITSpid) != 0);
366     
367     // accept only tracks which are TPC+ITS or ITS standalone
368     if (!isTPC && !isITSSA) continue;
369     
370     // check specific cuts for TPC and ITS-SA tracks
371     if (isTPC)
372     {
373       if (!fESDtrackCutsTPC.IsSelected(track)) continue;
374     }
375     else
376     {
377        if (!fESDtrackCutsITS.IsSelected(track)) continue;
378        itsCluMap = track->GetITSClusterMap();
379        nITS      = 0;
380        for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
381        if (nITS < 3) continue;
382     }
383
384     // if all checks are passed, add the track index in one of the
385     // charged tracks arrays
386     charge = (Int_t)track->Charge();
387     if (charge > 0)
388     {
389       pos[npos] = i;
390       if (isITSSA) itspos[npos] = 1; else itspos[npos] = 0;
391       npos++;
392     }
393     else if (charge < 0)
394     {
395       neg[nneg] = i;
396       if (isITSSA) itsneg[nneg] = 1; else itsneg[nneg] = 0;
397       nneg++;
398     }
399   }
400   
401   // resize arrays accordingly
402   pos.Set(npos);
403   neg.Set(nneg);
404   itspos.Set(npos);
405   itsneg.Set(nneg);
406
407   // loop on unlike-sign pairs to compute invariant mass signal
408   Int_t           ip, in, lp, ln;
409   AliPID          pid;
410   Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
411   Double_t        phimass = 1.019455;
412   TParticle      *partp = 0x0, *partn = 0x0;
413   AliESDtrack    *tp = 0x0, *tn = 0x0;
414   TLorentzVector  vp, vn, vsum, vref;
415   for (ip = 0; ip < npos; ip++)
416   {
417     tp = esd->GetTrack(pos[ip]);
418     lp = TMath::Abs(tp->GetLabel());
419     if (stack) partp = stack->Particle(lp);
420
421     for (in = 0; in < nneg; in++)
422     {
423       if (pos[ip] == neg[in]) 
424       {
425         AliError("POS = NEG");
426         continue;
427       }
428       tn = esd->GetTrack(neg[in]);
429       ln = TMath::Abs(tn->GetLabel());
430       if (stack) partn = stack->Particle(ln);
431
432       fPDG = 0;
433       if (partp && partn)
434       {
435         if (partp->GetFirstMother() == partn->GetFirstMother())
436         {
437           if (partp->GetFirstMother() > 0)
438           {
439             TParticle *mum = stack->Particle(partp->GetFirstMother());
440             fPDG = mum->GetPdgCode();
441           }
442         }
443       }
444       fPDG = TMath::Abs(fPDG);
445
446       vp.SetXYZM(tp->Px(), tp->Py(), tp->Pz(), kmass);
447       vn.SetXYZM(tn->Px(), tn->Py(), tn->Pz(), kmass);
448       vsum = vp + vn;
449       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
450
451       fCh     = 0;
452       fIM     = (Float_t)vsum.M();
453       fPt     = (Float_t)vsum.Perp();
454       fEta    = (Float_t)vsum.Eta();
455       fY      = (Float_t)vref.Rapidity();
456       fITS[0] = itspos[ip];
457       fITS[1] = itsneg[in];
458
459       if (fIM < 0.9 || fIM >  5.0) continue;
460       if (fPt < 0.0 || fPt > 20.0) continue;
461       
462       // PID signal for track #1
463       // here it is enough to check if it is a TPC track
464       // since we excluded the case that it is neither a TPC+ITS nor an ITS-SA
465       if ((tp->GetStatus() & AliESDtrack::kTPCin) != 0)
466       {
467         fP        [0] = tp->P();
468         fPTPC     [0] = tp->GetInnerParam()->P();
469         fTPCnsigma[0] = TMath::Abs(fESDpid->NumberOfSigmasTPC(tp, AliPID::kKaon));
470         fITSnsigma[0] = 1E6;
471         fTOFdiff  [0] = 1E6;
472         // check TOF (only if momentum is large than function asymptote and flags are OK)
473         if (((tp->GetStatus() & AliESDtrack::kTOFout) != 0) && ((tp->GetStatus() & AliESDtrack::kTIME) != 0))
474         {
475           tp->GetIntegratedTimes(times);
476           tofTime  = (Double_t)tp->GetTOFsignal();
477           tofRef   = times[AliPID::kKaon];
478           if (tofRef > 0.0) fTOFdiff[0] = (tofTime - tofRef) / tofRef;
479         }
480       }
481       else
482       {
483         fP        [0] = tp->P();
484         fPTPC     [0] = 1E6;
485         fTPCnsigma[0] = 1E6;
486         fTOFdiff  [0] = 1E6;
487         // check dE/dx
488         itsSignal = tp->GetITSsignal();
489         itsCluMap = tp->GetITSClusterMap();
490         nITS      = 0;
491         for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
492         fITSnsigma[0] = itsrsp.GetNumberOfSigmas(fP[0], itsSignal, AliPID::kKaon, nITS, kTRUE);
493       }
494       
495       // PID signal for track #2
496       // here it is enough to check if it is a TPC track
497       // since we excluded the case that it is neither a TPC+ITS nor an ITS-SA
498       if ((tp->GetStatus() & AliESDtrack::kTPCin) != 0)
499       {
500         fP        [1] = tn->P();
501         fPTPC     [1] = tn->GetInnerParam()->P();
502         fTPCnsigma[1] = TMath::Abs(fESDpid->NumberOfSigmasTPC(tn, AliPID::kKaon));
503         fITSnsigma[1] = 1E6;
504         fTOFdiff  [1] = 1E6;
505         // check TOF (only if momentum is large than function asymptote and flags are OK)
506         if (((tn->GetStatus() & AliESDtrack::kTOFout) != 0) && ((tn->GetStatus() & AliESDtrack::kTIME) != 0))
507         {
508           tn->GetIntegratedTimes(times);
509           tofTime  = (Double_t)tn->GetTOFsignal();
510           tofRef   = times[AliPID::kKaon];
511           if (tofRef > 0.0) fTOFdiff[1] = (tofTime - tofRef) / tofRef;
512         }
513       }
514       else
515       {
516         fP        [1] = tn->P();
517         fPTPC     [1] = 1E6;
518         fTPCnsigma[1] = 1E6;
519         fTOFdiff  [1] = 1E6;
520         // check dE/dx
521         itsSignal = tn->GetITSsignal();
522         itsCluMap = tn->GetITSClusterMap();
523         nITS      = 0;
524         for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
525         fITSnsigma[1] = itsrsp.GetNumberOfSigmas(fP[1], itsSignal, AliPID::kKaon, nITS, kTRUE);
526       }
527       
528       fRsnTreeComp->Fill();
529     }
530   }
531   
532   // loop on like-sign pairs to compute invariant mass background
533   Int_t           i1, i2;
534   AliESDtrack    *t1 = 0x0, *t2 = 0x0;
535   TLorentzVector  v1, v2;
536   
537   // pos-pos
538   for (i1 = 0; i1 < npos; i1++)
539   {
540     t1 = esd->GetTrack(pos[i1]);
541
542     for (i2 = i1+1; i2 < npos; i2++)
543     {
544       t2 = esd->GetTrack(pos[i2]);
545
546       v1.SetXYZM(t1->Px(), t1->Py(), t1->Pz(), kmass);
547       v2.SetXYZM(t2->Px(), t2->Py(), t2->Pz(), kmass);
548       vsum = v1 + v2;
549       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
550
551       fPDG    = 0;
552       fCh     = 1;
553       fIM     = (Float_t)vsum.M();
554       fPt     = (Float_t)vsum.Perp();
555       fEta    = (Float_t)vsum.Eta();
556       fY      = (Float_t)vref.Rapidity();
557       fITS[0] = itspos[i1];
558       fITS[1] = itspos[i2];
559
560       if (fIM < 0.9 || fIM >  5.0) continue;
561       if (fPt < 0.0 || fPt > 20.0) continue;
562       
563       // PID signal for track #1
564       // here it is enough to check if it is a TPC track
565       // since we excluded the case that it is neither a TPC+ITS nor an ITS-SA
566       if ((t1->GetStatus() & AliESDtrack::kTPCin) != 0)
567       {
568         fP        [0] = t1->P();
569         fPTPC     [0] = t1->GetInnerParam()->P();
570         fTPCnsigma[0] = TMath::Abs(fESDpid->NumberOfSigmasTPC(t1, AliPID::kKaon));
571         fITSnsigma[0] = 1E6;
572         fTOFdiff  [0] = 1E6;
573         // check TOF (only if momentum is large than function asymptote and flags are OK)
574         if (((t1->GetStatus() & AliESDtrack::kTOFout) != 0) && ((t1->GetStatus() & AliESDtrack::kTIME) != 0))
575         {
576           t1->GetIntegratedTimes(times);
577           tofTime  = (Double_t)t1->GetTOFsignal();
578           tofRef   = times[AliPID::kKaon];
579           if (tofRef > 0.0) fTOFdiff[0] = (tofTime - tofRef) / tofRef;
580         }
581       }
582       else
583       {
584         fP        [0] = t1->P();
585         fPTPC     [0] = 1E6;
586         fTPCnsigma[0] = 1E6;
587         fTOFdiff  [0] = 1E6;
588         // check dE/dx
589         itsSignal = t1->GetITSsignal();
590         itsCluMap = t1->GetITSClusterMap();
591         nITS      = 0;
592         for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
593         fITSnsigma[0] = itsrsp.GetNumberOfSigmas(fP[0], itsSignal, AliPID::kKaon, nITS, kTRUE);
594       }
595       
596       // PID signal for track #2
597       // here it is enough to check if it is a TPC track
598       // since we excluded the case that it is neither a TPC+ITS nor an ITS-SA
599       if ((t1->GetStatus() & AliESDtrack::kTPCin) != 0)
600       {
601         fP        [1] = t2->P();
602         fPTPC     [1] = t2->GetInnerParam()->P();
603         fTPCnsigma[1] = TMath::Abs(fESDpid->NumberOfSigmasTPC(t2, AliPID::kKaon));
604         fITSnsigma[1] = 1E6;
605         fTOFdiff  [1] = 1E6;
606         // check TOF (only if momentum is large than function asymptote and flags are OK)
607         if (((t2->GetStatus() & AliESDtrack::kTOFout) != 0) && ((t2->GetStatus() & AliESDtrack::kTIME) != 0))
608         {
609           t2->GetIntegratedTimes(times);
610           tofTime  = (Double_t)t2->GetTOFsignal();
611           tofRef   = times[AliPID::kKaon];
612           if (tofRef > 0.0) fTOFdiff[1] = (tofTime - tofRef) / tofRef;
613         }
614       }
615       else
616       {
617         fP        [1] = t2->P();
618         fPTPC     [1] = 1E6;
619         fTPCnsigma[1] = 1E6;
620         fTOFdiff  [1] = 1E6;
621         // check dE/dx
622         itsSignal = t2->GetITSsignal();
623         itsCluMap = t2->GetITSClusterMap();
624         nITS      = 0;
625         for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
626         fITSnsigma[1] = itsrsp.GetNumberOfSigmas(fP[1], itsSignal, AliPID::kKaon, nITS, kTRUE);
627       }
628       
629       //fRsnTreeComp->Fill();
630     }
631   }
632   // neg-neg
633   for (i1 = 0; i1 < nneg; i1++)
634   {
635     t1 = esd->GetTrack(neg[i1]);
636
637     for (i2 = i1+1; i2 < nneg; i2++)
638     {
639       t2 = esd->GetTrack(neg[i2]);
640
641       v1.SetXYZM(t1->Px(), t1->Py(), t1->Pz(), kmass);
642       v2.SetXYZM(t2->Px(), t2->Py(), t2->Pz(), kmass);
643       vsum = v1 + v2;
644       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
645
646       fPDG    = 0;
647       fCh     = -1;
648       fIM     = (Float_t)vsum.M();
649       fPt     = (Float_t)vsum.Perp();
650       fEta    = (Float_t)vsum.Eta();
651       fY      = (Float_t)vref.Rapidity();
652       fITS[0] = itsneg[i1];
653       fITS[1] = itsneg[i2];
654
655       if (fIM < 0.9 || fIM >  5.0) continue;
656       if (fPt < 0.0 || fPt > 20.0) continue;
657       
658       // PID signal for track #1
659       // here it is enough to check if it is a TPC track
660       // since we excluded the case that it is neither a TPC+ITS nor an ITS-SA
661       if ((t1->GetStatus() & AliESDtrack::kTPCin) != 0)
662       {
663         fP        [0] = t1->P();
664         fPTPC     [0] = t1->GetInnerParam()->P();
665         fTPCnsigma[0] = TMath::Abs(fESDpid->NumberOfSigmasTPC(t1, AliPID::kKaon));
666         fITSnsigma[0] = 1E6;
667         fTOFdiff  [0] = 1E6;
668         // check TOF (only if momentum is large than function asymptote and flags are OK)
669         if (((t1->GetStatus() & AliESDtrack::kTOFout) != 0) && ((t1->GetStatus() & AliESDtrack::kTIME) != 0))
670         {
671           t1->GetIntegratedTimes(times);
672           tofTime  = (Double_t)t1->GetTOFsignal();
673           tofRef   = times[AliPID::kKaon];
674           if (tofRef > 0.0) fTOFdiff[0] = (tofTime - tofRef) / tofRef;
675         }
676       }
677       else
678       {
679         fP        [0] = t1->P();
680         fPTPC     [0] = 1E6;
681         fTPCnsigma[0] = 1E6;
682         fTOFdiff  [0] = 1E6;
683         // check dE/dx
684         itsSignal = t1->GetITSsignal();
685         itsCluMap = t1->GetITSClusterMap();
686         nITS      = 0;
687         for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
688         fITSnsigma[0] = itsrsp.GetNumberOfSigmas(fP[0], itsSignal, AliPID::kKaon, nITS, kTRUE);
689       }
690       
691       // PID signal for track #2
692       // here it is enough to check if it is a TPC track
693       // since we excluded the case that it is neither a TPC+ITS nor an ITS-SA
694       if ((t1->GetStatus() & AliESDtrack::kTPCin) != 0)
695       {
696         fP        [1] = t2->P();
697         fPTPC     [1] = t2->GetInnerParam()->P();
698         fTPCnsigma[1] = TMath::Abs(fESDpid->NumberOfSigmasTPC(t2, AliPID::kKaon));
699         fITSnsigma[1] = 1E6;
700         fTOFdiff  [1] = 1E6;
701         // check TOF (only if momentum is large than function asymptote and flags are OK)
702         if (((t2->GetStatus() & AliESDtrack::kTOFout) != 0) && ((t2->GetStatus() & AliESDtrack::kTIME) != 0))
703         {
704           t2->GetIntegratedTimes(times);
705           tofTime  = (Double_t)t2->GetTOFsignal();
706           tofRef   = times[AliPID::kKaon];
707           if (tofRef > 0.0) fTOFdiff[1] = (tofTime - tofRef) / tofRef;
708         }
709       }
710       else
711       {
712         fP        [1] = t2->P();
713         fPTPC     [1] = 1E6;
714         fTPCnsigma[1] = 1E6;
715         fTOFdiff  [1] = 1E6;
716         // check dE/dx
717         itsSignal = t2->GetITSsignal();
718         itsCluMap = t2->GetITSClusterMap();
719         nITS      = 0;
720         for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
721         fITSnsigma[1] = itsrsp.GetNumberOfSigmas(fP[1], itsSignal, AliPID::kKaon, nITS, kTRUE);
722       }
723       
724       //fRsnTreeComp->Fill();
725     }
726   }
727
728   PostData(1, fRsnTreeComp);
729 }
730
731 //__________________________________________________________________________________________________
732 void AliRsnAnalysisPhi7TeVNoPID::ProcessMC(AliStack *stack)
733 {
734 //
735 // Function to process stack only
736 //
737
738   if (!stack) return;
739   Int_t nPart = stack->GetNtrack();
740
741   // loop to compute invariant mass
742   Int_t           ip, in;
743   AliPID          pid;
744   Double_t        kmass = pid.ParticleMass(AliPID::kKaon);
745   Double_t        phimass = 1.019455;
746   TParticle      *partp = 0x0, *partn = 0x0;
747   TLorentzVector  vp, vn, vsum, vref;
748
749   for (ip = 0; ip < nPart; ip++)
750   {
751     partp = stack->Particle(ip);
752     if (partp->GetPdgCode() != 321) continue;
753
754     for (in = 0; in < nPart; in++)
755     {
756       partn = stack->Particle(in);
757       if (partn->GetPdgCode() != -321) continue;
758
759       fPDG = 0;
760       if (partp->GetFirstMother() == partn->GetFirstMother())
761       {
762         if (partp->GetFirstMother() > 0)
763         {
764           TParticle *mum = stack->Particle(partp->GetFirstMother());
765           fPDG = mum->GetPdgCode();
766         }
767       }
768       fPDG = TMath::Abs(fPDG);
769       if (fPDG != 333) continue;
770
771       vp.SetXYZM(partp->Px(), partp->Py(), partp->Pz(), kmass);
772       vn.SetXYZM(partn->Px(), partn->Py(), partn->Pz(), kmass);
773       vsum = vp + vn;
774       vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
775
776       fIM  = (Float_t)vsum.M();
777       fPt  = (Float_t)vsum.Perp();
778       fEta = (Float_t)vsum.Eta();
779       fY   = (Float_t)vref.Rapidity();
780
781       fRsnTreeTrue->Fill();
782     }
783   }
784
785   PostData(2, fRsnTreeTrue);
786 }