2 // Implementation file for implementation of data analysis aft 900 GeV
4 // Author: A. Pulvirenti
14 #include "TParticle.h"
16 #include "TLorentzVector.h"
19 #include "AliESDpid.h"
20 #include "AliESDEvent.h"
21 #include "AliESDVertex.h"
22 #include "AliESDtrack.h"
23 #include "AliESDtrackCuts.h"
25 #include "AliMCEvent.h"
26 #include "AliTOFT0maker.h"
27 #include "AliTOFcalib.h"
28 #include "AliCDBManager.h"
29 #include "AliITSPIDResponse.h"
31 #include "AliRsnEvent.h"
32 #include "AliRsnTarget.h"
33 #include "AliRsnDaughter.h"
34 #include "AliRsnCutESD2010.h"
36 #include "AliRsnAnalysisPhi7TeV.h"
38 //__________________________________________________________________________________________________
39 AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const char *name, Bool_t isMC) :
40 AliAnalysisTaskSE(name),
65 fTOFcalibrateESD(!isMC),
66 fTOFcorrectTExp(kTRUE),
69 fTOFresolution(100.0),
78 DefineOutput(1, TList::Class());
81 //__________________________________________________________________________________________________
82 AliRsnAnalysisPhi7TeV::AliRsnAnalysisPhi7TeV(const AliRsnAnalysisPhi7TeV& copy) :
83 AliAnalysisTaskSE(copy),
85 fCheckITS(copy.fCheckITS),
86 fCheckTPC(copy.fCheckTPC),
87 fCheckTOF(copy.fCheckTOF),
88 fAddITSSA(copy.fAddITSSA),
90 fMaxITSband(copy.fMaxITSband),
91 fMaxITSmom(copy.fMaxITSmom),
92 fTPCpLimit(copy.fTPCpLimit),
93 fMinTPCband(copy.fMinTPCband),
94 fMaxTPCband(copy.fMaxTPCband),
95 fMinTOF(copy.fMinTOF),
96 fMaxTOF(copy.fMaxTOF),
103 fESDtrackCutsTPC(copy.fESDtrackCutsTPC),
104 fESDtrackCutsITS(copy.fESDtrackCutsITS),
108 fTOFcalibrateESD(kFALSE),
109 fTOFcorrectTExp(kFALSE),
120 SetUseMC(copy.fUseMC);
123 //__________________________________________________________________________________________________
124 AliRsnAnalysisPhi7TeV& AliRsnAnalysisPhi7TeV::operator=(const AliRsnAnalysisPhi7TeV& copy)
127 // Assignment operator
130 fUseMC = copy.fUseMC;
131 fCheckITS = copy.fCheckITS;
132 fCheckTPC = copy.fCheckTPC;
133 fCheckTOF = copy.fCheckTOF;
134 fAddITSSA = copy.fAddITSSA;
136 fMaxVz = copy.fMaxVz;
137 fMaxITSband = copy.fMaxITSband;
138 fMaxITSmom = copy.fMaxITSmom;
140 fTPCpLimit = copy.fTPCpLimit;
141 fMinTPCband = copy.fMinTPCband;
142 fMaxTPCband = copy.fMaxTPCband;
144 fESDtrackCutsTPC = copy.fESDtrackCutsTPC;
145 fESDtrackCutsITS = copy.fESDtrackCutsITS;
147 fTOFcalibrateESD = copy.fTOFcalibrateESD;
148 fTOFcorrectTExp = copy.fTOFcorrectTExp;
149 fTOFuseT0 = copy.fTOFuseT0;
150 fTOFtuneMC = copy.fTOFtuneMC;
151 fTOFresolution = copy.fTOFresolution;
153 SetUseMC(copy.fUseMC);
158 //__________________________________________________________________________________________________
159 AliRsnAnalysisPhi7TeV::~AliRsnAnalysisPhi7TeV()
166 //_________________________________________________________________________________________________
167 void AliRsnAnalysisPhi7TeV::SetUseMC(Bool_t isMC)
170 // Sets some aspects of cuts depending on the fact that runs on MC or not
177 SetTPCpar(2.15898 / 50.0, 1.75295E1, 3.40030E-9, 1.96178, 3.91720);
178 SetTOFcalibrateESD(kFALSE);
183 SetTPCpar(1.41543 / 50.0, 2.63394E1, 5.0411E-11, 2.12543, 4.88663);
184 SetTOFcalibrateESD(kTRUE);
185 SetTOFtuneMC(kFALSE);
189 //__________________________________________________________________________________________________
190 void AliRsnAnalysisPhi7TeV::UserCreateOutputObjects()
193 // Create the output data container
197 //gRandom->SetSeed(0);
199 // create output trees
201 fOutList = new TList;
202 fHEvents = new TH1I("hEvents", "Event details", 6, 0, 6);
203 fVertexX[0] = new TH1F("hVertexTracksX", "X position of primary vertex (tracks)", 200, -2, 2);
204 fVertexY[0] = new TH1F("hVertexTracksY", "Y position of primary vertex (tracks)", 200, -2, 2);
205 fVertexZ[0] = new TH1F("hVertexTracksZ", "Z position of primary vertex (tracks)", 400, -40, 40);
206 fVertexX[1] = new TH1F("hVertexSPDX", "X position of primary vertex (SPD)", 200, -2, 2);
207 fVertexY[1] = new TH1F("hVertexSPDY", "Y position of primary vertex (SPD)", 200, -2, 2);
208 fVertexZ[1] = new TH1F("hVertexSPDZ", "Z position of primary vertex (SPD)", 400, -40, 40);
209 //fUnlike = new TH3F("hPM", "+- pairs", 500, 0.9, 1.4, 100, 0.0, 10.0, 24, -1.2, 1.2);
210 //fLikePP = new TH3F("hPP", "++ pairs", 500, 0.9, 1.4, 100, 0.0, 10.0, 24, -1.2, 1.2);
211 //fLikeMM = new TH3F("hMM", "-- pairs", 500, 0.9, 1.4, 100, 0.0, 10.0, 24, -1.2, 1.2);
212 //fTrues = new TH3F("hTR", "True pairs", 500, 0.9, 1.4, 100, 0.0, 10.0, 24, -1.2, 1.2);
213 TFile *ffile = TFile::Open("template.root");
216 TH3F *tmp = (TH3F*)ffile->Get("template");
217 fUnlike = (TH3F*)tmp->Clone("hPM");
218 fLikePP = (TH3F*)tmp->Clone("hPP");
219 fLikeMM = (TH3F*)tmp->Clone("hMM");
220 fTrues = (TH3F*)tmp->Clone("hTR");
223 fHEvents->GetXaxis()->SetBinLabel(1, "Good vertex with tracks");
224 fHEvents->GetXaxis()->SetBinLabel(2, "Good vertex with SPD");
225 fHEvents->GetXaxis()->SetBinLabel(3, "Far vertex with tracks");
226 fHEvents->GetXaxis()->SetBinLabel(4, "Far vertex with SPD");
227 fHEvents->GetXaxis()->SetBinLabel(5, "No good vertex");
228 fHEvents->GetXaxis()->SetBinLabel(6, "Empty event");
230 fOutList->Add(fHEvents);
231 fOutList->Add(fVertexX[0]);
232 fOutList->Add(fVertexY[0]);
233 fOutList->Add(fVertexZ[0]);
234 fOutList->Add(fVertexX[1]);
235 fOutList->Add(fVertexY[1]);
236 fOutList->Add(fVertexZ[1]);
237 fOutList->Add(fUnlike);
238 fOutList->Add(fLikePP);
239 fOutList->Add(fLikeMM);
240 fOutList->Add(fTrues);
242 // setup RSN-related objects
244 fRsnCuts.SetMC (fUseMC);
245 fRsnCuts.SetCheckITS (fCheckITS);
246 fRsnCuts.SetCheckTPC (fCheckTPC);
247 fRsnCuts.SetCheckTOF (fCheckTOF);
248 fRsnCuts.SetUseITSTPC(kTRUE);
249 fRsnCuts.SetUseITSSA (fAddITSSA);
250 fRsnCuts.SetPID (AliPID::kKaon);
251 fRsnCuts.SetMaxITSPIDmom(fMaxITSmom);
252 fRsnCuts.SetITSband(fMaxITSband);
253 fRsnCuts.SetTPCpLimit(fTPCpLimit);
254 fRsnCuts.SetTPCrange(fMinTPCband, fMaxTPCband);
255 fRsnCuts.SetTPCpar(fTPCpar[0],fTPCpar[1],fTPCpar[2],fTPCpar[3],fTPCpar[4]);
256 //fRsnCuts.SetTOFcalibrateESD(fTOFcalibrateESD);
257 fRsnCuts.SetTOFcorrectTExp (fTOFcorrectTExp);
258 fRsnCuts.SetTOFuseT0 (fTOFuseT0);
259 fRsnCuts.SetTOFtuneMC (fTOFtuneMC);
260 fRsnCuts.SetTOFresolution (fTOFresolution);
261 fRsnCuts.SetTOFrange (fMinTOF, fMaxTOF);
262 fRsnCuts.GetCutsTPC()->Copy(fESDtrackCutsTPC);
263 fRsnCuts.GetCutsITS()->Copy(fESDtrackCutsITS);
266 //__________________________________________________________________________________________________
267 void AliRsnAnalysisPhi7TeV::UserExec(Option_t *)
270 // Main execution function.
271 // Fills the fHEvents data member with the following legenda:
272 // 0 -- event OK, prim vertex with tracks
273 // 1 -- event OK, prim vertex with SPD
274 // 2 -- event OK but vz large
278 // retrieve ESD event and related stack (if available)
279 AliESDEvent *esd = dynamic_cast<AliESDEvent*>(fInputEvent);
280 AliStack *stack = (fMCEvent ? fMCEvent->Stack() : 0x0);
283 Int_t eval = EventEval(esd);
284 fHEvents->Fill(eval);
286 // if the event is good for analysis, process it
287 if (eval == kGoodTracksPrimaryVertex || eval == kGoodSPDPrimaryVertex)
289 ProcessESD(esd, stack);
292 // update histogram container
293 PostData(1, fOutList);
296 //__________________________________________________________________________________________________
297 void AliRsnAnalysisPhi7TeV::Terminate(Option_t *)
304 //__________________________________________________________________________________________________
305 Int_t AliRsnAnalysisPhi7TeV::EventEval(AliESDEvent *esd)
308 // Checks if the event is good for analysis.
309 // Returns one of the flag values defined in the header
312 static Int_t evNum = 0;
315 // reject empty events
316 Int_t ntracks = esd->GetNumberOfTracks();
317 if (!ntracks) return kEmptyEvent;
319 // get the best primary vertex:
320 // first try the one with tracks
321 const AliESDVertex *vTrk = esd->GetPrimaryVertexTracks();
322 const AliESDVertex *vSPD = esd->GetPrimaryVertexSPD();
323 Double_t vzTrk = 1000000.0;
324 Double_t vzSPD = 1000000.0;
327 if (vTrk) ncTrk = (Int_t)vTrk->GetNContributors();
328 if (vSPD) ncSPD = (Int_t)vSPD->GetNContributors();
329 if (vTrk) vzTrk = TMath::Abs(vTrk->GetZv());
330 if (vSPD) vzSPD = TMath::Abs(vSPD->GetZv());
331 if(vTrk && ncTrk > 0)
333 // fill the histograms
334 fVertexX[0]->Fill(vTrk->GetXv());
335 fVertexY[0]->Fill(vTrk->GetYv());
336 fVertexZ[0]->Fill(vTrk->GetZv());
340 return kGoodTracksPrimaryVertex;
342 return kFarTracksPrimaryVertex;
344 else if (vSPD && ncSPD > 0)
346 // fill the histograms
347 fVertexX[1]->Fill(vSPD->GetXv());
348 fVertexY[1]->Fill(vSPD->GetYv());
349 fVertexZ[1]->Fill(vSPD->GetZv());
353 return kGoodSPDPrimaryVertex;
355 return kFarSPDPrimaryVertex;
358 return kNoGoodPrimaryVertex;
361 //__________________________________________________________________________________________________
362 void AliRsnAnalysisPhi7TeV::ProcessESD
363 (AliESDEvent *esd, AliStack *stack)
366 // This function works with the ESD object
369 static Int_t lastRun = -1;
370 static Int_t evnum = 0;
374 Int_t run = esd->GetRunNumber();
376 // if absent, initialize ESD pid response
379 AliITSPIDResponse itsresponse(fUseMC);
381 fESDpid = new AliESDpid;
382 fESDpid->GetTPCResponse().SetBetheBlochParameters(fTPCpar[0],fTPCpar[1],fTPCpar[2],fTPCpar[3],fTPCpar[4]);
383 fESDpid->GetITSResponse() = itsresponse;
386 // if the run number has changed, update it now and give a message
389 AliInfo("============================================================================================");
390 AliInfo(Form("*** CHANGING RUN NUMBER: PREVIOUS = %d --> CURRENT = %d ***", lastRun, run));
391 AliInfo("============================================================================================");
394 AliCDBManager::Instance()->SetDefaultStorage("raw://");
395 AliCDBManager::Instance()->SetRun(lastRun);
397 if (fTOFmaker) delete fTOFmaker;
398 if (fTOFcalib) delete fTOFcalib;
400 fTOFcalib = new AliTOFcalib();
401 if (fTOFcorrectTExp) fTOFcalib->SetCorrectTExp(kTRUE);
404 fTOFmaker = new AliTOFT0maker(fESDpid, fTOFcalib);
405 fTOFmaker->SetTimeResolution(fTOFresolution);
408 if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(esd);
409 if (fTOFtuneMC) fTOFmaker->TuneForMC(esd);
412 fTOFmaker->ComputeT0TOF(esd);
413 fTOFmaker->ApplyT0TOF(esd);
414 fESDpid->MakePID(esd, kFALSE, 0.);
420 event.SetRefMC(fMCEvent);
421 //fRsnCuts.ProcessEvent(esd);
423 // loop on all tracks
424 Int_t i1, i2, label, pdg, ntracks = esd->GetNumberOfTracks();
427 Double_t kmass = pid.ParticleMass(AliPID::kKaon);
428 Double_t phimass = 1.019455;
429 Double_t angle, cosangle;
430 AliMCParticle *p1 = 0x0, *p2 = 0x0;
431 AliESDtrack *t1 = 0x0, *t2 = 0x0;
432 TLorentzVector v1, v2, vsum, vref;
434 // external loop (T1)
435 for (i1 = 0; i1 < ntracks; i1++)
437 t1 = esd->GetTrack(i1);
441 fDaughter.SetRef(t1);
444 p1 = (AliMCParticle*)fMCEvent->GetTrack(TMath::Abs(t1->GetLabel()));
445 if (p1) fDaughter.SetRefMC(p1);
447 fDaughter.SetMass(kmass);
448 v1.SetXYZM(t1->Px(), t1->Py(), t1->Pz(), kmass);
451 okTrack = OkTrack(t1, AliPID::kKaon);
453 // internal loop (T2)
454 for (i2 = i1+1; i2 < ntracks; i2++)
456 t2 = esd->GetTrack(i2);
460 if (!OkTrack(t2, AliPID::kKaon)) continue;
462 // if unlike pair, check if it is true
464 if ((t1->Charge() == t2->Charge()) && p1 && p2)
466 TParticle *part1 = p1->Particle();
467 TParticle *part2 = p2->Particle();
468 if (TMath::Abs(part1->GetPdgCode()) == 321 && TMath::Abs(part2->GetPdgCode()) == 321 && part1->GetFirstMother() == part2->GetFirstMother())
470 label = part1->GetFirstMother();
471 if (label >= 0 && label < stack->GetNtrack())
473 TParticle *mum = stack->Particle(label);
474 pdg = mum->GetPdgCode();
478 pdg = TMath::Abs(pdg);
481 v1.SetXYZM(t1->Px(), t1->Py(), t1->Pz(), kmass);
482 v2.SetXYZM(t2->Px(), t2->Py(), t2->Pz(), kmass);
483 angle = v1.Angle(v2.Vect());
484 cosangle = TMath::Abs(TMath::Cos(angle));
485 if (cosangle <= 0.02) continue;
487 vref.SetXYZM(vsum.X(), vsum.Y(), vsum.Z(), phimass);
489 // fill appropriate histogram
490 if (t1->Charge() != t2->Charge())
492 fUnlike->Fill(vsum.M(), vsum.Perp(), vref.Rapidity());
493 if (pdg == 333) fTrues->Fill(vsum.M(), vsum.Perp(), vref.Rapidity());
495 else if (t1->Charge() > 0)
497 fLikePP->Fill(vsum.M(), vsum.Perp(), vref.Rapidity());
501 fLikeMM->Fill(vsum.M(), vsum.Perp(), vref.Rapidity());
506 PostData(1, fOutList);
509 //______________________________________________________________________________
510 Bool_t AliRsnAnalysisPhi7TeV::OkQuality(AliESDtrack *track)
513 // Check track quality parameters.
514 // Rejects all tracks which are not either TPC+ITS nor ITS standalone.
515 // If tracks of any type are not flagged to be used, they are rejected anyway.
519 return fESDtrackCutsTPC.IsSelected(track);
520 else if (IsITSSA (track))
523 return fESDtrackCutsITS.IsSelected(track);
531 //______________________________________________________________________________
532 Bool_t AliRsnAnalysisPhi7TeV::OkITSPID (AliESDtrack *track, AliPID::EParticleType pid)
535 // Check ITS particle identification with 3sigma cut
538 // reject not ITS standalone tracks
539 if (!IsITSSA(track)) return kFALSE;
541 // count PID layers and reject if they are too few
542 Int_t k, nITSpidLayers = 0;
543 UChar_t itsCluMap = track->GetITSClusterMap();
544 for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITSpidLayers;
545 if (nITSpidLayers < 3)
547 AliDebug(AliLog::kDebug+2, "Rejecting track with too few ITS pid layers");
551 // check the track type (ITS+TPC or ITS standalone)
552 // and reject it if it is of none of the allowed types
553 Bool_t isSA = kFALSE;
554 if (IsITSTPC(track)) isSA = kFALSE;
555 else if (IsITSSA(track)) isSA = kTRUE;
558 AliWarning("Track is neither ITS+TPC nor ITS standalone");
562 // create the PID response object and compute nsigma
563 AliITSPIDResponse &itsrsp = fESDpid->GetITSResponse();
564 Double_t mom = track->P();
565 Double_t nSigma = itsrsp.GetNumberOfSigmas(mom, track->GetITSsignal(), pid, nITSpidLayers, isSA);
568 Bool_t ok = (TMath::Abs(nSigma) <= fMaxITSband);
571 AliDebug(AliLog::kDebug + 2, Form("ITS nsigma = %f -- max = %f -- cut %s", nSigma, fMaxITSband, (ok ? "passed" : "failed")));
577 //______________________________________________________________________________
578 Bool_t AliRsnAnalysisPhi7TeV::OkTPCPID (AliESDtrack *track, AliPID::EParticleType pid)
581 // Check TPC particle identification with {3|5}sigmacut,
582 // depending on the track total momentum.
585 // reject not TPC tracks
586 if (!IsITSTPC(track)) return kFALSE;
588 // setup TPC PID response
589 AliTPCPIDResponse &tpcrsp = fESDpid->GetTPCResponse();
590 tpcrsp.SetBetheBlochParameters(fTPCpar[0],fTPCpar[1],fTPCpar[2],fTPCpar[3],fTPCpar[4]);
592 // get momentum and number of sigmas and choose the reference band
593 Double_t mom = track->GetInnerParam()->P();
594 Double_t nSigma = tpcrsp.GetNumberOfSigmas(mom, track->GetTPCsignal(), track->GetTPCsignalN(), pid);
595 Double_t maxNSigma = fMaxTPCband;
596 if (mom < fTPCpLimit) maxNSigma = fMinTPCband;
599 Bool_t ok = (TMath::Abs(nSigma) <= maxNSigma);
602 AliDebug(AliLog::kDebug + 2, Form("TPC nsigma = %f -- max = %f -- cut %s", nSigma, maxNSigma, (ok ? "passed" : "failed")));
608 //______________________________________________________________________________
609 Bool_t AliRsnAnalysisPhi7TeV::OkTOFPID (AliESDtrack *track, AliPID::EParticleType pid)
612 // Check TOF particle identification if matched there.
615 // reject not TOF-matched tracks
616 if (!MatchTOF(track)) return kFALSE;
618 // setup TOF PID response
619 AliTOFPIDResponse &tofrsp = fESDpid->GetTOFResponse();
621 // get info for computation
622 Double_t momentum = track->P();
623 Double_t time = track->GetTOFsignal();
624 Double_t timeint[AliPID::kSPECIES];
625 tofrsp.GetStartTime(momentum);
626 track->GetIntegratedTimes(timeint);
629 Double_t timeDiff = time - timeint[(Int_t)pid];
630 Double_t sigmaRef = tofrsp.GetExpectedSigma(momentum, timeint[(Int_t)pid], AliPID::ParticleMass(pid));
631 Double_t nSigma = timeDiff / sigmaRef;
634 Bool_t ok = (nSigma >= fMinTOF && nSigma <= fMaxTOF);
637 AliDebug(AliLog::kDebug + 2, Form("TOF nsigma = %f -- range = %f - %f -- cut %s", nSigma, fMinTOF, fMaxTOF, (ok ? "passed" : "failed")));
638 //if (print) cout << Form("**PHI** TOF nsigma = %f -- range = %f - %f -- cut %s", nSigma, fMinTOF, fMaxTOF, (ok ? "passed" : "failed")) << endl;
644 //______________________________________________________________________________
645 Bool_t AliRsnAnalysisPhi7TeV::OkTrack(AliESDtrack *track, AliPID::EParticleType pid)
648 // Global track cut check
651 Bool_t okITS, okTPC, okTOF;
653 // check quality and track type and reject tracks not passing this step
654 if (!OkQuality(track))
656 AliDebug(AliLog::kDebug+2, "Failed quality cut");
657 //printf("[PHI] Track %4d: REJECTED\n", i);
661 // ITS PID can be checked always
662 // if PID is not required, the flag is sed as
663 // if the cut was alsways passed
664 okITS = OkITSPID(track, AliPID::kKaon);
665 if (!fCheckITS) okITS = kTRUE;
667 // TPC PID can be checked only for TPC+ITS tracks
668 // if PID is not required, the flag is sed as
669 // if the cut was alsways passed
671 if (IsITSTPC(track)) okTPC = OkTPCPID(track, AliPID::kKaon);
672 if (!fCheckTPC) okTPC = kTRUE;
674 // TOF PID can be checked only if TOF is matched
675 // if PID is not required, the flag is sed as
676 // if the cut was alsways passed
678 if (IsITSTPC(track) && MatchTOF(track)) okTOF = OkTOFPID(track, AliPID::kKaon);
679 if (!fCheckTOF) okTOF = kTRUE;
681 // now combine all outcomes according to the different possibilities:
682 // -- ITS standalone:
683 // --> only ITS PID, always
685 // --> ITS PID, only for momenta lower than 'fMaxITSPIDmom' and when the ITSpid flag is active
686 // --> TPC PID, always --> MASTER (first to be checked, if fails, track is rejected)
687 // --> TOF PID, only if matched
692 AliDebug(AliLog::kDebug+2, "ITS standalone track --> ITS PID failed");
693 //printf("[PHI] Track %4d: REJECTED\n", i);
697 else // checking IsITSTPC() is redundant due to OkQuality() cut check
701 AliDebug(AliLog::kDebug+2, "ITS+TPC track --> TPC PID failed");
702 //printf("[PHI] Track %4d: REJECTED\n", i);
705 else if (MatchTOF(track) && !okTOF)
707 AliDebug(AliLog::kDebug+2, "ITS+TPC track --> TOF matched but TOF PID failed");
708 //printf("[PHI] Track %4d: REJECTED\n", i);
711 else if (track->IsOn(AliESDtrack::kITSpid) && track->P() <= fMaxITSmom && !okITS)
713 AliDebug(AliLog::kDebug+2, Form("ITS+TPC track --> Momentum lower than limit (%.2f) and ITS PID failed", fMaxITSmom));
714 //printf("[PHI] Track %4d: REJECTED\n", i);
719 // this point is reached only if function didn't exit before due to somecheck not passed