]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/RESONANCES/AliRsnAnalysisPhi7TeV.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
fixed sig.segv
[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnAnalysisPhi7TeV.cxx
CommitLineData
6279d59e 1//
2// Implementation file for implementation of data analysis aft 900 GeV
3//
4// Author: A. Pulvirenti
5//
6
7#include "Riostream.h"
4c06ce33 8#include <iomanip>
6279d59e 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"
4c06ce33 26#include "AliITSPIDResponse.h"
6279d59e 27
28#include "AliRsnAnalysisPhi7TeV.h"
29
30//__________________________________________________________________________________________________
31AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const char *name) :
32 AliAnalysisTaskSE(name),
33 fUseMC(kFALSE),
0d73200d 34 fCheckITS(kTRUE),
35 fCheckTPC(kTRUE),
36 fCheckTOF(kTRUE),
6279d59e 37 fPDG(0),
4c06ce33 38 fCh(0),
6279d59e 39 fIM(0.0),
40 fPt(0.0),
41 fY(0.0),
42 fEta(0.0),
4c06ce33 43 fMaxVz(1E6),
44 fMaxITSband(1E6),
45 fTPCpLimit(0.35),
6279d59e 46 fMinTPCband(-1E6),
47 fMaxTPCband( 1E6),
48 fRsnTreeComp(0x0),
49 fRsnTreeTrue(0x0),
50 fOutList(0x0),
51 fHEvents(0x0),
4c06ce33 52 fESDtrackCutsTPC(),
53 fESDtrackCutsITS(),
6279d59e 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),
0d73200d 77 fCheckITS(copy.fCheckITS),
78 fCheckTPC(copy.fCheckTPC),
79 fCheckTOF(copy.fCheckTOF),
6279d59e 80 fPDG(0),
4c06ce33 81 fCh(0),
6279d59e 82 fIM(0.0),
83 fPt(0.0),
84 fY(0.0),
85 fEta(0.0),
4c06ce33 86 fMaxVz(copy.fMaxVz),
87 fMaxITSband(copy.fMaxITSband),
88 fTPCpLimit(copy.fTPCpLimit),
6279d59e 89 fMinTPCband(copy.fMinTPCband),
90 fMaxTPCband(copy.fMaxTPCband),
91 fRsnTreeComp(0x0),
92 fRsnTreeTrue(0x0),
93 fOutList(0x0),
94 fHEvents(0x0),
4c06ce33 95 fESDtrackCutsTPC(copy.fESDtrackCutsTPC),
96 fESDtrackCutsITS(copy.fESDtrackCutsITS),
6279d59e 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;
0d73200d 119 fCheckITS = copy.fCheckITS;
120 fCheckTPC = copy.fCheckTPC;
121 fCheckTOF = copy.fCheckTOF;
6279d59e 122
4c06ce33 123 fMaxVz = copy.fMaxVz;
124 fMaxITSband = copy.fMaxITSband;
125
126 fTPCpLimit = copy.fTPCpLimit;
6279d59e 127 fMinTPCband = copy.fMinTPCband;
128 fMaxTPCband = copy.fMaxTPCband;
129
4c06ce33 130 fESDtrackCutsTPC = copy.fESDtrackCutsTPC;
131 fESDtrackCutsITS = copy.fESDtrackCutsITS;
132
6279d59e 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;
6279d59e 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
4c06ce33 178 fRsnTreeComp->Branch("pdg", &fPDG, "pdg/S" );
37dcaaad 179 fRsnTreeComp->Branch("ch" , &fCh , "ch/S" );
4c06ce33 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" );
37dcaaad 184 fRsnTreeComp->Branch("its", &fITS, "its[2]/S");
6279d59e 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);
4c06ce33 195 fOutList = new TList;
196 fHEvents = new TH1I("hEvents", "Event details", 5, 0, 5);
37dcaaad 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);
4c06ce33 199 fVertexZ[0] = new TH1F("hVertexTracksZ", "Z position of primary vertex (tracks)", 400, -40, 40);
37dcaaad 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);
4c06ce33 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
6279d59e 210 fOutList->Add(fHEvents);
4c06ce33 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]);
6279d59e 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);
6279d59e 237
4c06ce33 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)
6279d59e 244 {
4c06ce33 245 ProcessESD(esd, stack);
246 ProcessMC(stack);
6279d59e 247 }
248
4c06ce33 249 // update histogram container
6279d59e 250 PostData(3, fOutList);
251}
252
253//__________________________________________________________________________________________________
254void AliRsnAnalysisPhi7TeV::Terminate(Option_t *)
255{
256//
257// Terminate
258//
259}
260
4c06ce33 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++;
4c06ce33 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;
0d73200d 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();
4c06ce33 282 if (vTrk) vzTrk = TMath::Abs(vTrk->GetZv());
283 if (vSPD) vzSPD = TMath::Abs(vSPD->GetZv());
0d73200d 284 if(vTrk && ncTrk > 0)
4c06ce33 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 }
0d73200d 297 else if (vSPD && ncSPD > 0)
4c06ce33 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
6279d59e 314//__________________________________________________________________________________________________
315void AliRsnAnalysisPhi7TeV::ProcessESD
4c06ce33 316(AliESDEvent *esd, AliStack *stack)
6279d59e 317{
318//
319// This function works with the ESD object
320//
321
0d73200d 322 static Int_t lastRun = -1;
323
4c06ce33 324 // ITS stuff #1 create the response function
325 Bool_t isMC = (stack != 0x0);
326 AliITSPIDResponse itsrsp(isMC);
327
6279d59e 328 // TOF stuff #1: init OCDB
329 Int_t run = esd->GetRunNumber();
0d73200d 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);
6279d59e 343 // TOF stuff #2: init calibration
0d73200d 344 //fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
345 //fTOFcalib->Init();
6279d59e 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 }
4c06ce33 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;
6279d59e 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);
6aecf4fd 366 TArrayI itspos(ntracks);
367 TArrayI itsneg(ntracks);
6279d59e 368
6279d59e 369 // loop on all tracks
4c06ce33 370 ULong_t status;
371 Int_t i, k, charge, npos = 0, nneg = 0, nITS;
0d73200d 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;
4c06ce33 374 UChar_t itsCluMap;
6279d59e 375 for (i = 0; i < ntracks; i++)
376 {
377 AliESDtrack *track = esd->GetTrack(i);
4c06ce33 378 if (!track) continue;
6279d59e 379
4c06ce33 380 // get commonly used variables
6aecf4fd 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);
0d73200d 385 //cout << isTPC << ' ' << isITSSA << endl;
6aecf4fd 386
387 // accept only tracks which are TPC+ITS or ITS standalone
388 if (!isTPC && !isITSSA) continue;
4c06ce33 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
37dcaaad 393 okTrack = kTRUE;
4c06ce33 394
6aecf4fd 395 if (isTPC)
6279d59e 396 {
4c06ce33 397 // check standard ESD cuts
0d73200d 398 okQuality = fESDtrackCutsTPC.IsSelected(track);
399 //cout << "GLOBAL -- quality = " << (okQuality ? "GOOD" : "BAD") << endl;
400 if (!okQuality) continue;
4c06ce33 401
402 // check TPC dE/dx
0d73200d 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 }
4c06ce33 415
416 // check TOF (only if momentum is large than function asymptote and flags are OK)
0d73200d 417 if (fCheckTOF)
6279d59e 418 {
0d73200d 419 if (((status & AliESDtrack::kTOFout) != 0) && ((status & AliESDtrack::kTIME) != 0) && mom > TMath::Max(b1, b2))
4c06ce33 420 {
0d73200d 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;
4c06ce33 444 }
6279d59e 445 }
0d73200d 446 else
447 {
448 okTOF = kTRUE;
449 }
450
451 okTrack = okQuality && okTPC && okTOF;
452 //cout << "GLOBAL -- overall = " << (okTrack ? "ACCEPTED" : "REJECTED") << endl;
6279d59e 453 }
6aecf4fd 454 else
4c06ce33 455 {
456 // check standard ESD cuts
0d73200d 457 okQuality = fESDtrackCutsITS.IsSelected(track);
458 //cout << "ITSSA -- quality = " << (okQuality ? "GOOD" : "BAD") << endl;
459 if (!okQuality) continue;
4c06ce33 460
4c06ce33 461 // check dE/dx
0d73200d 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;
4c06ce33 487 }
37dcaaad 488
489 // skip tracks not passing cuts
490 if (!okTrack) continue;
4c06ce33 491
492 // if all checks are passed, add the track index in one of the
493 // charged tracks arrays
6279d59e 494 charge = (Int_t)track->Charge();
495 if (charge > 0)
6aecf4fd 496 {
497 pos[npos] = i;
498 if (isITSSA) itspos[npos] = 1; else itspos[npos] = 0;
499 npos++;
500 }
6279d59e 501 else if (charge < 0)
6aecf4fd 502 {
503 neg[nneg] = i;
504 if (isITSSA) itsneg[nneg] = 1; else itsneg[nneg] = 0;
505 nneg++;
506 }
6279d59e 507 }
508
509 // resize arrays accordingly
510 pos.Set(npos);
511 neg.Set(nneg);
6aecf4fd 512 itspos.Set(npos);
513 itsneg.Set(nneg);
0d73200d 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
4c06ce33 533 // loop on unlike-sign pairs to compute invariant mass signal
6279d59e 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 {
4c06ce33 549 if (pos[ip] == neg[in])
550 {
551 AliError("POS = NEG");
552 continue;
553 }
6279d59e 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
37dcaaad 577 fCh = 0;
4c06ce33 578 fIM = (Float_t)vsum.M();
579 fPt = (Float_t)vsum.Perp();
580 fEta = (Float_t)vsum.Eta();
581 fY = (Float_t)vref.Rapidity();
6aecf4fd 582 fITS[0] = itspos[ip];
583 fITS[1] = itsneg[in];
4c06ce33 584
0d73200d 585 //if (fIM < 0.9 || fIM > 1.4) continue;
586 //if (fPt < 0.0 || fPt > 20.0) continue;
37dcaaad 587
4c06ce33 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;
37dcaaad 612 fCh = 1;
4c06ce33 613 fIM = (Float_t)vsum.M();
614 fPt = (Float_t)vsum.Perp();
615 fEta = (Float_t)vsum.Eta();
616 fY = (Float_t)vref.Rapidity();
6aecf4fd 617 fITS[0] = itspos[i1];
618 fITS[1] = itspos[i2];
4c06ce33 619
0d73200d 620 //if (fIM < 0.9 || fIM > 1.4) continue;
621 //if (fPt < 0.0 || fPt > 20.0) continue;
37dcaaad 622
4c06ce33 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;
37dcaaad 641 fCh = -1;
4c06ce33 642 fIM = (Float_t)vsum.M();
643 fPt = (Float_t)vsum.Perp();
644 fEta = (Float_t)vsum.Eta();
645 fY = (Float_t)vref.Rapidity();
6aecf4fd 646 fITS[0] = itsneg[i1];
647 fITS[1] = itsneg[i2];
6279d59e 648
0d73200d 649 //if (fIM < 0.9 || fIM > 1.4) continue;
650 //if (fPt < 0.0 || fPt > 20.0) continue;
37dcaaad 651
6279d59e 652 fRsnTreeComp->Fill();
653 }
654 }
655
656 PostData(1, fRsnTreeComp);
657}
658
0d73200d 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
6279d59e 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}