2 /**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////
20 // Analysis for identified particle spectra measured with TPC dE/dx. //
23 ///////////////////////////////////////////////////////////////////////////
25 #include "Riostream.h"
34 #include "TObjArray.h"
38 #include "AliAnalysisTaskSE.h"
39 #include "AliAnalysisManager.h"
41 #include "AliHeader.h"
42 #include "AliGenPythiaEventHeader.h"
43 #include "AliGenCocktailEventHeader.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliESDVertex.h"
48 #include "AliESDEvent.h"
49 #include "AliESDInputHandler.h"
50 #include "AliESDtrack.h"
51 #include "AliESDpid.h"
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
59 #include "AliAnalysisTaskChargedHadronSpectra.h"
62 ClassImp(AliAnalysisTaskChargedHadronSpectra)
64 //________________________________________________________________________
65 AliAnalysisTaskChargedHadronSpectra::AliAnalysisTaskChargedHadronSpectra()
66 : AliAnalysisTaskSE("TaskChargedHadron"), fESD(0), fListHist(0), fESDtrackCuts(0),fESDpid(0),
75 fHistPtEtaKaonNoKink(0),
77 fHistPtEtaProtonDCA(0),
79 fHistPtEtaElectron(0),
81 fHistClassicProton(0),
83 fHistClassicElectron(0),
85 fHistTrackPerEvent(0),
86 fHistTrackPerEventMC(0),
98 // default Constructor
102 //________________________________________________________________________
103 AliAnalysisTaskChargedHadronSpectra::AliAnalysisTaskChargedHadronSpectra(const char *name)
104 : AliAnalysisTaskSE(name), fESD(0), fListHist(0), fESDtrackCuts(0),fESDpid(0),
110 fHistPtMCElectron(0),
113 fHistPtEtaKaonNoKink(0),
115 fHistPtEtaProtonDCA(0),
117 fHistPtEtaElectron(0),
119 fHistClassicProton(0),
121 fHistClassicElectron(0),
123 fHistTrackPerEvent(0),
124 fHistTrackPerEventMC(0),
130 fHistEffProtonDCA(0),
137 // standard constructur which should be used - PID objects is initialized
140 fMCtrue = kTRUE; // change two things for processing real data: 1. set this to false! / 2. change ALEPH parameters
141 fAlephParameters[0] = 4.23232575531564326e+00;//50*0.76176e-1;
142 fAlephParameters[1] = 8.68482806165147636e+00;//10.632;
143 fAlephParameters[2] = 1.34000000000000005e-05;//0.13279e-4;
144 fAlephParameters[3] = 2.30445734159456084e+00;//1.8631;
145 fAlephParameters[4] = 2.25624744086878559e+00;//1.9479;
147 fESDpid = new AliESDpid();
148 fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0]/50.,fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
151 Printf("*** CONSTRUCTOR CALLED ****");
152 // Define input and output slots here
153 // Input slot #0 works with a TChain
154 //DefineInput(0, TChain::Class()); <-> not needed in AliAnalysisTaskSE
156 // Output slot #0 writes into a TList container
157 DefineOutput(1, TList::Class());
164 //________________________________________________________________________
165 void AliAnalysisTaskChargedHadronSpectra::UserCreateOutputObjects()
169 fListHist = new TList();
170 //fListHist->SetOwner(); // Whoever knows how the data handling is ...?
172 const Int_t kPtBins = 2*56;
173 const Double_t kPtMax = 16.0;
174 const Int_t kEtaBins = 4;
175 const Double_t kEtaMax = 0.8;
176 const Int_t kDeDxBins = 200;
177 const Double_t kDeDxMax = 1;
178 const Int_t kMultBins = 10;
179 const Int_t kMultMax = 300;
182 Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
183 Int_t indexes[kPtBins+1];
184 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
185 Double_t binsPt[kPtBins+1];
186 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
190 fHistPtMCKaon = new TH3F("HistPtMCKaon", "PtEtaKaon; mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
191 fHistPtMCProton = new TH3F("HistPtMCProton", "PtEtaProton;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
192 fHistPtMCPion = new TH3F("HistPtMCPion", "PtEtaPion;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
193 fHistPtMCElectron = new TH3F("HistPtMCElectron", "PtEtaElectron;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
194 fHistPtMCMuon = new TH3F("HistPtMCMuon", "PtEtaMuon;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
196 fHistPtMCKaon->GetZaxis()->Set(kPtBins, binsPt);
197 fHistPtMCProton->GetZaxis()->Set(kPtBins, binsPt);
198 fHistPtMCPion->GetZaxis()->Set(kPtBins, binsPt);
199 fHistPtMCElectron->GetZaxis()->Set(kPtBins, binsPt);
200 fHistPtMCMuon->GetZaxis()->Set(kPtBins, binsPt);
202 fListHist->Add(fHistPtMCKaon);
203 fListHist->Add(fHistPtMCProton);
204 fListHist->Add(fHistPtMCPion);
205 fListHist->Add(fHistPtMCElectron);
206 fListHist->Add(fHistPtMCMuon);
208 // reconstructed particle histograms
209 fHistPtEtaKaon = new TH3F("HistPtEtaKaon", "PtEtaKaon;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
210 fHistPtEtaKaonNoKink = new TH3F("HistPtEtaKaonNoKink", "PtEtaKaonNoKink;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
211 fHistPtEtaProton = new TH3F("HistPtEtaProton", "PtEtaProton;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
212 fHistPtEtaProtonDCA = new TH3F("HistPtEtaProtonDCA", "PtEtaProton;DCA (cm); #eta; p_{T} (GeV)",200,0,15,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
213 fHistPtEtaPion = new TH3F("HistPtEtaPion", "PtEtaPion;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
214 fHistPtEtaElectron = new TH3F("HistPtEtaElectron", "PtEtaElectron;mult; #eta; p_{T} (GeV)",kMultBins,-0.5,kMultMax,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
216 fHistPtEtaKaon->GetZaxis()->Set(kPtBins, binsPt);
217 fHistPtEtaKaonNoKink->GetZaxis()->Set(kPtBins, binsPt);
218 fHistPtEtaProton->GetZaxis()->Set(kPtBins, binsPt);
219 fHistPtEtaProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
220 fHistPtEtaPion->GetZaxis()->Set(kPtBins, binsPt);
221 fHistPtEtaElectron->GetZaxis()->Set(kPtBins, binsPt);
223 fListHist->Add(fHistPtEtaKaon);
224 fListHist->Add(fHistPtEtaKaonNoKink);
225 fListHist->Add(fHistPtEtaProton);
226 fListHist->Add(fHistPtEtaProtonDCA);
227 fListHist->Add(fHistPtEtaPion);
228 fListHist->Add(fHistPtEtaElectron);
230 // histograms for the classical analysis
231 fHistClassicKaon = new TH3F("HistClassicKaon", "PtEtaKaon;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
232 fHistClassicProton = new TH3F("HistClassicProton", "PtEtaProton;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
233 fHistClassicPion = new TH3F("HistClassicPion", "PtEtaPion;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
234 fHistClassicElectron = new TH3F("HistClassicElectron", "PtEtaElectron;p_{T} (GeV);#eta;delta dEdx",kPtBins,-kPtMax,kPtMax,kEtaBins,0,kEtaMax,kDeDxBins,-kDeDxMax,kDeDxMax);
236 fHistClassicKaon->GetXaxis()->Set(kPtBins, binsPt);
237 fHistClassicProton->GetXaxis()->Set(kPtBins, binsPt);
238 fHistClassicPion->GetXaxis()->Set(kPtBins, binsPt);
239 fHistClassicElectron->GetXaxis()->Set(kPtBins, binsPt);
241 fListHist->Add(fHistClassicKaon);
242 fListHist->Add(fHistClassicProton);
243 fListHist->Add(fHistClassicPion);
244 fListHist->Add(fHistClassicElectron);
246 // histograms of general interest
247 fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,20.,1000,0.,500);
249 fListHist->Add(fDeDx);
250 fHistTrackPerEvent = new TH1F("HistTrackPerEvent", "Tracks per event;Number of Tracks;Number of Events",101,-0.5,100);
251 fListHist->Add(fHistTrackPerEvent);
252 fHistTrackPerEventMC = new TH2F("HistTrackPerEventMC", "Tracks per event;Number of Tracks;Number of Events",10,-0.5,9,101,-0.5,100);
253 fListHist->Add(fHistTrackPerEventMC);
254 fSecProtons = new TH2F("SecProtons", "xy;x;y",100,-10,10,100,-10,10);
255 fListHist->Add(fSecProtons);
256 fVertexZ = new TH1F("VertexZ", "vertex position;z (cm);counts",200,-50,50);
257 fListHist->Add(fVertexZ);
259 fHistEtaNcls = new TH2F("HistEtaNcls", "EtaNcls;#eta;Ncls",100,-1.5,1.5,80,0,160);
260 fListHist->Add(fHistEtaNcls);
261 fHistEtaPhi = new TH2F("HistEtaPhi", "EtaNcls;#eta;#phi",100,-1.5,1.5,100,0,2*TMath::Pi());
262 fListHist->Add(fHistEtaPhi);
264 // histograms for a refined efficiency study
265 fHistEffProton = new TH3F("HistEffProton", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
266 fListHist->Add(fHistEffProton);
267 fHistEffProtonDCA = new TH3F("HistEffProtonDCA", "Eff;p_{T} (GeV); code",10,-0.5,9.5,200,0,15,kPtBins,-kPtMax,kPtMax);
268 fListHist->Add(fHistEffProtonDCA);
269 fHistEffPion = new TH3F("HistEffPion", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
270 fListHist->Add(fHistEffPion);
271 fHistEffKaon = new TH3F("HistEffKaon", "Eff;p_{T} (GeV); code",10,-0.5,9.5,kEtaBins,0,kEtaMax,kPtBins,-kPtMax,kPtMax);
272 fListHist->Add(fHistEffKaon);
274 fHistEffProton->GetZaxis()->Set(kPtBins, binsPt);
275 fHistEffProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
276 fHistEffPion->GetZaxis()->Set(kPtBins, binsPt);
277 fHistEffKaon->GetZaxis()->Set(kPtBins, binsPt);
281 // histograms for dN/dpT
282 fHighPtElectrons = new TH1F("HistHighPtElectrons", "backgr;p_{T} (GeV); electron contamination (%)",15,2,15);
283 fListHist->Add(fHighPtElectrons);
284 fHighPtHadrons = new TH1F("HistHighPtHadrons", "backgr;p_{T} (GeV); electron contamination (%)",15,2,15);
285 fListHist->Add(fHighPtHadrons);
288 // Create Objects which are needed for the rest of the analysis
290 fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
291 fESDtrackCuts->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2); // BEWARE STANDARD VALUES ARE: 2, 2, 0.5, 0.5, 2
292 fESDtrackCuts->SetMaxNsigmaToVertex(3);
293 fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
294 fESDtrackCuts->SetRequireTPCRefit(kFALSE);
295 fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
296 fESDtrackCuts->SetMinNClustersTPC(100);
297 fESDtrackCuts->SetMaxChi2PerClusterTPC(3.5);
298 //fESDtrackCuts->SetMaxDCAToVertex(3);
300 fListHist->Add(fESDtrackCuts);
306 //________________________________________________________________________
307 void AliAnalysisTaskChargedHadronSpectra::UserExec(Option_t *)
311 AliLog::SetGlobalLogLevel(AliLog::kError);
313 // Check Monte Carlo information and other access first:
315 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
317 //Printf("ERROR: Could not retrieve MC event handler");
321 AliMCEvent* mcEvent = eventHandler->MCEvent();
323 //Printf("ERROR: Could not retrieve MC event");
327 fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
329 //Printf("ERROR: fESD not available");
333 if (!fESDtrackCuts) {
334 Printf("ERROR: fESDtrackCuts not available");
338 Int_t trackCounter = 0;
340 // monitor vertex position
341 const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
343 fHistTrackPerEvent->Fill(trackCounter);
344 PostData(1, fListHist);
347 if (vertex->GetZv() != 0) fVertexZ->Fill(vertex->GetZv());
348 if (TMath::Abs(vertex->GetZv()) > 10) {
349 fHistTrackPerEvent->Fill(trackCounter);
350 PostData(1, fListHist);
357 // 1st track loop to determine multiplicities
358 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
359 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
360 if (!track) continue;
361 if (fESDtrackCuts->AcceptTrack(track) && TMath::Abs(track->Eta())< 0.9) {
370 AliStack* stack = mcEvent->Stack();
376 AliHeader * header = mcEvent->Header();
377 Int_t processtype = GetPythiaEventProcessType(header);
379 if (processtype !=92 && processtype !=93 && processtype != 94) fHistTrackPerEventMC->Fill(1., trackCounter);
380 // single diffractive
381 if ((processtype == 92 || processtype == 93)) fHistTrackPerEventMC->Fill(2., trackCounter);
382 // double diffractive
383 if (processtype == 94) fHistTrackPerEventMC->Fill(3., trackCounter);
385 Int_t mult = trackCounter;//stack->GetNprimary();
386 fHistTrackPerEventMC->Fill(0., trackCounter);
388 for(Int_t i = 0; i < stack->GetNtrack(); i++) {
389 TParticle * trackMC = stack->Particle(i);
390 Int_t pdg = trackMC->GetPdgCode();
392 Double_t xv = trackMC->Vx();
393 Double_t yv = trackMC->Vy();
394 Double_t zv = trackMC->Vz();
396 dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
398 dz = TMath::Abs(zv); // so stupid to avoid warnings
400 // vertex cut - selection of primaries
402 //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
403 if (pdg == 11) fHistPtMCElectron->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select e- only; before vertex
404 if (pdg == -11) fHistPtMCElectron->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select e+ only; before vertex
406 if (!stack->IsPhysicalPrimary(i)) continue;
408 if (pdg == 321) fHistPtMCKaon->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select K+ only
409 if (pdg == 211) fHistPtMCPion->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select Pi+ only
410 if (pdg == 2212) fHistPtMCProton->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select p+ only
411 if (pdg == -13) fHistPtMCMuon->Fill(mult, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // select mu+ only
413 if (pdg == -321) fHistPtMCKaon->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select K- only
414 if (pdg == -211) fHistPtMCPion->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select Pi- only
415 if (pdg == -2212) fHistPtMCProton->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select p- only
416 if (pdg == 13) fHistPtMCMuon->Fill(mult, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // select mu- only
418 // special Kaon checks - those which decay before entering the TPC fiducial volume
420 if (TMath::Abs(pdg)==321) {
421 TParticle * trackDaughterMC = stack->Particle(TMath::Abs(trackMC->GetFirstDaughter()));
422 Int_t pdgDaughter = trackDaughterMC->GetPdgCode();
423 if (pdgDaughter == TMath::Abs(13) && pdg == 321) fHistEffKaon->Fill(5, TMath::Abs(trackMC->Eta()), trackMC->Pt()); // Kaon control hist
424 if (pdgDaughter == TMath::Abs(13) && pdg == -321) fHistEffKaon->Fill(5, TMath::Abs(trackMC->Eta()), -trackMC->Pt()); // Kaon control hist
433 // definition of PID functions --> to be put to CreateOutputObjects()
435 TF1 foProton("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
436 TF1 foPion("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
437 TF1 foElec("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
438 TF1 foKaon("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
439 TF1 foMuon("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
441 foProton.SetParameters(fAlephParameters);
442 foPion.SetParameters(fAlephParameters);
443 foElec.SetParameters(fAlephParameters);
444 foKaon.SetParameters(fAlephParameters);
445 foMuon.SetParameters(fAlephParameters);
447 const Float_t k2sigmaCorr = 1/0.9545;
449 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
451 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
453 // 1. select TPC tracks only and propagate them to the primary vertex determined with the TPC standalone
455 AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
456 if (!track) continue;
457 AliESDtrack *trackDummy = fESD->GetTrack(i);
458 if (!trackDummy->GetInnerParam()) { // DEBUG DEBUG
462 AliExternalTrackParam trackIn(*trackDummy->GetInnerParam());
463 Double_t ptot = trackIn.GetP(); // momentum for dEdx determination
464 Double_t pT = track->Pt();
465 Double_t eta = TMath::Abs(track->Eta());
467 fHistEtaNcls->Fill(track->Eta(),track->GetTPCNcls());
468 fHistEtaPhi->Fill(track->Eta(),track->Phi());
469 // cut for dead regions in the detector
470 // if (track->Eta() > 0.1 && (track->Eta() < 0.2 && track->Phi() > 0.1 && track->Phi() < 0.1) continue;
472 // 2.a) apply some standard track cuts according to general recommendations
475 if (!fESDtrackCuts->AcceptTrack(track)) {
478 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
479 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
480 if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(6,eta,track->GetSign()*pT); // Kaon control hist
490 Double_t sign = track->GetSign();
492 Double_t tpcSignal = track->GetTPCsignal();
494 fDeDx->Fill(ptot, tpcSignal);
498 fHistClassicKaon->Fill(sign*pT,eta,(tpcSignal-foKaon.Eval(ptot))/foKaon.Eval(ptot));
499 fHistClassicProton->Fill(sign*pT,eta,(tpcSignal-foProton.Eval(ptot))/foProton.Eval(ptot));
500 fHistClassicPion->Fill(sign*pT,eta,(tpcSignal-foPion.Eval(ptot))/foPion.Eval(ptot));
501 fHistClassicElectron->Fill(sign*pT,eta,(tpcSignal-foElec.Eval(ptot))/foElec.Eval(ptot));
503 // fill histograms for dNdpT
506 fHighPtHadrons->Fill(pT);
507 Double_t delta = (tpcSignal - foElec.Eval(ptot))/foElec.Eval(ptot);
508 if (delta > 0) fHighPtElectrons->Fill(pT);
511 /* 2sigma PID with 2sigma eff correction! */
512 Double_t tpcMom = track->GetP();
513 const AliExternalTrackParam *in = track->GetTPCInnerParam();
517 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kPion)) < 2) {
518 fHistPtEtaPion->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
520 if (trackCounter < 300 && fMCtrue) {
521 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
522 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
523 if (pdg == 211 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffPion->Fill(0,eta,sign*pT,k2sigmaCorr);
524 if (pdg == 211 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
525 fHistEffPion->Fill(1,eta,sign*pT,k2sigmaCorr);
526 TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
527 Int_t code = TMath::Abs(trackMother->GetPdgCode());
528 if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) fHistEffPion->Fill(3,eta,sign*pT,k2sigmaCorr);
530 if (TMath::Abs(trackMC->GetPdgCode()) != 211) fHistEffPion->Fill(2,eta,sign*pT,k2sigmaCorr);
531 if (TMath::Abs(trackMC->GetPdgCode()) == 13) fHistEffPion->Fill(4,eta,sign*pT,k2sigmaCorr);
535 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon)) < 2) {
536 fHistPtEtaKaon->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
538 if (trackCounter < 300 && fMCtrue) {
539 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
540 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
541 if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(0,eta,sign*pT,k2sigmaCorr);
542 if (pdg == 321 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
543 fHistEffKaon->Fill(1,eta,sign*pT,k2sigmaCorr);
544 TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
545 Int_t code = TMath::Abs(trackMother->GetPdgCode());
546 if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) fHistEffKaon->Fill(3,eta,sign*pT,k2sigmaCorr);
548 if (TMath::Abs(trackMC->GetPdgCode()) != 321) fHistEffKaon->Fill(2,eta,sign*pT,k2sigmaCorr);
552 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon)) < 2 && track->GetKinkIndex(0)==0) fHistPtEtaKaonNoKink->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
554 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kProton)) < 2) {
555 fHistPtEtaProton->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
557 track->GetImpactParameters(dca,cov);
558 Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0] + dca[1]*dca[1]);
559 fHistPtEtaProtonDCA->Fill(primVtxDCA, eta, sign*pT, k2sigmaCorr);
561 if (trackCounter < 300 && fMCtrue) {
562 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
563 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
564 if (pdg == 2212 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
565 fHistEffProton->Fill(0.,eta,sign*pT,k2sigmaCorr);
566 if (eta < 0.15) fHistEffProtonDCA->Fill(0.,primVtxDCA,sign*pT,k2sigmaCorr);
568 if (pdg == 2212 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
569 fHistEffProton->Fill(1,eta,sign*pT,k2sigmaCorr);
570 if (eta < 0.15) fHistEffProtonDCA->Fill(1,primVtxDCA,sign*pT,k2sigmaCorr);
571 TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
572 Int_t code = TMath::Abs(trackMother->GetPdgCode());
573 if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) {
574 fHistEffProton->Fill(3,eta,sign*pT,k2sigmaCorr);
575 if (eta < 0.15) fHistEffProtonDCA->Fill(3,primVtxDCA,sign*pT,k2sigmaCorr);
577 if (code!=3122&&code!=3222&&code!=3212&&code!=3112&&code!=3322&&code!=3312&&code!=3332&&code!=130&&code!=310) fSecProtons->Fill(trackMC->Vx(),trackMC->Vy());
579 if (TMath::Abs(trackMC->GetPdgCode()) != 2212) fHistEffProton->Fill(2,eta,sign*pT,k2sigmaCorr);
583 if (TMath::Abs(fESDpid->NumberOfSigmasTPC(track,AliPID::kElectron))) fHistPtEtaElectron->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
587 } // end of track loop
589 fHistTrackPerEvent->Fill(trackCounter);
594 PostData(1, fListHist);
599 //________________________________________________________________________
600 void AliAnalysisTaskChargedHadronSpectra::Terminate(Option_t *)
602 // Draw result to the screen
603 // Called once at the end of the query
604 cout << "*** TERMINATE ***" << endl;
611 //________________________________________________________________________
612 void AliAnalysisTaskChargedHadronSpectra::BinLogX(const TH1 *h) {
614 // Method for the correct logarithmic binning of histograms
616 TAxis *axis = h->GetXaxis();
617 int bins = axis->GetNbins();
619 Double_t from = axis->GetXmin();
620 Double_t to = axis->GetXmax();
621 Double_t *newBins = new Double_t[bins + 1];
624 Double_t factor = pow(to/from, 1./bins);
626 for (int i = 1; i <= bins; i++) {
627 newBins[i] = factor * newBins[i-1];
629 axis->Set(bins, newBins);
635 //________________________________________________________________________
636 Int_t AliAnalysisTaskChargedHadronSpectra::GetPythiaEventProcessType(const AliHeader* aHeader, const Bool_t adebug) const {
638 // get the process type of the event.
641 // can only read pythia headers, either directly or from cocktalil header
642 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
644 if (!pythiaGenHeader) {
646 AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
647 if (!genCocktailHeader) {
648 //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
652 TList* headerList = genCocktailHeader->GetHeaders();
657 for (Int_t i=0; i<headerList->GetEntries(); i++) {
658 pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
663 if (!pythiaGenHeader) {
664 //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Could not find Pythia header. \n");
670 //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
673 return pythiaGenHeader->ProcessType();
677 //________________________________________________________________________
678 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicProton(const TH3F * input, Int_t EtaBin, const Double_t * AlephParams) {
680 // The multiple Gauss fitting is happening here...
681 // input histogram is of the form (Pt, eta, difference in dEdx)
684 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
685 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
686 TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
687 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
688 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
690 foProton->SetParameters(AlephParams);
691 foPion->SetParameters(AlephParams);
692 foElec->SetParameters(AlephParams);
693 foKaon->SetParameters(AlephParams);
694 foMuon->SetParameters(AlephParams);
697 const Double_t kSigma = 0.06; // expected resolution (roughly)
699 const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
700 const Double_t kPtMax = input->GetXaxis()->GetXmax();
701 const Int_t kEtaBins = input->GetYaxis()->GetNbins();
703 TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
705 Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
706 Int_t indexes[kPtBins+1];
707 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
708 Double_t binsPt[kPtBins+1];
709 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
710 histPt->GetXaxis()->Set(kPtBins, binsPt);
713 TCanvas * canvMany = new TCanvas("canvManyProton","canvManyProton");
714 canvMany->Print("canvManyProton.ps[");
716 for(Int_t x=1; x < kPtBins+1; x++) {
717 Double_t pT = input->GetXaxis()->GetBinCenter(x);
718 for(Int_t y=1; y < kEtaBins+1; y++) {
719 Double_t eta = input->GetYaxis()->GetBinCenter(y);
720 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
721 Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
722 TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
723 histDeDx->SetTitle(Form("pt_%f",pT));
724 Float_t maxYield = histDeDx->GetEntries()/(TMath::Sqrt(2*TMath::Pi())*0.04/histDeDx->GetBinWidth(1));
726 TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-1/*-7*kSigma*/,1/*7*kSigma*/);
727 Double_t parameters[12] = {maxYield/2.,0,kSigma,maxYield/2.,(foPion->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma,maxYield/2.,(foElec->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma,maxYield/2.,(foKaon->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot),kSigma};
728 gausFit->SetParameters(parameters);
729 gausFit->SetParLimits(0,0.,maxYield);
730 gausFit->SetParLimits(1,-0.05,0.05);
731 gausFit->SetParLimits(2,0.04,0.08);
733 Double_t pionPosition = (foPion->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
734 gausFit->SetParLimits(0,0.,maxYield);
735 gausFit->SetParLimits(4, 0.95*pionPosition, 1.05*pionPosition);
736 gausFit->SetParLimits(5,0.8*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)));
738 Double_t elecPosition = (foElec->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
739 gausFit->SetParLimits(0,0.,maxYield);
740 gausFit->SetParLimits(7, 0.95*elecPosition, 1.05*elecPosition);
741 gausFit->SetParLimits(8,0.8*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)));
743 Double_t kaonPosition = (foKaon->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
744 gausFit->SetParLimits(0,0.,maxYield);
745 gausFit->SetParLimits(10, 0.95*kaonPosition, 1.05*kaonPosition);
746 gausFit->SetParLimits(11,0.8*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)));
747 histDeDx->Fit("gausFit","QNR");
748 gausFit->GetParameters(parameters);
753 histDeDx->SetMarkerStyle(21);
754 histDeDx->SetMarkerSize(0.7);
756 gausFit->SetLineWidth(2);
757 gausFit->Draw("same");
759 TF1 gausFit0("gausFit0", "gaus(0)",-1,1);
760 TF1 gausFit1("gausFit1", "gaus(0)",-1,1);
761 TF1 gausFit2("gausFit2", "gaus(0)",-1,1);
762 TF1 gausFit3("gausFit3", "gaus(0)",-1,1);
763 gausFit0.SetParameters(parameters[0],parameters[1],parameters[2]);
764 gausFit1.SetParameters(parameters[3],parameters[4],parameters[5]);
765 gausFit2.SetParameters(parameters[6],parameters[7],parameters[8]);
766 gausFit3.SetParameters(parameters[9],parameters[10],parameters[11]);
768 gausFit0.SetLineColor(4);
769 gausFit1.SetLineColor(2);
770 gausFit2.SetLineColor(6);
771 gausFit3.SetLineColor(8);
773 gausFit0.SetLineWidth(1);
774 gausFit1.SetLineWidth(1);
775 gausFit2.SetLineWidth(1);
776 gausFit3.SetLineWidth(1);
778 gausFit0.Draw("same");
779 gausFit1.Draw("same");
780 gausFit2.Draw("same");
781 gausFit3.Draw("same");
782 canvMany->Print("canvManyProton.ps");
785 Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2)/histDeDx->GetBinWidth(1); // area of the gaus fit
788 // stupid solution --> remove as soon as possible
790 TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
791 gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
792 for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
793 Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
794 yield += gausYield->Eval(delta);
797 if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
801 canvMany->Print("canvManyProton.ps]");
803 TCanvas * canvPt = new TCanvas();
812 //________________________________________________________________________
813 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicPion(const TH3F * input, Int_t EtaBin,const Double_t * AlephParams) {
815 // The multiple Gauss fitting is happening here...
818 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
819 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
820 TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
821 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
822 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
824 foProton->SetParameters(AlephParams);
825 foPion->SetParameters(AlephParams);
826 foElec->SetParameters(AlephParams);
827 foKaon->SetParameters(AlephParams);
828 foMuon->SetParameters(AlephParams);
830 // input histogram is of the form (Pt, eta, difference in dEdx)
832 const Double_t kSigma = 0.06; // expected resolution (roughly)
834 const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
835 const Double_t kPtMax = input->GetXaxis()->GetXmax();
836 const Int_t kEtaBins = input->GetYaxis()->GetNbins();
838 TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
840 Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
841 Int_t indexes[kPtBins+1];
842 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
843 Double_t binsPt[kPtBins+1];
844 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
845 histPt->GetXaxis()->Set(kPtBins, binsPt);
848 TCanvas * canvMany = new TCanvas();
849 canvMany->Divide(8,5);
851 for(Int_t x=1; x < kPtBins+1; x++) {
852 Double_t pT = input->GetXaxis()->GetBinCenter(x);
853 for(Int_t y=1; y < kEtaBins+1; y++) {
854 Double_t eta = input->GetYaxis()->GetBinCenter(y);
855 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
856 Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
857 TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
859 TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-5*kSigma,5*kSigma);
860 //Double_t parameters[12] = {100,0,kSigma,100,TMath::Log(foProton->Eval(ptot)/foPion->Eval(ptot)),kSigma,100,TMath::Log(foElec->Eval(ptot)/foPion->Eval(ptot)),kSigma,100,TMath::Log(foKaon->Eval(ptot)/foPion->Eval(ptot)),kSigma};
861 Double_t parameters[12] = {100,0,kSigma,100,(foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma,100,(foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma,100,(foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot),kSigma};
862 gausFit->SetParameters(parameters);
863 gausFit->SetParLimits(1,-0.05,0.05);//gausFit->FixParameter(1,0); DEBUG DEBUG DEBUG
864 gausFit->SetParLimits(2,0.04,0.08);//gausFit->FixParameter(2,kSigma); DEBUG DEBUG DEBUG
866 gausFit->FixParameter(4,(foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
867 gausFit->FixParameter(5,kSigma);
868 gausFit->FixParameter(7,(foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
869 gausFit->FixParameter(8,kSigma);
870 gausFit->FixParameter(10,(foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
871 gausFit->FixParameter(11,kSigma);
872 histDeDx->Fit("gausFit","Q0");
873 if (y == EtaBin && pT < 0) {
877 gausFit->Draw("same");
879 Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2); // area of the gaus fit
880 // stupid solution --> remove as soon as possible
882 TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
883 gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
884 for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
885 Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
886 yield += gausYield->Eval(delta);
889 if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
895 TCanvas * canvPt = new TCanvas();
904 //________________________________________________________________________
905 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicKaon(const TH3F * input, Int_t EtaBin,const Double_t * AlephParams) {
907 // The multiple Gauss fitting is happening here...
910 TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20);
911 TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
912 TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
913 TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
914 TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
916 foProton->SetParameters(AlephParams);
917 foPion->SetParameters(AlephParams);
918 foElec->SetParameters(AlephParams);
919 foKaon->SetParameters(AlephParams);
920 foMuon->SetParameters(AlephParams);
922 const Double_t kSigma = 0.06; // expected resolution (roughly)
924 const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
925 const Double_t kPtMax = input->GetXaxis()->GetXmax();
926 const Int_t kEtaBins = input->GetYaxis()->GetNbins();
928 TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
930 Double_t binsPtDummy[kPtBins+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, -0.05, -0.1, -0.15, -0.2, -0.25, -0.3, -0.35, -0.4, -0.45, -0.5, -0.55, -0.6, -0.65, -0.7, -0.75, -0.8, -0.85, -0.9, -0.95, -1.0, -1.1, -1.2, -1.3, -1.4, -1.5, -1.6, -1.7, -1.8, -1.9, -2.0, -2.2, -2.4, -2.6, -2.8, -3.0, -3.2, -3.4, -3.6, -3.8, -4.0, -4.5, -5.0, -5.5, -6.0, -6.5, -7.0, -7.5, -8.0, -9.0, -10.0, -11.0, -12.0, -13.0, -14.0, -15.0, -16.0};
931 Int_t indexes[kPtBins+1];
932 TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
933 Double_t binsPt[kPtBins+1];
934 for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
935 histPt->GetXaxis()->Set(kPtBins, binsPt);
938 TCanvas * canvMany = new TCanvas();
939 canvMany->Divide(8,5);
941 for(Int_t x=1; x < kPtBins+1; x++) {
942 Double_t pT = input->GetXaxis()->GetBinCenter(x);
943 for(Int_t y=1; y < kEtaBins+1; y++) {
944 Double_t eta = input->GetYaxis()->GetBinCenter(y);
945 Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
946 Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
947 TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
949 TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-5*kSigma,5*kSigma);
950 //Double_t parameters[12] = {100,0,kSigma,100,TMath::Log(foProton->Eval(ptot)/foKaon->Eval(ptot)),kSigma,100,TMath::Log(foElec->Eval(ptot)/foKaon->Eval(ptot)),kSigma,100,TMath::Log(foPion->Eval(ptot)/foPion->Eval(ptot)),kSigma};
951 Double_t parameters[12] = {100,0,kSigma,100,(foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma,100,(foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot),kSigma,100,(foPion->Eval(ptot)-foKaon->Eval(ptot))/foPion->Eval(ptot),kSigma};
952 gausFit->SetParameters(parameters);
953 gausFit->SetParLimits(1,-0.05,0.05);//gausFit->FixParameter(1,0); DEBUG DEBUG DEBUG
954 gausFit->SetParLimits(2,0.04,0.08);//gausFit->FixParameter(2,kSigma); DEBUG DEBUG DEBUG
956 gausFit->FixParameter(4,(foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot));
957 gausFit->FixParameter(5,kSigma);
958 gausFit->FixParameter(7,(foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot));
959 gausFit->FixParameter(8,kSigma);
960 gausFit->FixParameter(10,(foPion->Eval(ptot)- foKaon->Eval(ptot))/foKaon->Eval(ptot));
961 gausFit->FixParameter(11,kSigma);
962 histDeDx->Fit("gausFit","Q0");
963 if (y == EtaBin && pT < 0) {
967 gausFit->Draw("same");
969 Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2); // area of the gaus fit
970 // stupid solution --> remove as soon as possible
972 TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
973 gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
974 for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
975 Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
976 yield += gausYield->Eval(delta);
979 if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
985 TCanvas * canvPt = new TCanvas();
994 //________________________________________________________________________
995 void AliAnalysisTaskChargedHadronSpectra::Postprocess(const TList * ListOfHistogramsMC,const TList * ListOfHistogramsData, const Char_t *filename) {
998 const Int_t kMultLow = 1;
999 const Int_t kMultHigh = 10;
1001 const Int_t kEtaHigh = 4;
1003 // Extract histograms for the MC efficiency study
1004 TH3F* histPtMCKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCKaon");
1005 TH3F* histPtMCProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCProton");
1006 TH3F* histPtMCPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCPion");
1008 TH3F* histPtEtaKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaKaon");
1009 TH3F* histPtEtaProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaProton");
1010 TH3F* histPtEtaPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaPion");
1012 TH3F* histEffProton = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffProton");
1013 TH3F* histEffKaon = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffKaon");
1014 TH3F* histEffPion = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffPion");
1016 // Extract data for the final spectra
1017 TH3F* histPtEtaKaon = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaKaon");
1018 TH3F* histPtEtaProton = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaProton");
1019 TH3F* histPtEtaPion = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaPion");
1021 // "MC" of the real data for consistency checks
1022 TH3F* histPtMCKaonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCKaon");
1023 TH3F* histPtMCProtonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCProton");
1024 TH3F* histPtMCPionData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCPion");
1026 TFile spectraFile(filename,"RECREATE");
1029 for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1030 TH1D *protonRaw = histPtEtaProtonEff->ProjectionZ(Form("ProtonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1032 TH1D *protonMC = histPtMCProtonEff->ProjectionZ(Form("McProton_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1034 TH1D *protonRawPrim = histEffProton->ProjectionZ(Form("ProtonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1035 protonRawPrim->Sumw2();
1036 TH1D *protonSecReac = histEffProton->ProjectionZ(Form("ProtonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1037 protonSecReac->Sumw2();
1038 TH1D *protonSecWeak = histEffProton->ProjectionZ(Form("ProtonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1039 protonSecWeak->Sumw2();
1040 TH1D *protonMis = histEffProton->ProjectionZ(Form("ProtonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1042 protonSecReac->Add(protonSecWeak,-1);
1044 // 1. total (physical) efficiency
1046 TH1D *protonEffTot = new TH1D(*protonRaw);
1047 protonEffTot->Divide(protonRaw,protonMC);
1048 protonEffTot->SetName(Form("ProtonEffTot_%i",iEta));
1049 protonEffTot->Write();
1051 // 2. contributions from secondaries created by interactions
1053 TH1D *protonEffSecReac = new TH1D(*protonRaw);
1054 TH1D *dummy = new TH1D(*protonRaw);
1055 dummy->Add(protonSecReac,-1);
1056 protonEffSecReac->Divide(protonRaw,dummy);
1057 protonEffSecReac->SetName(Form("ProtonEffSecReac_%i",iEta));
1058 protonEffSecReac->Write();
1061 // 3. contributions from secondaries from weak decays
1063 TH1D *protonEffSecWeak = new TH1D(*protonRaw);
1064 TH1D *dummy2 = new TH1D(*protonRaw);
1065 dummy2->Add(protonSecWeak,-1);
1066 protonEffSecWeak->Divide(protonRaw,dummy2);
1067 protonEffSecWeak->SetName(Form("ProtonEffSecWeak_%i",iEta));
1068 protonEffSecWeak->Write();
1071 // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1073 TH1D *protonEffAbs = new TH1D(*protonRaw);
1074 protonEffAbs->Divide(protonRawPrim,protonMC);
1075 protonEffAbs->SetName(Form("ProtonEffAbs_%i",iEta));
1076 protonEffAbs->Write();
1080 TH1D *protonEffMis = new TH1D(*protonMis);
1081 protonEffMis->Divide(protonMis,protonRaw);
1082 protonEffMis->SetName(Form("ProtonContam_%i",iEta));
1083 protonEffMis->Write();
1086 // PROCESS THE REAL DATA FROM HERE ON
1088 TH1D *histPtProtonEta = histPtEtaProton->ProjectionZ(Form("HistPtProtonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1089 histPtProtonEta->Write();
1091 // Apply corrections as you would like them to have...
1093 TH1D *protonFinal = new TH1D(*histPtProtonEta);
1094 protonFinal->SetName(Form("HistPtProtonEtaFinal_%i",iEta));
1095 protonFinal->Sumw2();
1096 protonFinal->Divide(protonEffTot); // enable or disable specific corrections here...
1097 protonFinal->Write();
1099 // if available save also the MC truth for consistency checks
1101 TH1D *histPtProtonEtaMcData = histPtMCProtonData->ProjectionZ(Form("HistPtProtonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1102 histPtProtonEtaMcData->Write();
1106 for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1107 TH1D *kaonRaw = histPtEtaKaonEff->ProjectionZ(Form("KaonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1109 TH1D *kaonMC = histPtMCKaonEff->ProjectionZ(Form("McKaon_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1111 TH1D *kaonRawPrim = histEffKaon->ProjectionZ(Form("KaonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1112 kaonRawPrim->Sumw2();
1113 TH1D *kaonSecReac = histEffKaon->ProjectionZ(Form("KaonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1114 kaonSecReac->Sumw2();
1115 TH1D *kaonSecWeak = histEffKaon->ProjectionZ(Form("KaonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1116 kaonSecWeak->Sumw2();
1117 TH1D *kaonMis = histEffKaon->ProjectionZ(Form("KaonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1119 kaonSecReac->Add(kaonSecWeak,-1);
1121 // 1. total (physical) efficiency
1123 TH1D *kaonEffTot = new TH1D(*kaonRaw);
1124 kaonEffTot->Divide(kaonRaw,kaonMC);
1125 kaonEffTot->SetName(Form("KaonEffTot_%i",iEta));
1126 kaonEffTot->Write();
1128 // 2. contributions from secondaries created by interactions
1130 TH1D *kaonEffSecReac = new TH1D(*kaonRaw);
1131 TH1D *dummy = new TH1D(*kaonRaw);
1132 dummy->Add(kaonSecReac,-1);
1133 kaonEffSecReac->Divide(kaonRaw,dummy);
1134 kaonEffSecReac->SetName(Form("KaonEffSecReac_%i",iEta));
1135 kaonEffSecReac->Write();
1138 // 3. contributions from secondaries from weak decays
1140 TH1D *kaonEffSecWeak = new TH1D(*kaonRaw);
1141 TH1D *dummy2 = new TH1D(*kaonRaw);
1142 dummy2->Add(kaonSecWeak,-1);
1143 kaonEffSecWeak->Divide(kaonRaw,dummy2);
1144 kaonEffSecWeak->SetName(Form("KaonEffSecWeak_%i",iEta));
1145 kaonEffSecWeak->Write();
1148 // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1150 TH1D *kaonEffAbs = new TH1D(*kaonRaw);
1151 kaonEffAbs->Divide(kaonRawPrim,kaonMC);
1152 kaonEffAbs->SetName(Form("KaonEffAbs_%i",iEta));
1153 kaonEffAbs->Write();
1157 TH1D *kaonEffMis = new TH1D(*kaonMis);
1158 kaonEffMis->Divide(kaonMis,kaonRaw);
1159 kaonEffMis->SetName(Form("KaonContam_%i",iEta));
1160 kaonEffMis->Write();
1163 // PROCESS THE REAL DATA FROM HERE ON
1165 TH1D *histPtKaonEta = histPtEtaKaon->ProjectionZ(Form("HistPtKaonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1166 histPtKaonEta->Write();
1168 // Apply corrections as you would like them to have...
1170 TH1D *kaonFinal = new TH1D(*histPtKaonEta);
1171 kaonFinal->SetName(Form("HistPtKaonEtaFinal_%i",iEta));
1173 kaonFinal->Divide(kaonEffTot); // enable or disable specific corrections here...
1176 // if available save also the MC truth for consistency checks
1178 TH1D *histPtKaonEtaMcData = histPtMCKaonData->ProjectionZ(Form("HistPtKaonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1179 histPtKaonEtaMcData->Write();
1184 for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1185 TH1D *pionRaw = histPtEtaPionEff->ProjectionZ(Form("PionRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1187 TH1D *pionMC = histPtMCPionEff->ProjectionZ(Form("McPion_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1189 TH1D *pionRawPrim = histEffPion->ProjectionZ(Form("PionRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1190 pionRawPrim->Sumw2();
1191 TH1D *pionSecReac = histEffPion->ProjectionZ(Form("PionSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1192 pionSecReac->Sumw2();
1193 TH1D *pionSecWeak = histEffPion->ProjectionZ(Form("PionSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1194 pionSecWeak->Sumw2();
1195 TH1D *pionMis = histEffPion->ProjectionZ(Form("PionMis_%i",iEta),3,3,iEta,iEta); // thin black line
1197 pionSecReac->Add(pionSecWeak,-1);
1199 // 1. total (physical) efficiency
1201 TH1D *pionEffTot = new TH1D(*pionRaw);
1202 pionEffTot->Divide(pionRaw,pionMC);
1203 pionEffTot->SetName(Form("PionEffTot_%i",iEta));
1204 pionEffTot->Write();
1206 // 2. contributions from secondaries created by interactions
1208 TH1D *pionEffSecReac = new TH1D(*pionRaw);
1209 TH1D *dummy = new TH1D(*pionRaw);
1210 dummy->Add(pionSecReac,-1);
1211 pionEffSecReac->Divide(pionRaw,dummy);
1212 pionEffSecReac->SetName(Form("PionEffSecReac_%i",iEta));
1213 pionEffSecReac->Write();
1216 // 3. contributions from secondaries from weak decays
1218 TH1D *pionEffSecWeak = new TH1D(*pionRaw);
1219 TH1D *dummy2 = new TH1D(*pionRaw);
1220 dummy2->Add(pionSecWeak,-1);
1221 pionEffSecWeak->Divide(pionRaw,dummy2);
1222 pionEffSecWeak->SetName(Form("PionEffSecWeak_%i",iEta));
1223 pionEffSecWeak->Write();
1226 // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1228 TH1D *pionEffAbs = new TH1D(*pionRaw);
1229 pionEffAbs->Divide(pionRawPrim,pionMC);
1230 pionEffAbs->SetName(Form("PionEffAbs_%i",iEta));
1231 pionEffAbs->Write();
1235 TH1D *pionEffMis = new TH1D(*pionMis);
1236 pionEffMis->Divide(pionMis,pionRaw);
1237 pionEffMis->SetName(Form("PionContam_%i",iEta));
1238 pionEffMis->Write();
1241 // PROCESS THE REAL DATA FROM HERE ON
1243 TH1D *histPtPionEta = histPtEtaPion->ProjectionZ(Form("HistPtPionEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1244 histPtPionEta->Write();
1246 // Apply corrections as you would like them to have...
1248 TH1D *pionFinal = new TH1D(*histPtPionEta);
1249 pionFinal->SetName(Form("HistPtPionEtaFinal_%i",iEta));
1251 pionFinal->Divide(pionEffTot); // enable or disable specific corrections here...
1254 // if available save also the MC truth for consistency checks
1256 TH1D *histPtPionEtaMcData = histPtMCPionData->ProjectionZ(Form("HistPtPionEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1257 histPtPionEtaMcData->Write();
1262 // Close the file after everthing is done
1264 spectraFile.Close();