]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWG2/RESONANCES/AliRsnAnalysisPhi7TeV.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
fixed sig.segv
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnAnalysisPhi7TeV.cxx
... / ...
CommitLineData
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//__________________________________________________________________________________________________
31AliRsnAnalysisPhi7TeV::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//__________________________________________________________________________________________________
74AliRsnAnalysisPhi7TeV::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//__________________________________________________________________________________________________
112AliRsnAnalysisPhi7TeV& 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//__________________________________________________________________________________________________
143AliRsnAnalysisPhi7TeV::~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//__________________________________________________________________________________________________
156void 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//__________________________________________________________________________________________________
220void 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//__________________________________________________________________________________________________
254void AliRsnAnalysisPhi7TeV::Terminate(Option_t *)
255{
256//
257// Terminate
258//
259}
260
261//__________________________________________________________________________________________________
262Int_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//__________________________________________________________________________________________________
315void 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//__________________________________________________________________________________________________
660void 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//__________________________________________________________________________________________________
717void 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}