8 #include "AliMCParticle.h"
9 //#include "AliStack.h"
11 #include "AliAnalysisTask.h"
12 #include "AliAnalysisManager.h"
14 #include "AliESDEvent.h"
15 #include "AliMCEvent.h"
16 #include "AliESDInputHandler.h"
17 #include "AliInputEventHandler.h"
19 #include "AliVVertex.h"
20 #include "AliAnalysisFilter.h"
22 #include "AliPIDResponse.h"
23 #include "AliTPCPIDResponse.h"
25 #include "AliTPCPIDEtaQA.h"
28 This task determines the eta dependence of the TPC signal.
29 For this purpose, only tracks fulfilling some standard quality cuts are taken into account.
30 The obtained data can be used to derive the functional behaviour of the eta dependence.
31 Such a function can be plugged into this task to correct for the eta dependence and to see
32 if there is then really no eta dependence left.
34 Class written by Benjamin Hess.
35 Contact: bhess@cern.ch
38 ClassImp(AliTPCPIDEtaQA)
40 //________________________________________________________________________
41 AliTPCPIDEtaQA::AliTPCPIDEtaQA()
43 , fPtThresholdForPhiCut(2.0)
44 , fPhiCutSecondBranchLow(0x0)
45 , fPhiCutSecondBranchHigh(0x0)
47 //OLD clusterQA, fhNumClustersPhiPrimePtBeforeCut(0x0)
48 //OLD clusterQA, fhNumClustersPhiPrimePtAfterCut(0x0)
49 , fhPhiPrimeCutEfficiency(0x0)
50 , fOutputContainer(0x0)
52 // default Constructor
54 // Question: Is this the right place to initialize these functions?
55 // Will it work on proof? i.e. will they be streamed to the workers?
56 // Also one should add getters and setters
57 //TODO fPhiCutLow = new TF1("StandardPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 50);
58 //TODO fPhiCutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 50);
60 //TODO NEW fPhiCutLow = new TF1("StandardPhiCutLow", "0.072/x+pi/18.0-0.035", 0, 50);
61 //TODO NEW fPhiCutHigh = new TF1("StandardPhiCutHigh", "0.07/x/x+0.1/x+pi/18.0+0.035", 0, 50);
63 fPhiCutSecondBranchLow = new TF1("StandardPhiCutSecondBranchLow", "NewStandardPhiCutLow - 2.*pi/18.", 0, 50);
64 fPhiCutSecondBranchHigh = new TF1("StandardPhiCutSecondBranchHigh", "0.07/x+pi/18.0+0.125 - 2.*pi/18.", 0, 50);
67 //________________________________________________________________________
68 AliTPCPIDEtaQA::AliTPCPIDEtaQA(const char *name)
70 , fPtThresholdForPhiCut(2.0)
71 , fPhiCutSecondBranchLow(0x0)
72 , fPhiCutSecondBranchHigh(0x0)
74 //OLD clusterQA, fhNumClustersPhiPrimePtBeforeCut(0x0)
75 //OLD clusterQA, fhNumClustersPhiPrimePtAfterCut(0x0)
76 , fhPhiPrimeCutEfficiency(0x0)
77 , fOutputContainer(0x0)
81 // Question: Is this the right place to initialize these functions?
82 // Will it work on proof? i.e. will they be streamed to the workers?
83 // Also one should add getters and setters
84 //TODO fPhiCutLow = new TF1("StandardPhiCutLow", "0.1/x/x+pi/18.0-0.025", 0, 50);
85 //TODO fPhiCutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 50);
87 //TODO NEW fPhiCutLow = new TF1("StandardPhiCutLow", "0.072/x+pi/18.0-0.035", 0, 50);
88 //TODO NEW fPhiCutHigh = new TF1("StandardPhiCutHigh", "0.07/x/x+0.1/x+pi/18.0+0.035", 0, 50);
90 fPhiCutSecondBranchLow = new TF1("StandardPhiCutSecondBranchLow", "StandardPhiCutLow - 2.*pi/18.", 0, 50);
91 fPhiCutSecondBranchHigh = new TF1("StandardPhiCutSecondBranchHigh", "0.07/x+pi/18.0+0.125 - 2.*pi/18.", 0, 50);
93 // Define input and output slots here
94 // Input slot #0 works with a TChain
95 DefineInput(0, TChain::Class());
96 // Output slot #1 writes into a TObjArray container
97 DefineOutput(1, TObjArray::Class());
101 //________________________________________________________________________
102 AliTPCPIDEtaQA::~AliTPCPIDEtaQA()
106 delete fOutputContainer;
107 fOutputContainer = 0;
109 delete fPhiCutSecondBranchLow;
110 fPhiCutSecondBranchLow = 0;
112 delete fPhiCutSecondBranchHigh;
113 fPhiCutSecondBranchHigh = 0;
117 //________________________________________________________________________
118 void AliTPCPIDEtaQA::UserCreateOutputObjects()
123 AliTPCPIDBase::UserCreateOutputObjects();
127 fOutputContainer = new TObjArray(1);
128 fOutputContainer->SetName(GetName()) ;
129 fOutputContainer->SetOwner(kTRUE);
131 const Int_t nEta = TMath::Nint(40 * fEtaCut);
133 const Int_t nPtBins = 68;
134 Double_t binsPt[nPtBins+1] = {0. , 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
135 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95,
136 1.0, 1.1 , 1.2, 1.3 , 1.4, 1.5 , 1.6, 1.7 , 1.8, 1.9 ,
137 2.0, 2.2 , 2.4, 2.6 , 2.8, 3.0 , 3.2, 3.4 , 3.6, 3.8 ,
138 4.0, 4.5 , 5.0, 5.5 , 6.0, 6.5 , 7.0, 8.0 , 9.0, 10.0,
139 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0, 22.0, 24.0,
140 26.0, 28.0, 30.0, 32.0, 34.0, 36.0, 40.0, 45.0, 50.0 };
143 const Int_t nBins = 6;
144 // MC PID, SelectSpecies, P(TPC_inner), multiplicity, deltaPrimeSpecies, eta
147 { 9, 4, nPtBins, 40, 201, nEta };
148 Double_t xmin[nBins] =
149 { 0., 0., 0., 0., 0.5, -fEtaCut};
150 Double_t xmax[nBins] =
151 { 9., 4., 50.0, 20000, 2.0, fEtaCut };
154 fhPIDdataAll = new THnSparseI("hPIDdataAll","", nBins, bins, xmin, xmax);
155 SetUpHist(fhPIDdataAll, binsPt);
156 fOutputContainer->Add(fhPIDdataAll);
159 const Int_t nBinsQA = 3;
162 Int_t binsQA[nBinsQA] =
164 Double_t xminQA[nBinsQA] =
166 Double_t xmaxQA[nBinsQA] =
167 { 50.0, TMath::Pi() / 9., 160 };
169 fhNumClustersPhiPrimePtBeforeCut = new THnSparseI("hNumClustersPhiPrimePtBeforeCut", "", nBinsQA, binsQA, xminQA, xmaxQA);
170 fhNumClustersPhiPrimePtBeforeCut->SetBinEdges(0, binsPt);
171 fhNumClustersPhiPrimePtBeforeCut->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
172 fhNumClustersPhiPrimePtBeforeCut->GetAxis(1)->SetTitle("#phi'");
173 fhNumClustersPhiPrimePtBeforeCut->GetAxis(2)->SetTitle("Ncl");
174 fOutputContainer->Add(fhNumClustersPhiPrimePtBeforeCut);
176 fhNumClustersPhiPrimePtAfterCut = new THnSparseI("hNumClustersPhiPrimePtAfterCut", "", nBinsQA, binsQA, xminQA, xmaxQA);
177 fhNumClustersPhiPrimePtAfterCut->SetBinEdges(0, binsPt);
178 fhNumClustersPhiPrimePtAfterCut->GetAxis(0)->SetTitle("p_{T} (GeV/c)");
179 fhNumClustersPhiPrimePtAfterCut->GetAxis(1)->SetTitle("#phi'");
180 fhNumClustersPhiPrimePtAfterCut->GetAxis(2)->SetTitle("Ncl");
181 fOutputContainer->Add(fhNumClustersPhiPrimePtAfterCut);
183 fhPhiPrimeCutEfficiency = new TH1F("hPhiPrimeCutEfficiency", "Efficiency of #phi' cut;p_{T} (GeV/c);Cut efficiency", nPtBins, binsPt);
184 fOutputContainer->Add(fhPhiPrimeCutEfficiency);
187 PostData(1, fOutputContainer);
191 //________________________________________________________________________
192 void AliTPCPIDEtaQA::UserExec(Option_t *)
195 // Called for each event
197 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
199 Printf("ERROR: fESD not available");
203 fMC = dynamic_cast<AliMCEvent*>(MCEvent());
205 if (!fPIDResponse || !fV0KineCuts)
208 if (!GetVertexIsOk(fESD))
211 const AliVVertex* primaryVertex = fESD->GetPrimaryVertexTracks();
215 if(primaryVertex->GetNContributors() <= 0)
218 Double_t magField = fESD->GetMagneticField();
220 // Fill V0 arrays with V0s
223 Int_t multiplicity = fESD->GetNumberOfESDTracks();
225 // Track loop to fill a Train spectrum
226 for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
227 AliESDtrack* track = fESD->GetTrack(iTracks);
229 Printf("ERROR: Could not receive track %d", iTracks);
233 if (TMath::Abs(track->Eta()) >= fEtaCut) continue;
235 //if (TMath::Abs(track->Eta()) < 0.6 || TMath::Abs(track->Eta()) > 0.8) continue;
237 // Do not distinguish positively and negatively charged V0's
238 Char_t v0tagAllCharges = TMath::Abs(GetV0tag(iTracks));
239 if (v0tagAllCharges == -99) {
240 AliError(Form("Problem with V0 tag list (requested V0 track for track %d from %d (list status %d))!", iTracks, fESD->GetNumberOfTracks(),
245 Bool_t isV0el = v0tagAllCharges == 14;
246 Bool_t isV0pi = v0tagAllCharges == 15;
247 Bool_t isV0pr = v0tagAllCharges == 16;
248 Bool_t isV0 = isV0el || isV0pi || isV0pr;
250 // Apply track cut. Accept track, if from V0 or if accepted by cut.
251 // Do not apply the track cuts to V0s, since the track cuts are usually
252 // cuts on primaries, what will throw away almost all V0s.
253 if(!isV0 && (fTrackFilter && !fTrackFilter->IsSelected(track)))
256 if (GetUseTPCCutMIGeo()) {
257 // If cut on geometry is active, don't cut on number of clusters, since such a cut is already included.
259 if (!TPCCutMIGeo(track, InputEvent()))
263 // One should not cut on geometry for V0's since they have different topology. (Loosely) cut on num of clusters instead.
264 if (track->GetTPCsignalN() < 60)
269 // If cut on geometry is not active, always cut on num clusters
270 if (track->GetTPCsignalN() < 60)
274 Double_t pT = track->Pt();
275 //Double_t phiPrime = GetPhiPrime(track->Phi(), magField, track->Charge());
277 //OLD clusterQA Double_t entryQA[3] = { pT, phiPrime, track->GetTPCsignalN() };
278 //OLD clusterQA fhNumClustersPhiPrimePtBeforeCut->Fill(entryQA);
280 // Apply PhiPrimeCut only on high-pt tracks, in this case: Tracks with pt >= fPtThresholdForPhiCut
281 if(fUsePhiCut && pT >= fPtThresholdForPhiCut) {
284 if (!PhiPrimeCut(track, magField))
285 continue; // reject track
287 //Use cut for second branch?
288 if(!PhiPrimeCut(track, magField))
289 //if(!PhiPrimeCut(track, magField)
290 // || (phiPrime < fPhiCutSecondBranchHigh->Eval(pT) && phiPrime > fPhiCutSecondBranchLow->Eval(pT)))
291 continue; // reject track
294 //OLD clusterQA fhNumClustersPhiPrimePtAfterCut->Fill(entryQA);
296 Int_t label = track->GetLabel();
298 AliMCParticle* mcTrack = 0x0;
302 continue; // If MC is available, reject tracks with negative ESD label
303 mcTrack = dynamic_cast<AliMCParticle*>(fMC->GetTrack(TMath::Abs(label)));
305 Printf("ERROR: Could not receive mcTrack with label %d for track %d", label, iTracks);
309 /*// Only accept MC primaries
310 if (!fMC->Stack()->IsPhysicalPrimary(TMath::Abs(label))) {
316 Double_t pTPC = track->GetTPCmomentum();
319 Float_t eta = track->Eta();
322 Double_t dEdxTPC = track->GetTPCsignal();
324 Double_t dEdxExpectedEl = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxDefault,
325 fPIDResponse->UseTPCEtaCorrection(),
326 fPIDResponse->UseTPCMultiplicityCorrection());
327 Double_t dEdxExpectedKa = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kKaon, AliTPCPIDResponse::kdEdxDefault,
328 fPIDResponse->UseTPCEtaCorrection(),
329 fPIDResponse->UseTPCMultiplicityCorrection());
330 Double_t dEdxExpectedPi = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kPion, AliTPCPIDResponse::kdEdxDefault,
331 fPIDResponse->UseTPCEtaCorrection(),
332 fPIDResponse->UseTPCMultiplicityCorrection());
333 Double_t dEdxExpectedPr = fPIDResponse->GetTPCResponse().GetExpectedSignal(track, AliPID::kProton, AliTPCPIDResponse::kdEdxDefault,
334 fPIDResponse->UseTPCEtaCorrection(),
335 fPIDResponse->UseTPCMultiplicityCorrection());
337 Double_t deltaPrimeEl = (dEdxExpectedEl > 0) ? (dEdxTPC / dEdxExpectedEl) : (-1);
338 Double_t deltaPrimeKa = (dEdxExpectedKa > 0) ? (dEdxTPC / dEdxExpectedKa) : (-1);
339 Double_t deltaPrimePi = (dEdxExpectedPi > 0) ? (dEdxTPC / dEdxExpectedPi) : (-1);
340 Double_t deltaPrimePr = (dEdxExpectedPr > 0) ? (dEdxTPC / dEdxExpectedPr) : (-1);
344 Int_t pdg = mcTrack->PdgCode();
346 if (TMath::Abs(pdg) == 211) {//Pion
347 // If V0, only accept if correctly identified
357 else if (TMath::Abs(pdg) == 321) {//Kaon
358 // No kaons from V0 => Reject
364 else if (TMath::Abs(pdg) == 2212) {//Proton
365 // If V0, only accept if correctly identified
375 else if (TMath::Abs(pdg) == 11) {//Electron
376 // If V0, only accept if correctly identified
386 else if (TMath::Abs(pdg) == 13) {//Muon
387 // No muons from V0 => Reject
394 continue; // Reject all other particles
397 Bool_t isKaon = kFALSE;
398 Bool_t isPion = kFALSE;
399 Bool_t isProton = kFALSE;
400 Bool_t isElectron = kFALSE;
401 Bool_t isV0plusTOFel = kFALSE;
405 UInt_t status = track->GetStatus();
406 Bool_t hasTOFout = status&AliESDtrack::kTOFout;
407 Bool_t hasTOFtime = status&AliESDtrack::kTIME;
408 Bool_t hasTOF = kFALSE;
409 if (hasTOFout && hasTOFtime)
411 Float_t length = track->GetIntegratedLength();
412 // Length check only for primaries, not for V0's!
413 if (length < 350 && !isV0)
417 // Select Kaons, Pions and Protons in 3 sigma band (TOF and TPC). Below some momentum threshold, only use TPC cut,
418 // since TOF efficiency is really bad.
419 // In case of ambiguity, add entries for corresponding species
420 if ((pTPC <= 0.25 && TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < 6.0) ||
421 (pTPC > 0.25 && pTPC <= 0.3 && TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < 4.0) ||
422 (pTPC > 0.3 && pTPC <= 0.35 && TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < 3.0) ||
423 (pTPC > 0.35 && TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < 3.0 && hasTOF &&
424 TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon)) < 3.0)) {
427 if ((dEdxTPC >= 50. / TMath::Power(pTPC, 1.3)) && // Pattern recognition instead of TPC cut to be ~independent of old TPC expected dEdx
429 (pTPC >= 0.6 && hasTOF && TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) < 3.0))) {
433 if ((TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) < (pTPC < 0.3 ? 8.0 : 5.0)) &&
435 (pTPC >= 0.6 && hasTOF && TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton)) < 3.0))) {
438 if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) < 3.0 &&
440 (pTPC > 0.3 && hasTOF && TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion)) < 3.0))) {
443 if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron)) < 3.0 &&
445 (pTPC > 0.3 && hasTOF && TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kElectron)) < 3.0))) {
449 else if (!fMC && isV0el) {
450 // Special treatment for V0 electrons -> Look for V0+TOF electrons
452 UInt_t status = track->GetStatus();
453 Bool_t hasTOFout = status&AliESDtrack::kTOFout;
454 Bool_t hasTOFtime = status&AliESDtrack::kTIME;
455 Bool_t hasTOF = kFALSE;
456 if (hasTOFout && hasTOFtime)
459 // No length check for V0's
460 if (hasTOF && TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, AliPID::kElectron)) < 3.0) {
461 isV0plusTOFel = kTRUE;
465 // MC PID, SelectSpecies, P(TPC_inner), multiplicity, deltaPrimeSpecies, eta
466 Double_t entry[6] = { (Double_t)binMC, 0., pTPC, (Double_t)multiplicity, deltaPrimeEl, eta };
469 fhPIDdataAll->Fill(entry);
472 entry[4] = deltaPrimeKa;
473 fhPIDdataAll->Fill(entry);
476 entry[4] = deltaPrimePi;
477 fhPIDdataAll->Fill(entry);
480 entry[4] = deltaPrimePr;
481 fhPIDdataAll->Fill(entry);
484 for (Int_t i = 0; i < 8; i++) {
485 if (i == 0 && isKaon)
487 else if (i == 1 && isPion)
489 else if (i == 2 && isProton)
491 else if (i == 3 && isElectron)
493 else if (i == 4 && isV0plusTOFel)
495 else if (i == 5 && isV0el)
497 else if (i == 6 && isV0pi)
499 else if (i == 7 && isV0pr)
507 entry[4] = deltaPrimeEl;
508 fhPIDdataAll->Fill(entry);
511 entry[4] = deltaPrimeKa;
512 fhPIDdataAll->Fill(entry);
515 entry[4] = deltaPrimePi;
516 fhPIDdataAll->Fill(entry);
519 entry[4] = deltaPrimePr;
520 fhPIDdataAll->Fill(entry);
526 PostData(1, fOutputContainer);
528 // Clear the V0 PID arrays
532 //________________________________________________________________________
533 void AliTPCPIDEtaQA::Terminate(const Option_t *)
535 // Called once at the end of the query
537 TObjArray* output = (TObjArray*)GetOutputData(1);
541 TH1F* hPhiPrimeCutEfficiency = (TH1F*)(output->FindObject("hPhiPrimeCutEfficiency"));
542 if (!hPhiPrimeCutEfficiency)
545 THnSparse* hNumClustersPhiPrimePtBeforeCut = (THnSparse*)(output->FindObject("hNumClustersPhiPrimePtBeforeCut"));
546 if (!hNumClustersPhiPrimePtBeforeCut)
549 THnSparse* hNumClustersPhiPrimePtAfterCut = (THnSparse*)(output->FindObject("hNumClustersPhiPrimePtAfterCut"));
550 if (!hNumClustersPhiPrimePtAfterCut)
553 TH1D* hNumEntriesBefore = hNumClustersPhiPrimePtBeforeCut->Projection(0);
554 TH1D* hNumEntriesAfter = hNumClustersPhiPrimePtAfterCut->Projection(0);
556 for (Int_t bin = 1; bin <= hNumEntriesBefore->GetXaxis()->GetNbins(); bin++) {
557 Double_t numEntriesBefore = hNumEntriesBefore->GetBinContent(bin);
559 if (numEntriesBefore > 0) {
560 Double_t numEntriesAfter = hNumEntriesAfter->GetBinContent(bin);
561 hPhiPrimeCutEfficiency->SetBinContent(bin, numEntriesAfter / numEntriesBefore);
562 // Errors are 100%-correlated: Accepted tracks should be a constant fraction of all tracks in given bin.
563 // Thus: Only take statistical fluctuation from accepted tracks as error
564 hPhiPrimeCutEfficiency->SetBinError(bin, TMath::Sqrt(numEntriesAfter) / numEntriesBefore);
565 // (numEntriesAfter / numEntriesBefore) * TMath::Sqrt(1. / numEntriesAfter + 1. / numEntriesBefore));
569 hPhiPrimeCutEfficiency->SetDrawOption("e");
574 //________________________________________________________________________
575 void AliTPCPIDEtaQA::SetUpHist(THnSparse* hist, Double_t* binsPt) const
577 // Sets bin limits for axes which are not standard binned and the axes titles.
578 // MC PID, SelectSpecies, P(TPC_inner), multiplicity, deltaPrimeSpecies, eta
579 hist->SetBinEdges(2, binsPt);
582 hist->GetAxis(0)->SetTitle("MC PID");
583 hist->GetAxis(0)->SetBinLabel(1, "e");
584 hist->GetAxis(0)->SetBinLabel(2, "K");
585 hist->GetAxis(0)->SetBinLabel(3, "#mu");
586 hist->GetAxis(0)->SetBinLabel(4, "#pi");
587 hist->GetAxis(0)->SetBinLabel(5, "p");
588 hist->GetAxis(0)->SetBinLabel(6, "V0+TOF e");
589 hist->GetAxis(0)->SetBinLabel(7, "V0 e");
590 hist->GetAxis(0)->SetBinLabel(8, "V0 #pi");
591 hist->GetAxis(0)->SetBinLabel(9, "V0 p");
593 hist->GetAxis(1)->SetTitle("Select Species");
594 hist->GetAxis(1)->SetBinLabel(1, "e");
595 hist->GetAxis(1)->SetBinLabel(2, "K");
596 hist->GetAxis(1)->SetBinLabel(3, "#pi");
597 hist->GetAxis(1)->SetBinLabel(4, "p");
599 hist->GetAxis(2)->SetTitle("p_{TPC_inner} (GeV/c)");
601 hist->GetAxis(3)->SetTitle("Event multiplicity");
603 hist->GetAxis(4)->SetTitle("TPC #Delta'{species}");
605 hist->GetAxis(5)->SetTitle("#eta");