Taking into account that only 1 or 2 values may be present for the
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / AliAnalysisTaskChargedHadronSpectra.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
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  **************************************************************************/
16
17 ///////////////////////////////////////////////////////////////////////////
18 //                                                                       //
19 //                                                                       //
20 // Analysis for identified particle spectra measured with TPC dE/dx.     //
21 //                                                                       //
22 //                                                                       //
23 ///////////////////////////////////////////////////////////////////////////
24
25 #include "Riostream.h"
26 #include "TChain.h"
27 #include "TTree.h"
28 #include "TH1F.h"
29 #include "TH2F.h"
30 #include "TH3F.h"
31 #include "TList.h"
32 #include "TMath.h"
33 #include "TCanvas.h"
34 #include "TObjArray.h"
35 #include "TF1.h"
36 #include "TFile.h"
37
38 #include "AliAnalysisTaskSE.h"
39 #include "AliAnalysisManager.h"
40
41 #include "AliHeader.h"
42 #include "AliGenPythiaEventHeader.h"
43 #include "AliGenCocktailEventHeader.h"
44
45 #include "AliPID.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliESDVertex.h"
48 #include "AliESDEvent.h"
49 #include "AliESDInputHandler.h"
50 #include "AliESDtrack.h"
51 #include "AliTPCpidESD.h"
52
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
55 #include "AliStack.h"
56
57 #include "AliLog.h"
58
59 #include "AliAnalysisTaskChargedHadronSpectra.h"
60
61
62 ClassImp(AliAnalysisTaskChargedHadronSpectra)
63
64 //________________________________________________________________________
65 AliAnalysisTaskChargedHadronSpectra::AliAnalysisTaskChargedHadronSpectra() 
66   : AliAnalysisTaskSE("TaskChargedHadron"), fESD(0), fListHist(0), fESDtrackCuts(0),fPidObject(0),
67     fMCtrue(0),
68     fAlephParameters(),
69     fHistPtMCKaon(0),
70     fHistPtMCProton(0),
71     fHistPtMCPion(0),
72     fHistPtMCElectron(0),
73     fHistPtMCMuon(0),
74     fHistPtEtaKaon(0),
75     fHistPtEtaKaonNoKink(0),
76     fHistPtEtaProton(0),
77     fHistPtEtaProtonDCA(0),
78     fHistPtEtaPion(0),
79     fHistPtEtaElectron(0),
80     fHistClassicKaon(0),
81     fHistClassicProton(0),
82     fHistClassicPion(0),
83     fHistClassicElectron(0),
84     fDeDx(0),
85     fHistTrackPerEvent(0),
86     fHistTrackPerEventMC(0),
87     fSecProtons(0),
88     fVertexZ(0),
89     fHistEtaNcls(0),
90     fHistEtaPhi(0),
91     fHistEffProton(0),
92     fHistEffProtonDCA(0),
93     fHistEffPion(0),
94     fHistEffKaon(0),
95     fHighPtElectrons(0),
96     fHighPtHadrons(0)
97 {
98   // default Constructor
99 }
100
101
102 //________________________________________________________________________
103 AliAnalysisTaskChargedHadronSpectra::AliAnalysisTaskChargedHadronSpectra(const char *name) 
104   : AliAnalysisTaskSE(name), fESD(0), fListHist(0), fESDtrackCuts(0),fPidObject(0),
105     fMCtrue(0),
106     fAlephParameters(),
107     fHistPtMCKaon(0),
108     fHistPtMCProton(0),
109     fHistPtMCPion(0),
110     fHistPtMCElectron(0),
111     fHistPtMCMuon(0),
112     fHistPtEtaKaon(0),
113     fHistPtEtaKaonNoKink(0),
114     fHistPtEtaProton(0),
115     fHistPtEtaProtonDCA(0),
116     fHistPtEtaPion(0),
117     fHistPtEtaElectron(0),
118     fHistClassicKaon(0),
119     fHistClassicProton(0),
120     fHistClassicPion(0),
121     fHistClassicElectron(0),
122     fDeDx(0),
123     fHistTrackPerEvent(0),
124     fHistTrackPerEventMC(0),
125     fSecProtons(0),
126     fVertexZ(0),
127     fHistEtaNcls(0),
128     fHistEtaPhi(0),
129     fHistEffProton(0),
130     fHistEffProtonDCA(0),
131     fHistEffPion(0),
132     fHistEffKaon(0),
133     fHighPtElectrons(0),
134     fHighPtHadrons(0)
135 {
136   //
137   // standard constructur which should be used - PID objects is initialized
138   //
139
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;
146
147   fPidObject = new AliTPCpidESD();
148   fPidObject->SetBetheBlochParameters(fAlephParameters[0]/50.,fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
149   
150   // Constructor
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
155   
156   // Output slot #0 writes into a TList container
157   DefineOutput(1, TList::Class());
158
159  
160
161 }
162
163
164 //________________________________________________________________________
165 void AliAnalysisTaskChargedHadronSpectra::UserCreateOutputObjects() 
166 {
167   // Create histograms
168   // Called once
169   fListHist = new TList();
170   //fListHist->SetOwner(); // Whoever knows how the data handling is ...?
171
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;
180
181   // sort pT-bins ..
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]];
187   
188
189   // MC histograms
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);
195   //
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);
201   //
202   fListHist->Add(fHistPtMCKaon);
203   fListHist->Add(fHistPtMCProton);
204   fListHist->Add(fHistPtMCPion);
205   fListHist->Add(fHistPtMCElectron);
206   fListHist->Add(fHistPtMCMuon);
207
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);
215   //
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);
222   //
223   fListHist->Add(fHistPtEtaKaon);
224   fListHist->Add(fHistPtEtaKaonNoKink);
225   fListHist->Add(fHistPtEtaProton);
226   fListHist->Add(fHistPtEtaProtonDCA);
227   fListHist->Add(fHistPtEtaPion);
228   fListHist->Add(fHistPtEtaElectron);
229
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);
235   //
236   fHistClassicKaon->GetXaxis()->Set(kPtBins, binsPt);
237   fHistClassicProton->GetXaxis()->Set(kPtBins, binsPt);
238   fHistClassicPion->GetXaxis()->Set(kPtBins, binsPt);
239   fHistClassicElectron->GetXaxis()->Set(kPtBins, binsPt);
240   //
241   fListHist->Add(fHistClassicKaon);
242   fListHist->Add(fHistClassicProton);
243   fListHist->Add(fHistClassicPion);
244   fListHist->Add(fHistClassicElectron);
245
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);
248   BinLogX(fDeDx);
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);
258   //
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);
263
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);
273   //
274   fHistEffProton->GetZaxis()->Set(kPtBins, binsPt);
275   fHistEffProtonDCA->GetZaxis()->Set(kPtBins, binsPt);
276   fHistEffPion->GetZaxis()->Set(kPtBins, binsPt);
277   fHistEffKaon->GetZaxis()->Set(kPtBins, binsPt);
278   //
279
280   
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);
286
287   //
288   // Create Objects which are needed for the rest of the analysis
289   //
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);
299   //
300   fListHist->Add(fESDtrackCuts);
301   //
302   // 
303  
304 }
305
306 //________________________________________________________________________
307 void AliAnalysisTaskChargedHadronSpectra::UserExec(Option_t *) 
308 {
309   // main event loop
310   
311   AliLog::SetGlobalLogLevel(AliLog::kError);
312   //
313   // Check Monte Carlo information and other access first:
314   //
315   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
316   if (!eventHandler) {
317     //Printf("ERROR: Could not retrieve MC event handler");
318     if (fMCtrue) return;
319   }
320   
321   AliMCEvent* mcEvent = eventHandler->MCEvent();
322   if (!mcEvent) {
323     //Printf("ERROR: Could not retrieve MC event");
324     if (fMCtrue) return;
325   }
326   
327   fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
328   if (!fESD) {
329     //Printf("ERROR: fESD not available");
330     return;
331   }
332   
333   if (!fESDtrackCuts) {
334     Printf("ERROR: fESDtrackCuts not available");
335     return;
336   }
337   
338   Int_t trackCounter = 0;
339
340  // monitor vertex position
341  const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
342  if (!vertex) {
343    fHistTrackPerEvent->Fill(trackCounter);
344    PostData(1, fListHist);
345    return;
346  } else {
347    if (vertex->GetZv() != 0) fVertexZ->Fill(vertex->GetZv());
348    if (TMath::Abs(vertex->GetZv()) > 10) {
349      fHistTrackPerEvent->Fill(trackCounter);
350      PostData(1, fListHist);
351      return;
352    }   
353  }
354  
355  
356  
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) {
362      trackCounter++;
363      //
364    }
365    delete track;
366  }
367  
368  
369  // MC loop
370  AliStack* stack = mcEvent->Stack();
371  
372  if (fMCtrue) {
373    //
374    //
375    //
376    AliHeader * header = mcEvent->Header();
377    Int_t processtype = GetPythiaEventProcessType(header);
378    // non diffractive
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);
384    //
385    Int_t mult = trackCounter;//stack->GetNprimary();
386    fHistTrackPerEventMC->Fill(0., trackCounter);
387    //
388    for(Int_t i = 0; i < stack->GetNtrack(); i++) {
389      TParticle * trackMC = stack->Particle(i);
390      Int_t pdg = trackMC->GetPdgCode();
391      //
392      Double_t xv = trackMC->Vx();
393      Double_t yv = trackMC->Vy();
394      Double_t zv = trackMC->Vz();
395      Double_t dxy = 0;
396      dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
397      Double_t dz = 0;
398      dz = TMath::Abs(zv); // so stupid to avoid warnings
399      //
400      // vertex cut - selection of primaries
401      //
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
405      //
406      if (!stack->IsPhysicalPrimary(i)) continue;
407      //
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
412      //
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
417      //
418      // special Kaon checks - those which decay before entering the TPC fiducial volume
419      //
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
425      }
426    }
427  }
428  
429  //end MC tracks loop 
430  
431  
432  // ESD loop
433  // definition of PID functions --> to be put to CreateOutputObjects()
434  
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);
440
441  foProton.SetParameters(fAlephParameters);
442  foPion.SetParameters(fAlephParameters);
443  foElec.SetParameters(fAlephParameters);
444  foKaon.SetParameters(fAlephParameters);
445  foMuon.SetParameters(fAlephParameters);  
446  
447  const Float_t k2sigmaCorr = 1/0.9545;
448
449  Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
450  
451  for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
452    //
453    // 1. select TPC tracks only and propagate them to the primary vertex determined with the TPC standalone
454    //
455    AliESDtrack *track = AliESDtrackCuts::GetTPCOnlyTrack(fESD,i);
456    if (!track) continue;
457    AliESDtrack *trackDummy = fESD->GetTrack(i);
458    if (!trackDummy->GetInnerParam()) { // DEBUG DEBUG
459      delete track;
460      continue;
461    }
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());
466    //
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;
471    //
472    // 2.a) apply some standard track cuts according to general recommendations
473    //
474  
475    if (!fESDtrackCuts->AcceptTrack(track)) {
476      //
477      if (fMCtrue) {
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
481      }
482      //
483      delete track;
484      continue;
485    }
486    
487    //
488    // 3. make the PID
489    //
490    Double_t sign = track->GetSign();
491    
492    Double_t tpcSignal = track->GetTPCsignal();
493    //
494    fDeDx->Fill(ptot, tpcSignal);
495    //
496    //
497    //
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));
502    //
503    // fill histograms for dNdpT
504    //
505    if (eta < 0.15) {
506      fHighPtHadrons->Fill(pT);
507      Double_t delta = (tpcSignal - foElec.Eval(ptot))/foElec.Eval(ptot);
508      if (delta > 0) fHighPtElectrons->Fill(pT);
509    }
510    //
511    /* 2sigma PID with 2sigma eff correction! */
512    // PION
513    if (TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kPion)) < 2) {
514      fHistPtEtaPion->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
515      //
516      if (trackCounter < 300 && fMCtrue) {
517        TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
518        Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
519        if (pdg == 211 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffPion->Fill(0,eta,sign*pT,k2sigmaCorr);
520        if (pdg == 211 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
521          fHistEffPion->Fill(1,eta,sign*pT,k2sigmaCorr);
522          TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
523          Int_t code = TMath::Abs(trackMother->GetPdgCode());
524          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);
525        }
526        if (TMath::Abs(trackMC->GetPdgCode()) != 211) fHistEffPion->Fill(2,eta,sign*pT,k2sigmaCorr);
527        if (TMath::Abs(trackMC->GetPdgCode()) == 13)  fHistEffPion->Fill(4,eta,sign*pT,k2sigmaCorr);
528        }
529    }
530    // KAON
531    if (TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kKaon)) < 2) {
532      fHistPtEtaKaon->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
533      //
534      if (trackCounter < 300 && fMCtrue) {
535        TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
536        Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
537        if (pdg == 321 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) fHistEffKaon->Fill(0,eta,sign*pT,k2sigmaCorr);
538        if (pdg == 321 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
539          fHistEffKaon->Fill(1,eta,sign*pT,k2sigmaCorr);
540          TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
541          Int_t code = TMath::Abs(trackMother->GetPdgCode());
542          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);
543        }
544          if (TMath::Abs(trackMC->GetPdgCode()) != 321) fHistEffKaon->Fill(2,eta,sign*pT,k2sigmaCorr);
545      }
546    }
547    // KAON NO KINK
548    if (TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kKaon)) < 2 && track->GetKinkIndex(0)==0) fHistPtEtaKaonNoKink->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
549    // PROTON
550    if (TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kProton)) < 2) {
551      fHistPtEtaProton->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
552      //
553      track->GetImpactParameters(dca,cov);
554      Float_t primVtxDCA = TMath::Sqrt(dca[0]*dca[0] + dca[1]*dca[1]);
555      fHistPtEtaProtonDCA->Fill(primVtxDCA, eta, sign*pT, k2sigmaCorr);
556      //
557      if (trackCounter < 300 && fMCtrue) {
558        TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
559        Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
560        if (pdg == 2212 && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
561          fHistEffProton->Fill(0.,eta,sign*pT,k2sigmaCorr);
562          if (eta < 0.15) fHistEffProtonDCA->Fill(0.,primVtxDCA,sign*pT,k2sigmaCorr);
563        }
564        if (pdg == 2212 && !stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) {
565          fHistEffProton->Fill(1,eta,sign*pT,k2sigmaCorr);
566          if (eta < 0.15) fHistEffProtonDCA->Fill(1,primVtxDCA,sign*pT,k2sigmaCorr);
567          TParticle *trackMother = stack->Particle(trackMC->GetFirstMother());
568          Int_t code = TMath::Abs(trackMother->GetPdgCode());
569          if (code==3122||code==3222||code==3212||code==3112||code==3322||code==3312||code==3332||code==130||code==310) {
570            fHistEffProton->Fill(3,eta,sign*pT,k2sigmaCorr);
571            if (eta < 0.15) fHistEffProtonDCA->Fill(3,primVtxDCA,sign*pT,k2sigmaCorr);
572          }
573          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());
574        }
575        if (TMath::Abs(trackMC->GetPdgCode()) != 2212) fHistEffProton->Fill(2,eta,sign*pT,k2sigmaCorr);
576      }
577    }
578    // ELECTRON
579    if (TMath::Abs(fPidObject->GetNumberOfSigmas(track,AliPID::kElectron))) fHistPtEtaElectron->Fill(trackCounter, eta, sign*pT, k2sigmaCorr);
580    
581    delete track;
582    
583  } // end of track loop
584  
585  fHistTrackPerEvent->Fill(trackCounter);
586  
587  // Post output data
588  
589  
590  PostData(1, fListHist);
591
592 }      
593
594
595 //________________________________________________________________________
596 void AliAnalysisTaskChargedHadronSpectra::Terminate(Option_t *) 
597 {
598   // Draw result to the screen
599   // Called once at the end of the query
600   cout << "*** TERMINATE ***" << endl;
601
602 }
603
604
605
606
607 //________________________________________________________________________
608 void AliAnalysisTaskChargedHadronSpectra::BinLogX(const TH1 *h) {
609
610   // Method for the correct logarithmic binning of histograms
611
612   TAxis *axis = h->GetXaxis();
613   int bins = axis->GetNbins();
614
615   Double_t from = axis->GetXmin();
616   Double_t to = axis->GetXmax();
617   Double_t *newBins = new Double_t[bins + 1];
618    
619   newBins[0] = from;
620   Double_t factor = pow(to/from, 1./bins);
621   
622   for (int i = 1; i <= bins; i++) {
623    newBins[i] = factor * newBins[i-1];
624   }
625   axis->Set(bins, newBins);
626   delete newBins;
627   
628 }
629
630
631 //________________________________________________________________________
632 Int_t AliAnalysisTaskChargedHadronSpectra::GetPythiaEventProcessType(const AliHeader* aHeader, const Bool_t adebug) const {
633   //
634   // get the process type of the event.
635   //
636
637   // can only read pythia headers, either directly or from cocktalil header
638   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
639
640   if (!pythiaGenHeader) {
641
642     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
643     if (!genCocktailHeader) {
644       //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
645       return -1;
646     }
647
648     TList* headerList = genCocktailHeader->GetHeaders();
649     if (!headerList) {
650       return -1;
651     }
652
653     for (Int_t i=0; i<headerList->GetEntries(); i++) {
654       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
655       if (pythiaGenHeader)
656         break;
657     }
658
659     if (!pythiaGenHeader) {
660       //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Could not find Pythia header. \n");
661       return -1;
662     }
663   }
664
665   if (adebug) {
666     //printf("AliAnalysisTaskChargedHadronSpectra::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
667   }
668
669   return pythiaGenHeader->ProcessType();
670
671
672
673 //________________________________________________________________________
674 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicProton(const TH3F * input, Int_t EtaBin, const Double_t * AlephParams) {
675   //
676   // The multiple Gauss fitting is happening here...
677   // input histogram is of the form (Pt, eta, difference in dEdx)
678   //
679
680   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20); 
681   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
682   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
683   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
684   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
685   //
686   foProton->SetParameters(AlephParams);  
687   foPion->SetParameters(AlephParams);
688   foElec->SetParameters(AlephParams);  
689   foKaon->SetParameters(AlephParams);  
690   foMuon->SetParameters(AlephParams);
691
692   
693   const Double_t kSigma = 0.06; // expected resolution (roughly)
694
695   const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
696   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
697   const Int_t kEtaBins = input->GetYaxis()->GetNbins();
698
699   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
700   // sort pT-bins ..
701   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};
702   Int_t indexes[kPtBins+1];
703   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
704   Double_t binsPt[kPtBins+1];
705   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
706   histPt->GetXaxis()->Set(kPtBins, binsPt);
707   //
708
709   TCanvas * canvMany = new TCanvas("canvManyProton","canvManyProton");
710   canvMany->Print("canvManyProton.ps[");
711
712   for(Int_t x=1; x < kPtBins+1; x++) {
713     Double_t pT = input->GetXaxis()->GetBinCenter(x);
714     for(Int_t y=1; y < kEtaBins+1; y++) {
715       Double_t eta = input->GetYaxis()->GetBinCenter(y);
716       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
717       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
718       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
719       histDeDx->SetTitle(Form("pt_%f",pT));
720       Float_t maxYield = histDeDx->GetEntries()/(TMath::Sqrt(2*TMath::Pi())*0.04/histDeDx->GetBinWidth(1));
721       //
722       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-1/*-7*kSigma*/,1/*7*kSigma*/);
723       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};
724       gausFit->SetParameters(parameters);
725       gausFit->SetParLimits(0,0.,maxYield);
726       gausFit->SetParLimits(1,-0.05,0.05);
727       gausFit->SetParLimits(2,0.04,0.08);
728       //
729       Double_t pionPosition = (foPion->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
730       gausFit->SetParLimits(0,0.,maxYield);
731       gausFit->SetParLimits(4, 0.95*pionPosition, 1.05*pionPosition);
732       gausFit->SetParLimits(5,0.8*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foPion->Eval(ptot)/foProton->Eval(ptot)));
733       //
734       Double_t elecPosition = (foElec->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
735       gausFit->SetParLimits(0,0.,maxYield);
736       gausFit->SetParLimits(7, 0.95*elecPosition, 1.05*elecPosition);
737       gausFit->SetParLimits(8,0.8*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foElec->Eval(ptot)/foProton->Eval(ptot)));
738       //
739       Double_t kaonPosition = (foKaon->Eval(ptot)-foProton->Eval(ptot))/foProton->Eval(ptot);
740       gausFit->SetParLimits(0,0.,maxYield);
741       gausFit->SetParLimits(10, 0.95*kaonPosition, 1.05*kaonPosition);
742       gausFit->SetParLimits(11,0.8*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)),1.2*kSigma*(foKaon->Eval(ptot)/foProton->Eval(ptot)));
743       histDeDx->Fit("gausFit","QNR");
744       gausFit->GetParameters(parameters);
745       // visualization      
746       if (y == EtaBin) {
747         canvMany->cd(x);
748         gPad->SetLogy();  
749         histDeDx->SetMarkerStyle(21);
750         histDeDx->SetMarkerSize(0.7);
751         histDeDx->Draw("E");
752         gausFit->SetLineWidth(2);
753         gausFit->Draw("same");
754         //
755         TF1 gausFit0("gausFit0", "gaus(0)",-1,1);
756         TF1 gausFit1("gausFit1", "gaus(0)",-1,1);
757         TF1 gausFit2("gausFit2", "gaus(0)",-1,1);
758         TF1 gausFit3("gausFit3", "gaus(0)",-1,1);
759         gausFit0.SetParameters(parameters[0],parameters[1],parameters[2]);
760         gausFit1.SetParameters(parameters[3],parameters[4],parameters[5]);
761         gausFit2.SetParameters(parameters[6],parameters[7],parameters[8]);
762         gausFit3.SetParameters(parameters[9],parameters[10],parameters[11]);
763         //
764         gausFit0.SetLineColor(4);
765         gausFit1.SetLineColor(2);
766         gausFit2.SetLineColor(6);
767         gausFit3.SetLineColor(8);
768         //
769         gausFit0.SetLineWidth(1);
770         gausFit1.SetLineWidth(1);
771         gausFit2.SetLineWidth(1);
772         gausFit3.SetLineWidth(1);
773         //
774         gausFit0.Draw("same");
775         gausFit1.Draw("same"); 
776         gausFit2.Draw("same"); 
777         gausFit3.Draw("same");
778         canvMany->Print("canvManyProton.ps");
779       }
780       
781       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2)/histDeDx->GetBinWidth(1); // area of the gaus fit
782       delete gausFit;
783       //
784       // stupid solution --> remove as soon as possible
785       /*yield = 0;
786       TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
787       gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
788       for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
789         Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
790         yield += gausYield->Eval(delta);
791       }*/
792       //
793       if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
794       //delete gausYield;
795     }
796   }
797   canvMany->Print("canvManyProton.ps]");
798
799   TCanvas * canvPt = new TCanvas();
800   canvPt->cd();
801   histPt->Draw();
802   
803   return histPt;
804 }
805
806
807
808 //________________________________________________________________________
809 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicPion(const TH3F * input, Int_t EtaBin,const  Double_t * AlephParams) {
810   //
811   // The multiple Gauss fitting is happening here...
812   //
813   
814   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20); 
815   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
816   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
817   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
818   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
819   //
820   foProton->SetParameters(AlephParams);  
821   foPion->SetParameters(AlephParams);
822   foElec->SetParameters(AlephParams);  
823   foKaon->SetParameters(AlephParams);  
824   foMuon->SetParameters(AlephParams);
825
826   // input histogram is of the form (Pt, eta, difference in dEdx)
827
828   const Double_t kSigma = 0.06; // expected resolution (roughly)
829
830   const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
831   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
832   const Int_t kEtaBins = input->GetYaxis()->GetNbins();
833  
834   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
835   // sort pT-bins ..
836   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};
837   Int_t indexes[kPtBins+1];
838   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
839   Double_t binsPt[kPtBins+1];
840   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
841   histPt->GetXaxis()->Set(kPtBins, binsPt);
842   //
843
844   TCanvas * canvMany = new TCanvas();
845   canvMany->Divide(8,5);
846
847   for(Int_t x=1; x < kPtBins+1; x++) {
848     Double_t pT = input->GetXaxis()->GetBinCenter(x);
849     for(Int_t y=1; y < kEtaBins+1; y++) {
850       Double_t eta = input->GetYaxis()->GetBinCenter(y);
851       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
852       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
853       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
854       //
855       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-5*kSigma,5*kSigma);
856       //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};
857       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};
858       gausFit->SetParameters(parameters);
859       gausFit->SetParLimits(1,-0.05,0.05);//gausFit->FixParameter(1,0); DEBUG DEBUG DEBUG
860       gausFit->SetParLimits(2,0.04,0.08);//gausFit->FixParameter(2,kSigma); DEBUG DEBUG DEBUG
861       //
862       gausFit->FixParameter(4,(foProton->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
863       gausFit->FixParameter(5,kSigma);
864       gausFit->FixParameter(7,(foElec->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
865       gausFit->FixParameter(8,kSigma);
866       gausFit->FixParameter(10,(foKaon->Eval(ptot)-foPion->Eval(ptot))/foPion->Eval(ptot));
867       gausFit->FixParameter(11,kSigma);
868       histDeDx->Fit("gausFit","Q0");
869       if (y == EtaBin && pT < 0) {
870         canvMany->cd(x);        
871         gPad->SetLogy();
872         histDeDx->Draw();
873         gausFit->Draw("same");
874       }
875       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2); // area of the gaus fit
876       // stupid solution --> remove as soon as possible
877       yield = 0;
878       TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
879       gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
880       for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
881         Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
882         yield += gausYield->Eval(delta);
883       }
884       //
885       if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
886       //delete gausFit;
887       delete gausYield;      
888     }
889   }
890   
891   TCanvas * canvPt = new TCanvas();
892   canvPt->cd();
893   histPt->Draw();
894
895   return histPt;
896   
897 }
898
899
900 //________________________________________________________________________
901 TH1D * AliAnalysisTaskChargedHadronSpectra::AnalyseClassicKaon(const TH3F * input, Int_t EtaBin,const  Double_t * AlephParams) {
902   //
903   // The multiple Gauss fitting is happening here...
904   //
905
906   TF1 *foProton = new TF1("foProton", "AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20); 
907   TF1 *foPion = new TF1("foPion", "AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20);
908   TF1 *foElec = new TF1("foElec", "AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20);
909   TF1 *foKaon = new TF1("foKaon", "AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20);
910   TF1 *foMuon = new TF1("foMuon", "AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20);
911   //
912   foProton->SetParameters(AlephParams);  
913   foPion->SetParameters(AlephParams);
914   foElec->SetParameters(AlephParams);  
915   foKaon->SetParameters(AlephParams);  
916   foMuon->SetParameters(AlephParams);
917
918   const Double_t kSigma = 0.06; // expected resolution (roughly)
919
920   const Int_t kPtBins = 2*56/*input->GetXaxis()->GetNbins()*/;
921   const Double_t kPtMax  = input->GetXaxis()->GetXmax();
922   const Int_t kEtaBins = input->GetYaxis()->GetNbins();
923
924   TH1D * histPt = new TH1D("histPt", "Pt; pT (Gev); dN",kPtBins,-kPtMax,kPtMax);
925   // sort pT-bins ..
926   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};
927   Int_t indexes[kPtBins+1];
928   TMath::Sort(kPtBins+1,binsPtDummy,indexes,kFALSE);
929   Double_t binsPt[kPtBins+1];
930   for(Int_t i=0; i<kPtBins+1; i++) binsPt[i] = binsPtDummy[indexes[i]];
931   histPt->GetXaxis()->Set(kPtBins, binsPt);
932   //
933
934   TCanvas * canvMany = new TCanvas();
935   canvMany->Divide(8,5);
936
937   for(Int_t x=1; x < kPtBins+1; x++) {
938     Double_t pT = input->GetXaxis()->GetBinCenter(x);
939     for(Int_t y=1; y < kEtaBins+1; y++) {
940       Double_t eta = input->GetYaxis()->GetBinCenter(y);
941       Double_t theta = 2*TMath::ATan(TMath::Exp(-eta));
942       Double_t ptot = TMath::Abs(pT/TMath::Sin(theta));
943       TH1D *histDeDx = input->ProjectionZ(Form("dedx_%i_%i",x,y),x,x,y,y);
944       //
945       TF1 * gausFit = new TF1("gausFit", "gaus(0) + gaus(3) + gaus(6) + gaus(9)",-5*kSigma,5*kSigma);
946       //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};
947       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};
948       gausFit->SetParameters(parameters);
949       gausFit->SetParLimits(1,-0.05,0.05);//gausFit->FixParameter(1,0); DEBUG DEBUG DEBUG
950       gausFit->SetParLimits(2,0.04,0.08);//gausFit->FixParameter(2,kSigma); DEBUG DEBUG DEBUG
951       //
952       gausFit->FixParameter(4,(foProton->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot));
953       gausFit->FixParameter(5,kSigma);
954       gausFit->FixParameter(7,(foElec->Eval(ptot)-foKaon->Eval(ptot))/foKaon->Eval(ptot));
955       gausFit->FixParameter(8,kSigma);
956       gausFit->FixParameter(10,(foPion->Eval(ptot)- foKaon->Eval(ptot))/foKaon->Eval(ptot));
957       gausFit->FixParameter(11,kSigma);
958       histDeDx->Fit("gausFit","Q0");
959       if (y == EtaBin && pT < 0) {
960         canvMany->cd(x);
961         gPad->SetLogy();
962         histDeDx->Draw();
963         gausFit->Draw("same");
964       }
965       Double_t yield = gausFit->GetParameter(0)*TMath::Sqrt(2*TMath::Pi())*gausFit->GetParameter(2); // area of the gaus fit
966       // stupid solution --> remove as soon as possible
967       yield = 0;
968       TF1 * gausYield = new TF1("gausYield", "gaus(0)",-5*kSigma,5*kSigma);
969       gausYield->SetParameters(gausFit->GetParameter(0),gausFit->GetParameter(1),gausFit->GetParameter(2));
970       for(Int_t i=1; i < histDeDx->GetXaxis()->GetNbins(); i++) {
971         Double_t delta = histDeDx->GetXaxis()->GetBinCenter(i);
972         yield += gausYield->Eval(delta);
973       }
974       //
975       if (y == EtaBin && yield > 0) histPt->Fill(pT,yield);
976       //delete gausFit;
977       delete gausYield;
978     }
979   }
980
981   TCanvas * canvPt = new TCanvas();
982   canvPt->cd();
983   histPt->Draw();
984   
985   return histPt;
986
987 }
988
989
990 //________________________________________________________________________
991 void AliAnalysisTaskChargedHadronSpectra::Postprocess(const TList * ListOfHistogramsMC,const  TList * ListOfHistogramsData, const Char_t *filename) {
992
993   // Define ranges
994   const Int_t kMultLow  = 1;
995   const Int_t kMultHigh = 10;
996
997   const Int_t kEtaHigh  = 4;
998
999   // Extract histograms for the MC efficiency study
1000   TH3F* histPtMCKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCKaon");
1001   TH3F* histPtMCProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCProton");
1002   TH3F* histPtMCPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtMCPion");
1003   
1004   TH3F* histPtEtaKaonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaKaon");
1005   TH3F* histPtEtaProtonEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaProton");
1006   TH3F* histPtEtaPionEff = ( TH3F *)ListOfHistogramsMC->FindObject("HistPtEtaPion");
1007   
1008   TH3F* histEffProton = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffProton");
1009   TH3F* histEffKaon = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffKaon");
1010   TH3F* histEffPion = ( TH3F *)ListOfHistogramsMC->FindObject("HistEffPion");
1011
1012   // Extract data for the final spectra  
1013   TH3F* histPtEtaKaon = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaKaon");
1014   TH3F* histPtEtaProton = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaProton");
1015   TH3F* histPtEtaPion = ( TH3F *)ListOfHistogramsData->FindObject("HistPtEtaPion");
1016
1017   // "MC" of the real data for consistency checks
1018   TH3F* histPtMCKaonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCKaon");
1019   TH3F* histPtMCProtonData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCProton");
1020   TH3F* histPtMCPionData = ( TH3F *)ListOfHistogramsData->FindObject("HistPtMCPion");
1021
1022   TFile spectraFile(filename,"RECREATE");
1023
1024   // 1. Protons  
1025   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1026     TH1D *protonRaw     = histPtEtaProtonEff->ProjectionZ(Form("ProtonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1027     protonRaw->Sumw2();
1028     TH1D *protonMC      = histPtMCProtonEff->ProjectionZ(Form("McProton_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1029     protonMC->Sumw2();
1030     TH1D *protonRawPrim = histEffProton->ProjectionZ(Form("ProtonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1031     protonRawPrim->Sumw2();
1032     TH1D *protonSecReac = histEffProton->ProjectionZ(Form("ProtonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1033     protonSecReac->Sumw2();
1034     TH1D *protonSecWeak = histEffProton->ProjectionZ(Form("ProtonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1035     protonSecWeak->Sumw2();
1036     TH1D *protonMis     = histEffProton->ProjectionZ(Form("ProtonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1037     protonMis->Sumw2();
1038     protonSecReac->Add(protonSecWeak,-1);
1039     //
1040     // 1. total (physical) efficiency
1041     //
1042     TH1D *protonEffTot = new TH1D(*protonRaw);
1043     protonEffTot->Divide(protonRaw,protonMC);
1044     protonEffTot->SetName(Form("ProtonEffTot_%i",iEta));
1045     protonEffTot->Write();
1046     //
1047     // 2. contributions from secondaries created by interactions
1048     //
1049     TH1D *protonEffSecReac = new TH1D(*protonRaw);
1050     TH1D *dummy = new TH1D(*protonRaw);
1051     dummy->Add(protonSecReac,-1);
1052     protonEffSecReac->Divide(protonRaw,dummy);
1053     protonEffSecReac->SetName(Form("ProtonEffSecReac_%i",iEta));
1054     protonEffSecReac->Write();
1055     delete dummy;
1056     //
1057     // 3. contributions from secondaries from weak decays
1058     //
1059     TH1D *protonEffSecWeak = new TH1D(*protonRaw);
1060     TH1D *dummy2 = new TH1D(*protonRaw);
1061     dummy2->Add(protonSecWeak,-1);
1062     protonEffSecWeak->Divide(protonRaw,dummy2);
1063     protonEffSecWeak->SetName(Form("ProtonEffSecWeak_%i",iEta));
1064     protonEffSecWeak->Write();
1065     delete dummy2;
1066     //
1067     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1068     //
1069     TH1D *protonEffAbs = new TH1D(*protonRaw);
1070     protonEffAbs->Divide(protonRawPrim,protonMC);
1071     protonEffAbs->SetName(Form("ProtonEffAbs_%i",iEta));
1072     protonEffAbs->Write();
1073     //
1074     // 5. contamination
1075     //
1076     TH1D *protonEffMis = new TH1D(*protonMis);
1077     protonEffMis->Divide(protonMis,protonRaw);
1078     protonEffMis->SetName(Form("ProtonContam_%i",iEta));
1079     protonEffMis->Write();
1080
1081     //
1082     // PROCESS THE REAL DATA FROM HERE ON
1083     //
1084     TH1D *histPtProtonEta = histPtEtaProton->ProjectionZ(Form("HistPtProtonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1085     histPtProtonEta->Write();
1086     //
1087     // Apply corrections as you would like them to have...
1088     //
1089     TH1D *protonFinal = new TH1D(*histPtProtonEta);
1090     protonFinal->SetName(Form("HistPtProtonEtaFinal_%i",iEta));
1091     protonFinal->Sumw2();
1092     protonFinal->Divide(protonEffTot); // enable or disable specific corrections here...
1093     protonFinal->Write();
1094     //
1095     // if available save also the MC truth for consistency checks
1096     //
1097     TH1D *histPtProtonEtaMcData = histPtMCProtonData->ProjectionZ(Form("HistPtProtonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1098     histPtProtonEtaMcData->Write();
1099   }
1100   
1101   // 2. Kaons
1102   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1103     TH1D *kaonRaw     = histPtEtaKaonEff->ProjectionZ(Form("KaonRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1104     kaonRaw->Sumw2();
1105     TH1D *kaonMC      = histPtMCKaonEff->ProjectionZ(Form("McKaon_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1106     kaonMC->Sumw2();
1107     TH1D *kaonRawPrim = histEffKaon->ProjectionZ(Form("KaonRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1108     kaonRawPrim->Sumw2();
1109     TH1D *kaonSecReac = histEffKaon->ProjectionZ(Form("KaonSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1110     kaonSecReac->Sumw2();
1111     TH1D *kaonSecWeak = histEffKaon->ProjectionZ(Form("KaonSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1112     kaonSecWeak->Sumw2();
1113     TH1D *kaonMis     = histEffKaon->ProjectionZ(Form("KaonMis_%i",iEta),3,3,iEta,iEta); // thin black line
1114     kaonMis->Sumw2();
1115     kaonSecReac->Add(kaonSecWeak,-1);
1116     //
1117     // 1. total (physical) efficiency
1118     //
1119     TH1D *kaonEffTot = new TH1D(*kaonRaw);
1120     kaonEffTot->Divide(kaonRaw,kaonMC);
1121     kaonEffTot->SetName(Form("KaonEffTot_%i",iEta));
1122     kaonEffTot->Write();
1123     //
1124     // 2. contributions from secondaries created by interactions
1125     //
1126     TH1D *kaonEffSecReac = new TH1D(*kaonRaw);
1127     TH1D *dummy = new TH1D(*kaonRaw);
1128     dummy->Add(kaonSecReac,-1);
1129     kaonEffSecReac->Divide(kaonRaw,dummy);
1130     kaonEffSecReac->SetName(Form("KaonEffSecReac_%i",iEta));
1131     kaonEffSecReac->Write();
1132     delete dummy;
1133     //
1134     // 3. contributions from secondaries from weak decays
1135     //
1136     TH1D *kaonEffSecWeak = new TH1D(*kaonRaw);
1137     TH1D *dummy2 = new TH1D(*kaonRaw);
1138     dummy2->Add(kaonSecWeak,-1);
1139     kaonEffSecWeak->Divide(kaonRaw,dummy2);
1140     kaonEffSecWeak->SetName(Form("KaonEffSecWeak_%i",iEta));
1141     kaonEffSecWeak->Write();
1142     delete dummy2;
1143     //
1144     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1145     //
1146     TH1D *kaonEffAbs = new TH1D(*kaonRaw);
1147     kaonEffAbs->Divide(kaonRawPrim,kaonMC);
1148     kaonEffAbs->SetName(Form("KaonEffAbs_%i",iEta));
1149     kaonEffAbs->Write();
1150     //
1151     // 5. contamination
1152     //
1153     TH1D *kaonEffMis = new TH1D(*kaonMis);
1154     kaonEffMis->Divide(kaonMis,kaonRaw);
1155     kaonEffMis->SetName(Form("KaonContam_%i",iEta));
1156     kaonEffMis->Write();
1157     
1158     //
1159     // PROCESS THE REAL DATA FROM HERE ON
1160     //
1161     TH1D *histPtKaonEta = histPtEtaKaon->ProjectionZ(Form("HistPtKaonEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1162     histPtKaonEta->Write();
1163     //
1164     // Apply corrections as you would like them to have...
1165     //
1166     TH1D *kaonFinal = new TH1D(*histPtKaonEta);
1167     kaonFinal->SetName(Form("HistPtKaonEtaFinal_%i",iEta));
1168     kaonFinal->Sumw2();
1169     kaonFinal->Divide(kaonEffTot); // enable or disable specific corrections here...
1170     kaonFinal->Write();
1171     //
1172     // if available save also the MC truth for consistency checks
1173     //
1174     TH1D *histPtKaonEtaMcData = histPtMCKaonData->ProjectionZ(Form("HistPtKaonEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1175     histPtKaonEtaMcData->Write();
1176   }
1177
1178
1179   // 3. Pions
1180   for(Int_t iEta = 1; iEta < kEtaHigh+1; iEta++) {
1181     TH1D *pionRaw     = histPtEtaPionEff->ProjectionZ(Form("PionRawTot_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // red crosses
1182     pionRaw->Sumw2();
1183     TH1D *pionMC      = histPtMCPionEff->ProjectionZ(Form("McPion_%i",iEta),kMultLow,kMultHigh,iEta,iEta); // solid black line
1184     pionMC->Sumw2();
1185     TH1D *pionRawPrim = histEffPion->ProjectionZ(Form("PionRawPrim_%i",iEta),1,1,iEta,iEta); // thin red line
1186     pionRawPrim->Sumw2();
1187     TH1D *pionSecReac = histEffPion->ProjectionZ(Form("PionSecReac_%i",iEta),2,2,iEta,iEta); // blue line
1188     pionSecReac->Sumw2();
1189     TH1D *pionSecWeak = histEffPion->ProjectionZ(Form("PionSecWeak_%i",iEta),4,4,iEta,iEta); // green line
1190     pionSecWeak->Sumw2();
1191     TH1D *pionMis     = histEffPion->ProjectionZ(Form("PionMis_%i",iEta),3,3,iEta,iEta); // thin black line
1192     pionMis->Sumw2();
1193     pionSecReac->Add(pionSecWeak,-1);
1194     //
1195     // 1. total (physical) efficiency
1196     //
1197     TH1D *pionEffTot = new TH1D(*pionRaw);
1198     pionEffTot->Divide(pionRaw,pionMC);
1199     pionEffTot->SetName(Form("PionEffTot_%i",iEta));
1200     pionEffTot->Write();
1201     //
1202     // 2. contributions from secondaries created by interactions
1203     //
1204     TH1D *pionEffSecReac = new TH1D(*pionRaw);
1205     TH1D *dummy = new TH1D(*pionRaw);
1206     dummy->Add(pionSecReac,-1);
1207     pionEffSecReac->Divide(pionRaw,dummy);
1208     pionEffSecReac->SetName(Form("PionEffSecReac_%i",iEta));
1209     pionEffSecReac->Write();
1210     delete dummy;
1211     //
1212     // 3. contributions from secondaries from weak decays
1213     //
1214     TH1D *pionEffSecWeak = new TH1D(*pionRaw);
1215     TH1D *dummy2 = new TH1D(*pionRaw);
1216     dummy2->Add(pionSecWeak,-1);
1217     pionEffSecWeak->Divide(pionRaw,dummy2);
1218     pionEffSecWeak->SetName(Form("PionEffSecWeak_%i",iEta));
1219     pionEffSecWeak->Write();
1220     delete dummy2;
1221     //
1222     // 4. decay (Kaon,pion) and detector absorption (energy loss and multiple scattering)
1223     //
1224     TH1D *pionEffAbs = new TH1D(*pionRaw);
1225     pionEffAbs->Divide(pionRawPrim,pionMC);
1226     pionEffAbs->SetName(Form("PionEffAbs_%i",iEta));
1227     pionEffAbs->Write();
1228     //
1229     // 5. contamination
1230     //
1231     TH1D *pionEffMis = new TH1D(*pionMis);
1232     pionEffMis->Divide(pionMis,pionRaw);
1233     pionEffMis->SetName(Form("PionContam_%i",iEta));
1234     pionEffMis->Write();
1235     
1236     //
1237     // PROCESS THE REAL DATA FROM HERE ON
1238     //
1239     TH1D *histPtPionEta = histPtEtaPion->ProjectionZ(Form("HistPtPionEtaRaw_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1240     histPtPionEta->Write();
1241     //
1242     // Apply corrections as you would like them to have...
1243     //
1244     TH1D *pionFinal = new TH1D(*histPtPionEta);
1245     pionFinal->SetName(Form("HistPtPionEtaFinal_%i",iEta));
1246     pionFinal->Sumw2();
1247     pionFinal->Divide(pionEffTot); // enable or disable specific corrections here...
1248     pionFinal->Write();
1249     //
1250     // if available save also the MC truth for consistency checks
1251     //
1252     TH1D *histPtPionEtaMcData = histPtMCPionData->ProjectionZ(Form("HistPtPionEtaMcData_%i",iEta),kMultLow,kMultHigh,iEta,iEta);
1253     histPtPionEtaMcData->Write();
1254
1255   }
1256
1257   //
1258   // Close the file after everthing is done
1259   //
1260   spectraFile.Close();
1261
1262 }
1263
1264
1265
1266