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