Added protection for unidentified tracks (Nicolo' Jacazio)
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / ITSsa / AliAnalysisTaskSEITSsaSpectra.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 ///////////////////////////////////////////////////////////////////////////
17 // AliAnalysisTaskSE for the extraction of the various histograms to
18 // study the pt spectra of identified hadrons:
19 // - log(dEdx)-log(dEdxBB) distributions for pions, kaons and protons in pt bins
20 // - Pt distributions of pions, kaons and protons with nSigma PID
21 // Authors: 
22 // E. Biolcati, biolcati@to.infn.it
23 // L. Milano, milano@to.infn.it
24 // F. Prino, prino@to.infn.it
25 ///////////////////////////////////////////////////////////////////////////
26
27 #include <TH1F.h>
28 #include <TF1.h>
29 #include <TRandom3.h>
30 #include <TH2F.h>
31 #include <TChain.h>
32 #include <TNtuple.h>
33 #include <TParticle.h>
34 #include <Rtypes.h>
35 #include "AliAnalysisTaskSE.h"
36 #include "AliAnalysisManager.h"
37 #include "AliAnalysisDataContainer.h"
38 #include "AliESDEvent.h"
39 #include "AliESDInputHandler.h"
40 #include "AliESDtrack.h"
41 #include "AliStack.h"
42 #include "AliMCEventHandler.h"
43 #include "AliMCEvent.h"
44 #include "AliPhysicsSelection.h"
45 #include "AliAnalysisTaskSEITSsaSpectra.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliCentrality.h"
48 #include "AliMultiplicity.h"
49 #include "AliESDUtils.h"
50 #include "AliITSPIDResponse.h"
51
52 ClassImp(AliAnalysisTaskSEITSsaSpectra)
53 /* $Id$ */
54
55 //________________________________________________________________________
56 AliAnalysisTaskSEITSsaSpectra::AliAnalysisTaskSEITSsaSpectra():
57 AliAnalysisTaskSE("taskITSsaSpectra"),
58   fESD(0),
59   fOutput(0),
60   fListCuts(0),
61   fHistNEvents(0),
62   fHistMult(0),
63   fHistCen(0),
64   fHistNTracks(0),
65   fHistNTracksPos(0),
66   fHistNTracksNeg(0),
67   fHistDEDX(0),
68   fHistDEDXdouble(0),
69   fHistBeforeEvSel(0),
70   fHistAfterEvSel(0),
71   fDCAxyCutFunc(0x0),
72   fDCAzCutFunc(0x0),
73   fITSPIDResponse(0),
74   fMinSPDPts(1),
75   fMinNdEdxSamples(3),
76   fMindEdx(0.),
77   fMinNSigma(1.5),
78   fMaxY(0.5),
79   fMaxChi2Clu(2.5),
80   fNSigmaDCAxy(7.),
81   fNSigmaDCAz(7.),
82   fEtaRange(0.8),
83   fLowMult(-1),
84   fUpMult(-1),
85   fLowCentrality(-1.0),
86   fUpCentrality(-1.0),
87   fMultEstimator(0),
88   fHImode(0),
89   fYear(2010),
90   fMC(kFALSE), 
91   fSmearMC(kFALSE),
92   fSmearP(0.),
93   fSmeardEdx(0.),
94   fRandGener(0),
95   fFillNtuple(kFALSE),
96   fLowEnergypp(kFALSE),
97   fNtupleNSigma(0),
98   fNtupleMC(0)
99 {
100   // Constructor
101   Double_t xbins[kNbins+1]={0.08,0.10,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.65,0.70,0.75,0.80,0.85,0.90,0.95,1.0};
102   for(Int_t iBin=0; iBin<kNbins+1; iBin++) fPtBinLimits[iBin]=xbins[iBin];
103   fRandGener=new TRandom3(0);
104   DefineInput(0, TChain::Class());
105   DefineOutput(1, TList::Class());
106   DefineOutput(2,TList::Class());
107   Printf("end of AliAnalysisTaskSEITSsaSpectra");
108 }
109
110 //___________________________________________________________________________
111 AliAnalysisTaskSEITSsaSpectra::~AliAnalysisTaskSEITSsaSpectra(){
112   // Destructor
113   if (fOutput) {
114     delete fOutput;
115     fOutput = 0;
116   }
117   if (fListCuts) {
118     delete fListCuts;
119     fListCuts = 0;
120   }
121   if(fRandGener) delete fRandGener;
122   if(fITSPIDResponse) delete fITSPIDResponse;
123   delete fDCAxyCutFunc;
124   delete fDCAzCutFunc;
125 }
126
127 //________________________________________________________________________
128 Double_t AliAnalysisTaskSEITSsaSpectra::CookdEdx(Double_t *s) const {
129   // truncated mean for the dEdx
130   Int_t nc=0; 
131   Double_t dedx[4]={0.,0.,0.,0.};
132   for (Int_t il=0; il<4; il++) { // count good (>0) dE/dx values
133     if(s[il]>fMindEdx){
134       dedx[nc]= s[il];
135       nc++;
136     }
137   }
138   if(nc<fMinNdEdxSamples) return -1.;
139   
140   Double_t tmp;
141   Int_t swap; // sort in ascending order 
142   do {
143     swap=0;
144     for (Int_t i=0; i<nc-1; i++) {
145       if (dedx[i]<=dedx[i+1]) continue;
146       tmp=dedx[i];
147       dedx[i]=dedx[i+1];
148       dedx[i+1]=tmp;
149       swap++;
150     } 
151   } while (swap);
152   
153   Double_t sumamp=0,sumweight=0;
154   Double_t weight[4]={1.,1.,0.,0.};
155   if(nc==3) weight[1]=0.5;
156   else if(nc<3) weight[1]=0.;
157   for (Int_t i=0; i<nc; i++) {
158     sumamp+= dedx[i]*weight[i];
159     sumweight+=weight[i];
160   }
161   return sumamp/sumweight;
162 }
163
164
165 //________________________________________________________________________
166 Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt) const {
167   // cut on transverse impact parameter updaated on 20-5-2010
168   // from the study of L. Milano, F. Prino on the ITS standalone tracks
169   // using the common binning of the TPC tracks
170
171   Double_t xyMax = fDCAxyCutFunc->Eval(pt); //in micron
172   if((TMath::Abs(impactXY)*10000)>xyMax) return kFALSE;  
173   Double_t zMax = fDCAzCutFunc->Eval(pt); //in micron
174   if((TMath::Abs(impactZ)*10000)>zMax) return kFALSE;
175   
176   return kTRUE;
177 }
178
179 //________________________________________________________________________
180 Double_t AliAnalysisTaskSEITSsaSpectra::Eta2y(Double_t pt, Double_t m, Double_t eta) const {
181   // convert eta to y
182   Double_t mt = TMath::Sqrt(m*m + pt*pt);
183   return TMath::ASinH(pt/mt*TMath::SinH(eta));
184 }
185
186
187 //________________________________________________________________________
188 void AliAnalysisTaskSEITSsaSpectra::Init(){
189   //
190   // Initialization
191   //
192   fListCuts=new TList();
193   fListCuts->SetOwner();
194   Double_t xyP[3];
195   Double_t zP[3];
196   if(fMC){
197     if(fYear==2009){
198       xyP[0]=88.63; //MC LHC10a12
199       xyP[1]=19.57;
200       xyP[2]=1.65;
201       zP[0]=140.98;
202       zP[1]=62.33;
203       zP[2]=1.15;
204     }else{
205       xyP[0]=36.; //MC LHC10d1
206       xyP[1]=43.9;
207       xyP[2]=1.3;
208       zP[0]=111.9;
209       zP[1]=59.8;
210       zP[2]=1.2;
211     }
212   }else{
213     if(fYear==2009){
214       xyP[0]=85.28;//DATA 900 GeV pass6
215       xyP[1]=25.78;
216       xyP[2]=1.55;
217       zP[0]=146.80;
218       zP[1]=70.07;
219       zP[2]=1.11;
220     }else{
221       xyP[0]=32.7;//DATA 7 TeV pass2
222       xyP[1]=44.8;
223       xyP[2]=1.3;
224       zP[0]=117.3;
225       zP[1]=66.8;
226       zP[2]=1.2;
227     }
228   }
229   fDCAxyCutFunc = new TF1("fDCAxyCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
230   for(Int_t ipar=0; ipar<3; ipar++) fDCAxyCutFunc->SetParameter(ipar,xyP[ipar]);
231   fDCAxyCutFunc->SetParameter(3,fNSigmaDCAxy);
232
233
234                                                 
235   fDCAzCutFunc = new TF1("fDCAzCutFunc","[3]*([0]+[1]/TMath::Power(TMath::Abs(x),[2]))",0.05,10.);
236   for(Int_t ipar=0; ipar<3; ipar++) fDCAzCutFunc->SetParameter(ipar,zP[ipar]);
237   fDCAzCutFunc->SetParameter(3,fNSigmaDCAz);
238
239   fListCuts->Add(fDCAxyCutFunc);
240   fListCuts->Add(fDCAzCutFunc);
241
242   PostData(2,fListCuts);
243 }
244
245 //________________________________________________________________________
246 void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){
247   // Create a TList with histograms and a TNtuple
248   // Called once
249
250   fOutput = new TList();
251   fOutput->SetOwner();
252   fOutput->SetName("Spiderman");
253   
254   fHistNEvents = new TH1F("fHistNEvents", "Number of processed events",9,-1.5,7.5);
255   fHistNEvents->Sumw2();
256   fHistNEvents->SetMinimum(0);
257   fHistNEvents->GetXaxis()->SetBinLabel(1,"Read from ESD");
258   fHistNEvents->GetXaxis()->SetBinLabel(2,"Pass Phys. Sel. + Trig");
259   fHistNEvents->GetXaxis()->SetBinLabel(3,"SDD read out");
260   fHistNEvents->GetXaxis()->SetBinLabel(4,"In mult. range");
261   fHistNEvents->GetXaxis()->SetBinLabel(5,"With SPD vertex");
262   fHistNEvents->GetXaxis()->SetBinLabel(6,"Vertex contributors >0");
263   fHistNEvents->GetXaxis()->SetBinLabel(7,"|zVertex|<10");
264   fHistNEvents->GetXaxis()->SetBinLabel(8,"Error on zVertex<0.5");
265   fHistNEvents->GetXaxis()->SetBinLabel(9,"Good Z vertex");
266   fOutput->Add(fHistNEvents);
267   
268   fHistMult = new TH1F("fHistMult", "Event Multiplicity",3000,-0.5,2999.5);
269   fHistMult->Sumw2();
270   fHistMult->SetMinimum(0);
271   if(fMultEstimator==0) fHistMult->GetXaxis()->SetTitle("Multiplicity |#eta|<0.8");
272   else if(fMultEstimator==1) fHistMult->GetXaxis()->SetTitle("Tracklets |#eta|<0.8");
273   else if(fMultEstimator==2) fHistMult->GetXaxis()->SetTitle("Clusters on SPD1");
274   fOutput->Add(fHistMult);
275   
276   fHistCen = new TH1F("fHistCen", "Event Centrality",101,-0.5,100.5);
277   fHistCen->Sumw2();
278   fHistCen->SetMinimum(0);
279   fOutput->Add(fHistCen);
280   
281   fHistNTracks = new TH1F("fHistNTracks", "Number of ITSsa tracks",20,0.5,20.5);
282   fHistNTracks->Sumw2();
283   fHistNTracks->SetMinimum(0);
284   fOutput->Add(fHistNTracks);
285   
286   fHistNTracksPos = new TH1F("fHistNTracksPos", "Number of positive ITSsa tracks",20,0.5,20.5);
287   fHistNTracksPos->Sumw2();
288   fHistNTracksPos->SetMinimum(0);
289   fOutput->Add(fHistNTracksPos);
290   
291   fHistNTracksNeg = new TH1F("fHistNTracksNeg", "Number of negative ITSsa tracks",20,0.5,20.5);
292   fHistNTracksNeg->Sumw2();
293   fHistNTracksNeg->SetMinimum(0);
294   fOutput->Add(fHistNTracksNeg);
295   
296   //binning for the histogram
297   const Int_t hnbins=400;
298   Double_t hxmin = 0.01;
299   Double_t hxmax = 10;
300   Double_t hlogxmin = TMath::Log10(hxmin);
301   Double_t hlogxmax = TMath::Log10(hxmax);
302   Double_t hbinwidth = (hlogxmax-hlogxmin)/hnbins;
303   Double_t hxbins[hnbins+1];
304   hxbins[0] = 0.01; 
305   for (Int_t i=1;i<=hnbins;i++) {
306     hxbins[i] = hxmin + TMath::Power(10,hlogxmin+i*hbinwidth);
307   }
308   
309   fHistDEDX = new TH2F("fHistDEDX","",hnbins,hxbins,900,0,1000);
310   fOutput->Add(fHistDEDX);
311   
312   fHistDEDXdouble = new TH2F("fHistDEDXdouble","",500,-5,5,900,0,1000);
313   fOutput->Add(fHistDEDXdouble);
314   
315   fHistBeforeEvSel = new TH1F("fHistBeforeEvSel","fHistBeforeEvSel",kNbins,fPtBinLimits);
316   fHistAfterEvSel = new TH1F("fHistAfterEvSel","fHistAfterEvSel",kNbins,fPtBinLimits);
317   fOutput->Add(fHistBeforeEvSel);
318   fOutput->Add(fHistAfterEvSel);
319   
320   
321   
322   for(Int_t j=0;j<3;j++){
323     
324     fHistPosNSigmaSep[j] = new TH2F(Form("fHistPosNSigmaSep%d",j),"",hnbins,hxbins,1000,-10,10);
325     fOutput->Add(fHistPosNSigmaSep[j]);
326     fHistNegNSigmaSep[j] = new TH2F(Form("fHistNegNSigmaSep%d",j),"",hnbins,hxbins,1000,-10,10);
327     fOutput->Add(fHistNegNSigmaSep[j]);
328     
329     fHistPrimMCpos[j] = new TH1F(Form("fHistPrimMCpos%d",j),Form("fHistPrimMCpos%d",j),kNbins,fPtBinLimits);
330     fHistPrimMCneg[j] = new TH1F(Form("fHistPrimMCneg%d",j),Form("fHistPrimMCneg%d",j),kNbins,fPtBinLimits);
331     fOutput->Add(fHistPrimMCneg[j]);
332     fOutput->Add(fHistPrimMCpos[j]);
333     fHistSecStrMCpos[j] = new TH1F(Form("fHistSecStrMCpos%d",j),Form("fHistSecStrMCpos%d",j),kNbins,fPtBinLimits);
334     fHistSecStrMCneg[j] = new TH1F(Form("fHistSecStrMCneg%d",j),Form("fHistSecStrMCneg%d",j),kNbins,fPtBinLimits);
335     fOutput->Add(fHistSecStrMCneg[j]);
336     fOutput->Add(fHistSecStrMCpos[j]);
337     fHistSecMatMCpos[j] = new TH1F(Form("fHistSecMatMCpos%d",j),Form("fHistSecMatMCpos%d",j),kNbins,fPtBinLimits);
338     fHistSecMatMCneg[j] = new TH1F(Form("fHistSecMatMCneg%d",j),Form("fHistSecMatMCneg%d",j),kNbins,fPtBinLimits);
339     fOutput->Add(fHistSecMatMCneg[j]);
340     fOutput->Add(fHistSecMatMCpos[j]);
341     //
342     fHistPrimMCposBefEvSel[j] = new TH1F(Form("fHistPrimMCposBefEvSel%d",j),Form("fHistPrimMCposBefEvSel%d",j),kNbins,fPtBinLimits);
343     fHistPrimMCnegBefEvSel[j] = new TH1F(Form("fHistPrimMCnegBefEvSel%d",j),Form("fHistPrimMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
344     fOutput->Add(fHistPrimMCnegBefEvSel[j]);
345     fOutput->Add(fHistPrimMCposBefEvSel[j]);
346     fHistSecStrMCposBefEvSel[j] = new TH1F(Form("fHistSecStrMCposBefEvSel%d",j),Form("fHistSecStrMCposBefEvSel%d",j),kNbins,fPtBinLimits);
347     fHistSecStrMCnegBefEvSel[j] = new TH1F(Form("fHistSecStrMCnegBefEvSel%d",j),Form("fHistSecStrMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
348     fOutput->Add(fHistSecStrMCnegBefEvSel[j]);
349     fOutput->Add(fHistSecStrMCposBefEvSel[j]);
350     fHistSecMatMCposBefEvSel[j] = new TH1F(Form("fHistSecMatMCposBefEvSel%d",j),Form("fHistSecMatMCposBefEvSel%d",j),kNbins,fPtBinLimits);
351     fHistSecMatMCnegBefEvSel[j] = new TH1F(Form("fHistSecMatMCnegBefEvSel%d",j),Form("fHistSecMatMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
352     fOutput->Add(fHistSecMatMCnegBefEvSel[j]);
353     fOutput->Add(fHistSecMatMCposBefEvSel[j]);
354     //
355     fHistPrimMCposReco[j] = new TH1F(Form("fHistPrimMCposReco%d",j),Form("fHistPrimMCposReco%d",j),kNbins,fPtBinLimits);
356     fHistPrimMCnegReco[j] = new TH1F(Form("fHistPrimMCnegReco%d",j),Form("fHistPrimMCnegReco%d",j),kNbins,fPtBinLimits);
357     fOutput->Add(fHistPrimMCnegReco[j]);
358     fOutput->Add(fHistPrimMCposReco[j]);
359     fHistSecStrMCposReco[j] = new TH1F(Form("fHistSecStrMCposReco%d",j),Form("fHistSecStrMCposReco%d",j),kNbins,fPtBinLimits);
360     fHistSecStrMCnegReco[j] = new TH1F(Form("fHistSecStrMCnegReco%d",j),Form("fHistSecStrMCnegReco%d",j),kNbins,fPtBinLimits);
361     fOutput->Add(fHistSecStrMCnegReco[j]);
362     fOutput->Add(fHistSecStrMCposReco[j]);
363     fHistSecMatMCposReco[j] = new TH1F(Form("fHistSecMatMCposReco%d",j),Form("fHistSecMatMCposReco%d",j),kNbins,fPtBinLimits);
364     fHistSecMatMCnegReco[j] = new TH1F(Form("fHistSecMatMCnegReco%d",j),Form("fHistSecMatMCnegReco%d",j),kNbins,fPtBinLimits);
365     fOutput->Add(fHistSecMatMCnegReco[j]);
366     fOutput->Add(fHistSecMatMCposReco[j]);
367
368   }
369   
370   for(Int_t i=0; i<4; i++){
371     fHistCharge[i] = new TH1F(Form("fHistChargeLay%d",i),Form("fHistChargeLay%d",i),100,0,300);
372     fOutput->Add(fHistCharge[i]);
373   }
374   
375   for(Int_t i=0; i<kNbins; i++){
376     fHistPosPi[i] = new TH1F(Form("fHistPosPi%d",i),Form("fHistPosPi%d",i),175,-3.5,3.5);       
377     fHistPosK[i]  = new TH1F(Form("fHistPosK%d",i),Form("fHistPosK%d",i),175,-3.5,3.5); 
378     fHistPosP[i]  = new TH1F(Form("fHistPosP%d",i),Form("fHistPosP%d",i),175,-3.5,3.5); 
379     fHistNegPi[i] = new TH1F(Form("fHistNegPi%d",i),Form("fHistNegPi%d",i),175,-3.5,3.5);       
380     fHistNegK[i]  = new TH1F(Form("fHistNegK%d",i),Form("fHistNegK%d",i),175,-3.5,3.5); 
381     fHistNegP[i]  = new TH1F(Form("fHistNegP%d",i),Form("fHistNegP%d",i),175,-3.5,3.5); 
382     
383     fHistDCAPosPi[i] = new TH1F(Form("fHistDCAPosPi%d",i),Form("fHistDCAPosPi%d",i),2000,-1,1);  //DCA distr. with NSigma PID
384     fHistDCAPosK[i]  = new TH1F(Form("fHistDCAPosK%d",i),Form("fHistDCAPosK%d",i),2000,-1,1);   
385     fHistDCAPosP[i]  = new TH1F(Form("fHistDCAPosP%d",i),Form("fHistDCAPosP%d",i),2000,-1,1);   
386     fHistDCANegPi[i] = new TH1F(Form("fHistDCANegPi%d",i),Form("fHistDCANegPi%d",i),2000,-1,1); 
387     fHistDCANegK[i]  = new TH1F(Form("fHistDCANegK%d",i),Form("fHistDCANegK%d",i),2000,-1,1);   
388     fHistDCANegP[i]  = new TH1F(Form("fHistDCANegP%d",i),Form("fHistDCANegP%d",i),2000,-1,1);   
389     
390     fHistMCPrimDCAPosPi[i] = new TH1F(Form("fHistMCPrimDCAPosPi%d",i),Form("fHistMCPrimDCAPosPi%d",i),2000,-1,1);  //DCA distr. with MC truth
391     fHistMCPrimDCAPosK[i]  = new TH1F(Form("fHistMCPrimDCAPosK%d",i),Form("fHistMCPrimDCAPosK%d",i),2000,-1,1); 
392     fHistMCPrimDCAPosP[i]  = new TH1F(Form("fHistMCPrimDCAPosP%d",i),Form("fHistMCPrimDCAPosP%d",i),2000,-1,1); 
393     fHistMCPrimDCANegPi[i] = new TH1F(Form("fHistMCPrimDCANegPi%d",i),Form("fHistMCPrimDCANegPi%d",i),2000,-1,1);       
394     fHistMCPrimDCANegK[i]  = new TH1F(Form("fHistMCPrimDCANegK%d",i),Form("fHistMCPrimDCANegK%d",i),2000,-1,1); 
395     fHistMCPrimDCANegP[i]  = new TH1F(Form("fHistMCPrimDCANegP%d",i),Form("fHistMCPrimDCANegP%d",i),2000,-1,1); 
396     
397     fHistMCSecStDCAPosPi[i] = new TH1F(Form("fHistMCSecStDCAPosPi%d",i),Form("fHistMCSecStDCAPosPi%d",i),2000,-1,1);  //DCA distr. with MC truth
398     fHistMCSecStDCAPosK[i]  = new TH1F(Form("fHistMCSecStDCAPosK%d",i),Form("fHistMCSecStDCAPosK%d",i),2000,-1,1);      
399     fHistMCSecStDCAPosP[i]  = new TH1F(Form("fHistMCSecStDCAPosP%d",i),Form("fHistMCSecStDCAPosP%d",i),2000,-1,1);      
400     fHistMCSecStDCANegPi[i] = new TH1F(Form("fHistMCSecStDCANegPi%d",i),Form("fHistMCSecStDCANegPi%d",i),2000,-1,1);    
401     fHistMCSecStDCANegK[i]  = new TH1F(Form("fHistMCSecStDCANegK%d",i),Form("fHistMCSecStDCANegK%d",i),2000,-1,1);      
402     fHistMCSecStDCANegP[i]  = new TH1F(Form("fHistMCSecStDCANegP%d",i),Form("fHistMCSecStDCANegP%d",i),2000,-1,1);      
403     
404     fHistMCSecMatDCAPosPi[i] = new TH1F(Form("fHistMCSecMatDCAPosPi%d",i),Form("fHistMCSecMatDCAPosPi%d",i),2000,-1,1);  //DCA distr. with MC truth
405     fHistMCSecMatDCAPosK[i]  = new TH1F(Form("fHistMCSecMatDCAPosK%d",i),Form("fHistMCSecMatDCAPosK%d",i),2000,-1,1);   
406     fHistMCSecMatDCAPosP[i]  = new TH1F(Form("fHistMCSecMatDCAPosP%d",i),Form("fHistMCSecMatDCAPosP%d",i),2000,-1,1);   
407     fHistMCSecMatDCANegPi[i] = new TH1F(Form("fHistMCSecMatDCANegPi%d",i),Form("fHistMCSecMatDCANegPi%d",i),2000,-1,1); 
408     fHistMCSecMatDCANegK[i]  = new TH1F(Form("fHistMCSecMatDCANegK%d",i),Form("fHistMCSecMatDCANegK%d",i),2000,-1,1);   
409     fHistMCSecMatDCANegP[i]  = new TH1F(Form("fHistMCSecMatDCANegP%d",i),Form("fHistMCSecMatDCANegP%d",i),2000,-1,1);   
410     
411     fHistMCPosOtherHypPion[i] = new TH1F(Form("fHistMCPosOtherHypPion%d",i),Form("fHistMCPosOtherHypPion%d",i),175,-3.5,3.5);   //MC truth
412     fHistMCPosOtherHypKaon[i] = new TH1F(Form("fHistMCPosOtherHypKaon%d",i),Form("fHistMCPosOtherHypKaon%d",i),175,-3.5,3.5);
413     fHistMCPosOtherHypProton[i] = new TH1F(Form("fHistMCPosOtherHypProton%d",i),Form("fHistMCPosOtherHypProton%d",i),175,-3.5,3.5);
414     fHistMCPosElHypPion[i] = new TH1F(Form("fHistMCPosElHypPion%d",i),Form("fHistMCPosElHypPion%d",i),175,-3.5,3.5);    
415     fHistMCPosElHypKaon[i] = new TH1F(Form("fHistMCPosElHypKaon%d",i),Form("fHistMCPosElHypKaon%d",i),175,-3.5,3.5);
416     fHistMCPosElHypProton[i] = new TH1F(Form("fHistMCPosElHypProton%d",i),Form("fHistMCPosElHypProton%d",i),175,-3.5,3.5);
417     fHistMCPosPiHypPion[i] = new TH1F(Form("fHistMCPosPiHypPion%d",i),Form("fHistMCPosPiHypPion%d",i),175,-3.5,3.5);    
418     fHistMCPosPiHypKaon[i] = new TH1F(Form("fHistMCPosPiHypKaon%d",i),Form("fHistMCPosPiHypKaon%d",i),175,-3.5,3.5);
419     fHistMCPosPiHypProton[i] = new TH1F(Form("fHistMCPosPiHypProton%d",i),Form("fHistMCPosPiHypProton%d",i),175,-3.5,3.5);
420     fHistMCPosKHypPion[i] = new TH1F(Form("fHistMCPosKHypPion%d",i),Form("fHistMCPosKHypPion%d",i),175,-3.5,3.5);       
421     fHistMCPosKHypKaon[i]  = new TH1F(Form("fHistMCPosKHypKaon%d",i),Form("fHistMCPosKHypKaon%d",i),175,-3.5,3.5);      
422     fHistMCPosKHypProton[i] = new TH1F(Form("fHistMCPosKHypProton%d",i),Form("fHistMCPosKHypProton%d",i),175,-3.5,3.5);
423     fHistMCPosPHypPion[i] = new TH1F(Form("fHistMCPosPHypPion%d",i),Form("fHistMCPosPHypPion%d",i),175,-3.5,3.5);       
424     fHistMCPosPHypKaon[i] = new TH1F(Form("fHistMCPosPHypKaon%d",i),Form("fHistMCPosPHypKaon%d",i),175,-3.5,3.5);
425     fHistMCPosPHypProton[i]  = new TH1F(Form("fHistMCPosPHypProton%d",i),Form("fHistMCPosPHypProton%d",i),175,-3.5,3.5);        
426     
427     fHistMCNegOtherHypPion[i] = new TH1F(Form("fHistMCNegOtherHypPion%d",i),Form("fHistMCNegOtherHypPion%d",i),175,-3.5,3.5);   //MC truth
428     fHistMCNegOtherHypKaon[i] = new TH1F(Form("fHistMCNegOtherHypKaon%d",i),Form("fHistMCNegOtherHypKaon%d",i),175,-3.5,3.5);
429     fHistMCNegOtherHypProton[i] = new TH1F(Form("fHistMCNegOtherHypProton%d",i),Form("fHistMCNegOtherHypProton%d",i),175,-3.5,3.5);
430     fHistMCNegElHypPion[i] = new TH1F(Form("fHistMCNegElHypPion%d",i),Form("fHistMCNegElHypPion%d",i),175,-3.5,3.5);
431     fHistMCNegElHypKaon[i] = new TH1F(Form("fHistMCNegElHypKaon%d",i),Form("fHistMCNegElHypKaon%d",i),175,-3.5,3.5);
432     fHistMCNegElHypProton[i] = new TH1F(Form("fHistMCNegElHypProton%d",i),Form("fHistMCNegElHypProton%d",i),175,-3.5,3.5);
433     fHistMCNegPiHypPion[i] = new TH1F(Form("fHistMCNegPiHypPion%d",i),Form("fHistMCNegPiHypPion%d",i),175,-3.5,3.5);
434     fHistMCNegPiHypKaon[i] = new TH1F(Form("fHistMCNegPiHypKaon%d",i),Form("fHistMCNegPiHypKaon%d",i),175,-3.5,3.5);
435     fHistMCNegPiHypProton[i] = new TH1F(Form("fHistMCNegPiHypProton%d",i),Form("fHistMCNegPiHypProton%d",i),175,-3.5,3.5);
436     fHistMCNegKHypPion[i] = new TH1F(Form("fHistMCNegKHypPion%d",i),Form("fHistMCNegKHypPion%d",i),175,-3.5,3.5);       
437     fHistMCNegKHypKaon[i]  = new TH1F(Form("fHistMCNegKHypKaon%d",i),Form("fHistMCNegKHypKaon%d",i),175,-3.5,3.5);      
438     fHistMCNegKHypProton[i] = new TH1F(Form("fHistMCNegKHypProton%d",i),Form("fHistMCNegKHypProton%d",i),175,-3.5,3.5);
439     fHistMCNegPHypPion[i] = new TH1F(Form("fHistMCNegPHypPion%d",i),Form("fHistMCNegPHypPion%d",i),175,-3.5,3.5);       
440     fHistMCNegPHypKaon[i] = new TH1F(Form("fHistMCNegPHypKaon%d",i),Form("fHistMCNegPHypKaon%d",i),175,-3.5,3.5);
441     fHistMCNegPHypProton[i]  = new TH1F(Form("fHistMCNegPHypProton%d",i),Form("fHistMCNegPHypProton%d",i),175,-3.5,3.5);        
442     
443     
444     fOutput->Add(fHistPosPi[i]);
445     fOutput->Add(fHistPosK[i]);
446     fOutput->Add(fHistPosP[i]);
447     fOutput->Add(fHistNegPi[i]);
448     fOutput->Add(fHistNegK[i]);
449     fOutput->Add(fHistNegP[i]);
450     
451     fOutput->Add(fHistDCAPosPi[i]); //DCA distr
452     fOutput->Add(fHistDCAPosK[i]);
453     fOutput->Add(fHistDCAPosP[i]);
454     fOutput->Add(fHistDCANegPi[i]);
455     fOutput->Add(fHistDCANegK[i]);
456     fOutput->Add(fHistDCANegP[i]);
457     
458     fOutput->Add(fHistMCPrimDCAPosPi[i]);//DCA distr.
459     fOutput->Add(fHistMCPrimDCAPosK[i]);
460     fOutput->Add(fHistMCPrimDCAPosP[i]);
461     fOutput->Add(fHistMCPrimDCANegPi[i]);
462     fOutput->Add(fHistMCPrimDCANegK[i]);
463     fOutput->Add(fHistMCPrimDCANegP[i]);
464
465     fOutput->Add(fHistMCSecStDCAPosPi[i]);//DCA distr.
466     fOutput->Add(fHistMCSecStDCAPosK[i]);
467     fOutput->Add(fHistMCSecStDCAPosP[i]);
468     fOutput->Add(fHistMCSecStDCANegPi[i]);
469     fOutput->Add(fHistMCSecStDCANegK[i]);
470     fOutput->Add(fHistMCSecStDCANegP[i]);
471
472     fOutput->Add(fHistMCSecMatDCAPosPi[i]);//DCA distr.
473     fOutput->Add(fHistMCSecMatDCAPosK[i]);
474     fOutput->Add(fHistMCSecMatDCAPosP[i]);
475     fOutput->Add(fHistMCSecMatDCANegPi[i]);
476     fOutput->Add(fHistMCSecMatDCANegK[i]);
477     fOutput->Add(fHistMCSecMatDCANegP[i]);
478
479     fOutput->Add(fHistMCPosOtherHypPion[i]);//MC truth
480     fOutput->Add(fHistMCPosOtherHypKaon[i]);
481     fOutput->Add(fHistMCPosOtherHypProton[i]);
482     fOutput->Add(fHistMCPosElHypPion[i]);
483     fOutput->Add(fHistMCPosElHypKaon[i]);
484     fOutput->Add(fHistMCPosElHypProton[i]);
485     fOutput->Add(fHistMCPosPiHypPion[i]);
486     fOutput->Add(fHistMCPosPiHypKaon[i]);
487     fOutput->Add(fHistMCPosPiHypProton[i]);
488     fOutput->Add(fHistMCPosKHypPion[i]);
489     fOutput->Add(fHistMCPosKHypKaon[i]);
490     fOutput->Add(fHistMCPosKHypProton[i]);
491     fOutput->Add(fHistMCPosPHypPion[i]);
492     fOutput->Add(fHistMCPosPHypKaon[i]);
493     fOutput->Add(fHistMCPosPHypProton[i]);
494     
495     fOutput->Add(fHistMCNegOtherHypPion[i]);//MC truth
496     fOutput->Add(fHistMCNegOtherHypKaon[i]);
497     fOutput->Add(fHistMCNegOtherHypProton[i]);
498     fOutput->Add(fHistMCNegElHypPion[i]);
499     fOutput->Add(fHistMCNegElHypKaon[i]);
500     fOutput->Add(fHistMCNegElHypProton[i]);
501     fOutput->Add(fHistMCNegPiHypPion[i]);
502     fOutput->Add(fHistMCNegPiHypKaon[i]);
503     fOutput->Add(fHistMCNegPiHypProton[i]);
504     fOutput->Add(fHistMCNegKHypPion[i]);
505     fOutput->Add(fHistMCNegKHypKaon[i]);
506     fOutput->Add(fHistMCNegKHypProton[i]);
507     fOutput->Add(fHistMCNegPHypPion[i]);
508     fOutput->Add(fHistMCNegPHypKaon[i]);
509     fOutput->Add(fHistMCNegPHypProton[i]);
510     
511   }
512   
513   //NSigma Histos
514   for(Int_t j=0;j<3;j++){
515     
516     fHistPosNSigmaMean[j] = new TH1F(Form("hHistPosNSigmaMean%d",j),Form("hHistPosNSigmaMean%d",j),kNbins,fPtBinLimits);
517     fHistNegNSigmaMean[j] = new TH1F(Form("hHistNegNSigmaMean%d",j),Form("hHistNegNSigmaMean%d",j),kNbins,fPtBinLimits);
518     fHistPosNSigmaMCMean[j] = new TH1F(Form("hHistPosNSigmaMCMean%d",j),Form("hHistPosNSigmaMCMean%d",j),kNbins,fPtBinLimits);
519     fHistNegNSigmaMCMean[j] = new TH1F(Form("hHistNegNSigmaMCMean%d",j),Form("hHistNegNSigmaMCMean%d",j),kNbins,fPtBinLimits);
520     fHistPosNSigmaPrimMean[j] = new TH1F(Form("hHistPosNSigmaPrimMean%d",j),Form("hHistPosNSigmaPrimMean%d",j),kNbins,fPtBinLimits);
521     fHistNegNSigmaPrimMean[j] = new TH1F(Form("hHistNegNSigmaPrimMean%d",j),Form("hHistNegNSigmaPrimMean%d",j),kNbins,fPtBinLimits);
522     fHistPosNSigmaPrimMCMean[j] = new TH1F(Form("hHistPosNSigmaPrimMCMean%d",j),Form("hHistPosNSigmaPrimMCMean%d",j),kNbins,fPtBinLimits);
523     fHistNegNSigmaPrimMCMean[j] = new TH1F(Form("hHistNegNSigmaPrimMCMean%d",j),Form("hHistNegNSigmaPrimMCMean%d",j),kNbins,fPtBinLimits);
524     fOutput->Add(fHistPosNSigmaMean[j]);
525     fOutput->Add(fHistNegNSigmaMean[j]);
526     fOutput->Add(fHistPosNSigmaMCMean[j]);
527     fOutput->Add(fHistNegNSigmaMCMean[j]);
528     fOutput->Add(fHistPosNSigmaPrimMean[j]);
529     fOutput->Add(fHistNegNSigmaPrimMean[j]);
530     fOutput->Add(fHistPosNSigmaPrimMCMean[j]);
531     fOutput->Add(fHistNegNSigmaPrimMCMean[j]);
532
533     fHistPosNSigma[j] = new TH1F(Form("hHistPosNSigma%d",j),Form("hHistPosNSigma%d",j),kNbins,fPtBinLimits);
534     fHistNegNSigma[j] = new TH1F(Form("hHistNegNSigma%d",j),Form("hHistNegNSigma%d",j),kNbins,fPtBinLimits);
535     fHistPosNSigmaMC[j] = new TH1F(Form("hHistPosNSigmaMC%d",j),Form("hHistPosNSigmaMC%d",j),kNbins,fPtBinLimits);
536     fHistNegNSigmaMC[j] = new TH1F(Form("hHistNegNSigmaMC%d",j),Form("hHistNegNSigmaMC%d",j),kNbins,fPtBinLimits);
537     fHistPosNSigmaPrim[j] = new TH1F(Form("hHistPosNSigmaPrim%d",j),Form("hHistPosNSigmaPrim%d",j),kNbins,fPtBinLimits);
538     fHistNegNSigmaPrim[j] = new TH1F(Form("hHistNegNSigmaPrim%d",j),Form("hHistNegNSigmaPrim%d",j),kNbins,fPtBinLimits);
539     fHistPosNSigmaPrimMC[j] = new TH1F(Form("hHistPosNSigmaPrimMC%d",j),Form("hHistPosNSigmaPrimMC%d",j),kNbins,fPtBinLimits);
540     fHistNegNSigmaPrimMC[j] = new TH1F(Form("hHistNegNSigmaPrimMC%d",j),Form("hHistNegNSigmaPrimMC%d",j),kNbins,fPtBinLimits);
541     fOutput->Add(fHistPosNSigma[j]);
542     fOutput->Add(fHistNegNSigma[j]);
543     fOutput->Add(fHistPosNSigmaMC[j]);
544     fOutput->Add(fHistNegNSigmaMC[j]);
545     fOutput->Add(fHistPosNSigmaPrim[j]);
546     fOutput->Add(fHistNegNSigmaPrim[j]);
547     fOutput->Add(fHistPosNSigmaPrimMC[j]);
548     fOutput->Add(fHistNegNSigmaPrimMC[j]);
549
550   }
551   
552   fNtupleNSigma = new TNtuple("fNtupleNSigma","fNtupleNSigma","p:pt:dedx:dedx3:dedx4:dedx5:dedx6:ncls:nclspid:sign:run:eta:impactXY:impactZ:isph:pdgcode:mfl:uniqueID:chi2ncls");
553   fOutput->Add(fNtupleNSigma);
554   fNtupleMC = new TNtuple("fNtupleMC","fNtupleMC","ptMC:pdgcode:signMC:etaMC:yMC:isph:evSel:run");
555   fOutput->Add(fNtupleMC);
556
557   
558  
559   // Post output data.
560   PostData(1,fOutput);
561
562   Printf("end of CreateOutputObjects");
563 }
564
565 //________________________________________________________________________
566 void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
567   // Main loop
568   // Called for each event
569   
570   ///////////////////////////////////////
571   //variables
572   Float_t pdgmass[4]={0.13957,0.493677,0.938272,1.8756}; //mass for pi, K, P (Gev/c^2)
573   Int_t listcode[3]={211,321,2212};//code for pi, K, P (Gev/c^2)
574   Double_t dedxLay[4];
575   Float_t ptMC=-999;
576   Int_t code=-999, signMC=-999,isph=-999,mfl=-999,uniqueID=-999;
577   Float_t impactXY=-999, impactZ=-999;
578   Int_t evSel=1;
579   AliESDtrack* track;
580   UInt_t status; 
581   AliStack* stack=0;
582   TParticle *part=0;
583   TParticlePDG *pdgPart=0;
584         
585   /////////////////////
586
587   fESD=(AliESDEvent*)InputEvent();
588   if(!fESD) {
589     printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
590     return;
591   } 
592   fHistNEvents->Fill(-1);
593   
594   UInt_t maskPhysSel = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
595   TString firedTriggerClasses=fESD->GetFiredTriggerClasses();
596   //  if(!firedTriggerClasses.Contains("CINT1B")) return;
597   if((maskPhysSel & AliVEvent::kMB)==0) return;
598   fHistNEvents->Fill(0);
599
600   if(fLowEnergypp && !fMC){ // remove events without SDD in pp 2.76 TeV
601     if(!firedTriggerClasses.Contains("ALL")) return;
602   }
603   fHistNEvents->Fill(1);
604
605
606   if(fMC){
607     AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
608     if (!eventHandler) {
609       Printf("ERROR: Could not retrieve MC event handler");
610       return;
611     }
612     AliMCEvent* mcEvent = eventHandler->MCEvent();
613     if (!mcEvent) {
614       Printf("ERROR: Could not retrieve MC event");
615       return;
616     }
617     stack = mcEvent->Stack();
618     if (!stack) {
619       printf("ERROR: stack not available\n");
620       return;
621     }
622   }
623   if(!fITSPIDResponse){
624     fITSPIDResponse=new AliITSPIDResponse(fMC); 
625   }
626
627   //flags for MC
628   Int_t nTrackMC=0; 
629   if(stack) nTrackMC = stack->GetNtrack();      
630   const AliESDVertex *vtx =  fESD->GetPrimaryVertexSPD();
631
632   ///////////selection of the centrality or multiplicity bin
633
634   //selection on the event centrality
635   if(fHImode){
636     if(!(fLowCentrality<0.0)&&fUpCentrality>0.0)
637       {
638         AliCentrality *centrality = fESD->GetCentrality();
639         if(!centrality->IsEventInCentralityClass(fLowCentrality,fUpCentrality,"V0M"))
640           return;
641         Printf("Centrality of the event: %.1f",centrality->GetCentralityPercentile("V0M"));
642         Printf("Centrality cut: %.1f to %.1f",fLowCentrality,fUpCentrality);
643         fHistCen->Fill(centrality->GetCentralityPercentile("V0M"));
644       }
645   }
646   
647   //selection on the event multiplicity based on global tracks
648   Int_t multiplicity = -1;  
649   if(fMultEstimator==0){
650     // tracks+tracklets
651     multiplicity=AliESDtrackCuts::GetReferenceMultiplicity(fESD, AliESDtrackCuts::kTrackletsITSTPC, 0.8); 
652   }else if(fMultEstimator==1){
653     // tracklets
654     multiplicity=AliESDtrackCuts::GetReferenceMultiplicity(fESD, AliESDtrackCuts::kTracklets, 0.8);    
655   }else if(fMultEstimator==2){
656     // clusters in SPD1
657     const AliMultiplicity *mult = fESD->GetMultiplicity();
658     Float_t nClu1 = (Float_t)mult->GetNumberOfITSClusters(1);
659     multiplicity =(Int_t)(AliESDUtils::GetCorrSPD2(nClu1,vtx->GetZ())+0.5);
660   }
661   if(fLowMult>-1 && multiplicity<fLowMult) return;
662   if(fUpMult>-1 && multiplicity>fUpMult) return;
663   fHistMult->Fill(multiplicity);
664   fHistNEvents->Fill(2);
665   
666   //event selection
667   if(!vtx)evSel=0;
668   else{
669     fHistNEvents->Fill(3);
670     if(vtx->GetNContributors()<0) evSel=0;
671     else{
672       fHistNEvents->Fill(4);
673       if(TMath::Abs(vtx->GetZv())>10) evSel=0;
674       else{
675         fHistNEvents->Fill(5);
676         if(vtx->GetZRes()>0.5) evSel=0;
677         else{
678           fHistNEvents->Fill(6);
679           if(vtx->IsFromVertexerZ() && vtx->GetDispersion()>0.03) evSel=0;
680           else fHistNEvents->Fill(7);
681         }
682       }
683     }
684   }
685
686
687   /////first loop on stack, before event selection, filling MC ntuple
688   
689   for(Int_t imc=0; imc<nTrackMC; imc++){
690     part = stack->Particle(imc);
691     isph=1;    
692     if(!stack->IsPhysicalPrimary(imc)) isph=0;
693     pdgPart = part->GetPDG();
694     if(!pdgPart)continue;
695     if(pdgPart->Charge()==0) continue; //no neutral particles
696     Float_t yMC=-999.;
697     if(part->Energy() != TMath::Abs(part->Pz())) yMC = 0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
698     if(TMath::Abs(yMC) > fMaxY) continue; //rapidity cut
699     if(pdgPart->Charge()>0) signMC=1;
700     else signMC=-1;
701     ptMC=part->Pt();
702     code=pdgPart->PdgCode();
703     Int_t jpart=-1;
704     for(Int_t j=0; j<3; j++){
705       if(TMath::Abs(code)==listcode[j]){
706         jpart=j;
707         break;
708       }
709     }
710     Int_t indexMoth=part->GetFirstMother();
711     if(indexMoth>=0){
712       TParticle* moth = stack->Particle(indexMoth);
713       Float_t codemoth = TMath::Abs(moth->GetPdgCode());
714       mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
715     }
716     uniqueID = part->GetUniqueID();
717     
718     //filling MC ntuple
719     if(TMath::Abs(code)==211 || TMath::Abs(code)==321 || TMath::Abs(code)==2212){
720       Float_t xntMC[8];
721       Int_t indexMC=0;
722       xntMC[indexMC++]=(Float_t)ptMC;
723       xntMC[indexMC++]=(Float_t)code;
724       xntMC[indexMC++]=(Float_t)signMC;
725       xntMC[indexMC++]=(Float_t)part->Eta();
726       xntMC[indexMC++]=(Float_t)yMC;
727       xntMC[indexMC++]=(Float_t)isph;
728       xntMC[indexMC++]=(Float_t)evSel;
729       xntMC[indexMC++]=(Float_t)fESD->GetRunNumber();
730       
731       if(fFillNtuple) fNtupleMC->Fill(xntMC);
732     }
733     
734     if(jpart>=0){
735       if(stack->IsPhysicalPrimary(imc)){
736         if(signMC>0) fHistPrimMCposBefEvSel[jpart]->Fill(ptMC);
737         else  fHistPrimMCnegBefEvSel[jpart]->Fill(ptMC);
738         if(evSel==1){
739           if(signMC>0) fHistPrimMCpos[jpart]->Fill(ptMC);
740           else  fHistPrimMCneg[jpart]->Fill(ptMC);
741         }
742       }else{
743         if(mfl==3 && uniqueID == kPDecay){ // If a particle is not a physical primary, check if it comes from weak decay
744           if(signMC>0) fHistSecStrMCposBefEvSel[jpart]->Fill(ptMC);
745           else  fHistSecStrMCnegBefEvSel[jpart]->Fill(ptMC);        
746           if(evSel==1){
747             if(signMC>0) fHistSecStrMCpos[jpart]->Fill(ptMC);
748             else  fHistSecStrMCneg[jpart]->Fill(ptMC);      
749           }
750         }else{
751           if(signMC>0) fHistSecMatMCposBefEvSel[jpart]->Fill(ptMC);
752           else  fHistSecMatMCnegBefEvSel[jpart]->Fill(ptMC);        
753           if(evSel==1){
754             if(signMC>0) fHistSecMatMCpos[jpart]->Fill(ptMC);
755             else  fHistSecMatMCneg[jpart]->Fill(ptMC);      
756           }
757         }
758       }
759     }
760   }     
761   
762   
763   if(evSel==0)return;  //event selection
764   
765
766   //loop on tracks
767   for (Int_t iTrack=0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {  
768     isph=-999;
769     code=-999;
770     mfl=-999;
771     uniqueID=-999;
772     
773     track = (AliESDtrack*)fESD->GetTrack(iTrack);      
774     if (!track) continue;
775     
776     //track selection
777     Int_t countBinTrk=1;
778     TString label;
779     
780     label="no selection";
781     fHistNTracks->Fill(countBinTrk);
782     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
783     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
784     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
785     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
786     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
787     countBinTrk++;  
788     
789     status=track->GetStatus();
790     if((status&AliESDtrack::kITSpureSA)==0) continue; //its standalone
791   
792     label="ITSsa";
793     fHistNTracks->Fill(countBinTrk);
794     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
795     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
796     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
797     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
798     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
799     countBinTrk++;    
800     
801     if((status&AliESDtrack::kITSrefit)==0) continue; //its refit
802     
803     label="ITSrefit";
804     fHistNTracks->Fill(countBinTrk);
805     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
806     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
807     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
808     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
809     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
810     countBinTrk++;
811     
812     if(TMath::Abs(track->GetSign())<0.0001) continue; //no neutral particles
813     
814     label="neutral particle";
815     fHistNTracks->Fill(countBinTrk);
816     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
817     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
818     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
819     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
820     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
821     countBinTrk++;  
822     
823     //cluster in ITS
824     UInt_t clumap = track->GetITSClusterMap();
825     Int_t nSPD=0;
826     for(Int_t il=0; il<2; il++) if(TESTBIT(clumap,il)) nSPD++;
827     if(nSPD<fMinSPDPts) continue;
828     
829     label="SPDcls";
830     fHistNTracks->Fill(countBinTrk);
831     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
832     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
833     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
834     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
835     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
836     countBinTrk++;  
837     
838     Int_t nPtsForPid=0;
839     for(Int_t j=2;j<6;j++) if(TESTBIT(clumap,j)) nPtsForPid++;
840     if(nPtsForPid<fMinNdEdxSamples) continue; //at least 3 points on SSD/SDD
841     
842     label="SDD+SSD cls";
843     fHistNTracks->Fill(countBinTrk);
844     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
845     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
846     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
847     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
848     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
849     countBinTrk++;  
850     
851     //chisquare/nclusters       
852     Int_t nclu=nSPD+nPtsForPid;
853     if(track->GetITSchi2()/nclu > fMaxChi2Clu) continue; 
854     
855     label="chi2/ncls";
856     fHistNTracks->Fill(countBinTrk);
857     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
858     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
859     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
860     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
861     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
862     countBinTrk++;  
863     
864     //pseudorapidity and rapidity
865     if(TMath::Abs(track->Eta()) > fEtaRange) continue;
866     
867     label="eta";
868     fHistNTracks->Fill(countBinTrk);
869     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
870     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
871     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
872     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
873     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
874     countBinTrk++;  
875     
876     //truncated mean
877     track->GetITSdEdxSamples(dedxLay);
878     Double_t dedx = CookdEdx(dedxLay);
879     if(dedx<0) continue;
880
881     label="de/dx<0";
882     fHistNTracks->Fill(countBinTrk);
883     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
884     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
885     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
886     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
887     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
888     countBinTrk++;  
889     
890     Float_t pt = track->Pt();
891     Int_t theBin=-1;
892     for(Int_t m=0; m<kNbins; m++){
893       if(TMath::Abs(pt) > fPtBinLimits[m] && TMath::Abs(pt) < fPtBinLimits[m+1]){
894         theBin=m;
895         break;
896       }
897     }
898     track->GetImpactParameters(impactXY, impactZ);
899     
900     //Filling Ntuple
901     //information from the MC kinematics
902     if(fMC){
903       if(track->GetLabel()<0)isph=-1;
904       if(track->GetLabel()>=0){
905         part = (TParticle*)stack->Particle(track->GetLabel());
906         pdgPart = part->GetPDG();
907         code = pdgPart->PdgCode();
908         if(stack->IsPhysicalPrimary(track->GetLabel())) isph=1;
909         else{ 
910           isph=0;
911           Int_t indexMoth=part->GetFirstMother();
912           if(indexMoth>=0){
913             TParticle* moth = stack->Particle(indexMoth);
914             Float_t codemoth = TMath::Abs(moth->GetPdgCode());
915             mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
916           }
917           uniqueID = part->GetUniqueID();
918         }
919         
920         //Filling DCA distribution with MC truth
921         
922         if(theBin>=0 && theBin<kNbins){
923           if(isph==1){//primaries in MC
924             if(track->GetSign()>0){
925               if(TMath::Abs(code)==listcode[0]) fHistMCPrimDCAPosPi[theBin]->Fill(impactXY);
926               if(TMath::Abs(code)==listcode[1]) fHistMCPrimDCAPosK[theBin]->Fill(impactXY);
927               if(TMath::Abs(code)==listcode[2]) fHistMCPrimDCAPosP[theBin]->Fill(impactXY);
928             }else{
929               if(TMath::Abs(code)==listcode[0]) fHistMCPrimDCANegPi[theBin]->Fill(impactXY);
930               if(TMath::Abs(code)==listcode[1]) fHistMCPrimDCANegK[theBin]->Fill(impactXY);
931               if(TMath::Abs(code)==listcode[2]) fHistMCPrimDCANegP[theBin]->Fill(impactXY);
932             }
933           }
934           
935           if(isph==0){//primaries in MC
936             if(mfl==3 && uniqueID == kPDecay){
937               if(track->GetSign()>0){
938                 if(TMath::Abs(code)==listcode[0]) fHistMCSecStDCAPosPi[theBin]->Fill(impactXY);
939                 if(TMath::Abs(code)==listcode[1]) fHistMCSecStDCAPosK[theBin]->Fill(impactXY);
940                 if(TMath::Abs(code)==listcode[2]) fHistMCSecStDCAPosP[theBin]->Fill(impactXY);
941               }else{
942                 if(TMath::Abs(code)==listcode[0]) fHistMCSecStDCANegPi[theBin]->Fill(impactXY);
943                 if(TMath::Abs(code)==listcode[1]) fHistMCSecStDCANegK[theBin]->Fill(impactXY);
944                 if(TMath::Abs(code)==listcode[2]) fHistMCSecStDCANegP[theBin]->Fill(impactXY);
945               }
946             }else{
947               if(track->GetSign()>0){
948                 if(TMath::Abs(code)==listcode[0]) fHistMCSecMatDCAPosPi[theBin]->Fill(impactXY);
949                 if(TMath::Abs(code)==listcode[1]) fHistMCSecMatDCAPosK[theBin]->Fill(impactXY);
950                 if(TMath::Abs(code)==listcode[2]) fHistMCSecMatDCAPosP[theBin]->Fill(impactXY);
951               }else{
952                 if(TMath::Abs(code)==listcode[0]) fHistMCSecMatDCANegPi[theBin]->Fill(impactXY);
953                 if(TMath::Abs(code)==listcode[1]) fHistMCSecMatDCANegK[theBin]->Fill(impactXY);
954                 if(TMath::Abs(code)==listcode[2]) fHistMCSecMatDCANegP[theBin]->Fill(impactXY);
955               }
956             }
957           }
958         }
959       }
960     }
961     Float_t xnt[19];
962     Int_t index=0;
963     xnt[index++]=(Float_t)track->GetP();
964     xnt[index++]=(Float_t)track->Pt();
965     xnt[index++]=(Float_t)dedx;
966     xnt[index++]=(Float_t)dedxLay[0];
967     xnt[index++]=(Float_t)dedxLay[1];
968     xnt[index++]=(Float_t)dedxLay[2];
969     xnt[index++]=(Float_t)dedxLay[3];
970     xnt[index++]=(Float_t)nclu;
971     xnt[index++]=(Float_t)nPtsForPid;
972     xnt[index++]=(Float_t)track->GetSign();
973     xnt[index++]=(Float_t)fESD->GetRunNumber();
974     xnt[index++]=(Float_t)track->Eta();
975     xnt[index++]=(Float_t)impactXY;
976     xnt[index++]=(Float_t)impactZ;
977     xnt[index++]=(Float_t)isph;
978     xnt[index++]=(Float_t)code;
979     xnt[index++]=(Float_t)mfl;
980     xnt[index++]=(Float_t)uniqueID;
981     xnt[index]=(Float_t)track->GetITSchi2()/nclu;
982           
983     if(fFillNtuple) fNtupleNSigma->Fill(xnt);
984     
985     
986         
987     //Compute y and bb
988     Double_t y[4],bbtheo[4],logdiff[4];
989     Float_t p=track->GetP();
990     if(fMC && fSmearMC){
991       dedx=fRandGener->Gaus(dedx,fSmeardEdx*dedx);
992       p=fRandGener->Gaus(p,fSmearP*p);     
993     }
994     
995     //Nsigma Method
996     Float_t resodedx[4];
997     for(Int_t ires=0;ires<4;ires++){
998       resodedx[ires]=fITSPIDResponse->GetResolution(1,ires+1,kTRUE);
999     }
1000     
1001     for(Int_t i=0;i<4;i++){
1002       y[i] = Eta2y(pt,pdgmass[i],track->Eta());
1003       //bbtheo[i]=fITSPIDResponse->Bethe(p,pdgmass[i],kTRUE); //Pure PHOBOS BB
1004       bbtheo[i]=fITSPIDResponse->BetheITSsaHybrid(p,pdgmass[i]); //PHOBOS + polinomial at low pt (below beta*gamma = 0.76)
1005       logdiff[i]=TMath::Log(dedx) - TMath::Log(bbtheo[i]);
1006     }
1007     
1008     Int_t resocls=(Int_t)nPtsForPid-1;
1009     
1010     //NSigma Method, with asymmetric bands
1011     Int_t minPosMean=-1;
1012     for(Int_t isp=0; isp<3; isp++){
1013       if(dedx<bbtheo[0])continue;
1014       Double_t bb=bbtheo[isp];
1015       if(dedx<bb){
1016         Double_t bbdistance=TMath::Abs((bbtheo[isp]-bbtheo[isp-1])/2);
1017         Double_t nsigma=TMath::Abs((dedx-bb)/bbdistance);
1018         if(nsigma<1.)minPosMean=isp;
1019       }
1020       else{
1021         Double_t bbdistance=TMath::Abs((bbtheo[isp]-bbtheo[isp+1])/2);
1022         Double_t nsigma=TMath::Abs((dedx-bb)/bbdistance);
1023         if(nsigma<1.)minPosMean=isp;
1024       }
1025     }
1026     if(dedx<bbtheo[0] && TMath::Abs((dedx-bbtheo[0])/(resodedx[resocls]*bbtheo[0]))<2)minPosMean=0;
1027     
1028     //NSigma method with simmetric bands
1029     
1030     Double_t nsigmas[3];
1031     Double_t min=999999.;
1032     Int_t minPos=-1;
1033     for(Int_t isp=0; isp<3; isp++){
1034       Double_t bb=bbtheo[isp];
1035       nsigmas[isp]=TMath::Abs((dedx-bb)/(resodedx[resocls]*bb));
1036       if(nsigmas[isp]<min){
1037         min=nsigmas[isp];
1038         minPos=isp;
1039       }
1040       //Filling histos with nsigma separation
1041       if(track->GetSign()>0)fHistPosNSigmaSep[isp]->Fill(track->GetP(),((dedx-bb)/(resodedx[resocls]*bb)));
1042       else fHistNegNSigmaSep[isp]->Fill(track->GetP(),((dedx-bb)/(resodedx[resocls]*bb)));
1043     }
1044     
1045     // y calculation
1046     Double_t yPartMean=-999;
1047     Double_t yPart=-999;
1048     if(minPosMean>-1) yPartMean=y[minPosMean];
1049     if(minPos>-1) yPart=y[minPos];
1050     
1051     if(TMath::Abs(yPartMean)<fMaxY){     
1052       //DCA distributions, before the DCA cuts, based on asymmetrinc nsigma approach
1053       if(theBin>=0 && theBin<kNbins){
1054         if(track->GetSign()>0){
1055           if(minPosMean==0) fHistDCAPosPi[theBin]->Fill(impactXY);
1056           else if(minPosMean==1) fHistDCAPosK[theBin]->Fill(impactXY);
1057           else if(minPosMean==2) fHistDCAPosP[theBin]->Fill(impactXY);
1058         }else{
1059           if(minPosMean==0) fHistDCANegPi[theBin]->Fill(impactXY);
1060           else if(minPosMean==1) fHistDCANegK[theBin]->Fill(impactXY);
1061           else if(minPosMean==2) fHistDCANegP[theBin]->Fill(impactXY);
1062         }
1063       } 
1064     }
1065     
1066     //DCA cut on xy and z
1067     if(!DCAcut(impactXY,impactZ,pt)) continue;
1068     
1069     label="DCA";
1070     fHistNTracks->Fill(countBinTrk);
1071     if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
1072     if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
1073     fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
1074     fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
1075     fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
1076     countBinTrk++;  
1077
1078
1079     Int_t jpart=-1;
1080     
1081     //Filling Histos for Reco Efficiency
1082     //information from the MC kinematics
1083     if(fMC){
1084       if(track->GetLabel()<0)isph=-1;
1085       if(track->GetLabel()>=0){
1086         part = (TParticle*)stack->Particle(track->GetLabel());
1087         pdgPart = part->GetPDG();
1088         code = pdgPart->PdgCode();
1089         for(Int_t j=0; j<3; j++){
1090           if(TMath::Abs(code)==listcode[j]){
1091             jpart=j;
1092             break;
1093           }
1094         }
1095         if(jpart>=0){
1096           if(pdgPart->Charge()>0) signMC=1;
1097           else signMC=-1;
1098           ptMC=part->Pt();
1099           if(stack->IsPhysicalPrimary(track->GetLabel())){
1100             if(signMC>0) fHistPrimMCposReco[jpart]->Fill(ptMC);
1101             else  fHistPrimMCnegReco[jpart]->Fill(ptMC);
1102           }else{ 
1103             Int_t indexMoth=part->GetFirstMother();
1104             if(indexMoth>=0){
1105               TParticle* moth = stack->Particle(indexMoth);
1106               Float_t codemoth = TMath::Abs(moth->GetPdgCode());
1107               mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
1108             }
1109             uniqueID = part->GetUniqueID();
1110             if(mfl==3 && uniqueID == kPDecay){ // strangeness
1111               if(signMC>0) fHistSecStrMCposReco[jpart]->Fill(ptMC);
1112               else  fHistSecStrMCnegReco[jpart]->Fill(ptMC);        
1113             }else{
1114               if(signMC>0) fHistSecMatMCposReco[jpart]->Fill(ptMC);
1115               else  fHistSecMatMCnegReco[jpart]->Fill(ptMC);        
1116             }
1117           }
1118         }
1119       }
1120     }
1121     
1122     //Nsigma histos with MC truth
1123     
1124     //asymmetric approach
1125     if(TMath::Abs(yPartMean)<fMaxY && minPosMean>-1){
1126       //nsigma histos
1127       if(track->GetSign()>0) fHistPosNSigmaMean[minPosMean]->Fill(pt);
1128       else fHistNegNSigmaMean[minPosMean]->Fill(pt);
1129       if(fMC){
1130         //nsigma histos with MC truth on PID
1131         if(TMath::Abs(code)==listcode[minPosMean]){
1132           if(track->GetSign()>0) fHistPosNSigmaMCMean[minPosMean]->Fill(pt);
1133           else fHistNegNSigmaMCMean[minPosMean]->Fill(pt);
1134         }
1135         //nsigma histos with MC truth on IsPhysicalPrimary
1136         if(isph==1){
1137           if(track->GetSign()>0) fHistPosNSigmaPrimMean[minPosMean]->Fill(pt);
1138           else fHistNegNSigmaPrimMean[minPosMean]->Fill(pt);
1139           //nsigma histos with MC truth on IsPhysicalPrimary and PID
1140           if(TMath::Abs(code)==listcode[minPosMean]){
1141             if(track->GetSign()>0) fHistPosNSigmaPrimMCMean[minPosMean]->Fill(pt);
1142             else fHistNegNSigmaPrimMCMean[minPosMean]->Fill(pt);
1143           }
1144         }
1145       }
1146     }
1147     
1148     //symmetric bands
1149     if(min<fMinNSigma && TMath::Abs(yPart)<fMaxY){
1150       //nsigma histos
1151       if(track->GetSign()>0) fHistPosNSigma[minPos]->Fill(pt);
1152       else fHistNegNSigma[minPos]->Fill(pt);
1153       if(fMC){
1154         //nsigma histos with MC truth on PID
1155         if(TMath::Abs(code)==listcode[minPos]){
1156           if(track->GetSign()>0) fHistPosNSigmaMC[minPos]->Fill(pt);
1157           else fHistNegNSigmaMC[minPos]->Fill(pt);
1158         }
1159         //nsigma histos with MC truth on IsPhysicalPrimary
1160         if(isph==1){
1161           if(track->GetSign()>0) fHistPosNSigmaPrim[minPos]->Fill(pt);
1162           else fHistNegNSigmaPrim[minPos]->Fill(pt);
1163           //nsigma histos with MC truth on IsPhysicalPrimary and PID
1164           if(TMath::Abs(code)==listcode[minPos]){
1165             if(track->GetSign()>0) fHistPosNSigmaPrimMC[minPos]->Fill(pt);
1166             else fHistNegNSigmaPrimMC[minPos]->Fill(pt);
1167           }
1168         }
1169       }
1170     }
1171     
1172     
1173     //integral approach histograms
1174     if(theBin>=0 && theBin<kNbins){
1175       if(track->GetSign()>0){
1176         if(TMath::Abs(y[0]) < fMaxY)fHistPosPi[theBin]->Fill(logdiff[0]);
1177         if(TMath::Abs(y[1]) < fMaxY)fHistPosK[theBin]->Fill(logdiff[1]);
1178         if(TMath::Abs(y[2]) < fMaxY)fHistPosP[theBin]->Fill(logdiff[2]);
1179         if(fMC){
1180           if(TMath::Abs(y[0])<fMaxY){
1181             if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypPion[theBin]->Fill(logdiff[0]);
1182             if(TMath::Abs(code)==11)fHistMCPosElHypPion[theBin]->Fill(logdiff[0]);
1183             if(TMath::Abs(code)==211)fHistMCPosPiHypPion[theBin]->Fill(logdiff[0]);
1184             if(TMath::Abs(code)==321)fHistMCPosKHypPion[theBin]->Fill(logdiff[0]);
1185             if(TMath::Abs(code)==2212)fHistMCPosPHypPion[theBin]->Fill(logdiff[0]);
1186           }
1187           if(TMath::Abs(y[1])<fMaxY){
1188             if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypKaon[theBin]->Fill(logdiff[1]);
1189             if(TMath::Abs(code)==11)fHistMCPosElHypKaon[theBin]->Fill(logdiff[1]);
1190             if(TMath::Abs(code)==211)fHistMCPosPiHypKaon[theBin]->Fill(logdiff[1]);
1191             if(TMath::Abs(code)==321)fHistMCPosKHypKaon[theBin]->Fill(logdiff[1]);
1192             if(TMath::Abs(code)==2212)fHistMCPosPHypKaon[theBin]->Fill(logdiff[1]);
1193           }
1194           if(TMath::Abs(y[2])<fMaxY){
1195               if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypProton[theBin]->Fill(logdiff[2]);
1196               if(TMath::Abs(code)==11)fHistMCPosElHypProton[theBin]->Fill(logdiff[2]);
1197               if(TMath::Abs(code)==211)fHistMCPosPiHypProton[theBin]->Fill(logdiff[2]);
1198               if(TMath::Abs(code)==321)fHistMCPosKHypProton[theBin]->Fill(logdiff[2]);
1199               if(TMath::Abs(code)==2212)fHistMCPosPHypProton[theBin]->Fill(logdiff[2]);
1200           }
1201         }
1202       }else{
1203         if(TMath::Abs(y[0]) < fMaxY)fHistNegPi[theBin]->Fill(logdiff[0]);
1204         if(TMath::Abs(y[1]) < fMaxY)fHistNegK[theBin]->Fill(logdiff[1]);
1205         if(TMath::Abs(y[2]) < fMaxY)fHistNegP[theBin]->Fill(logdiff[2]);
1206         if(fMC){
1207           if(TMath::Abs(y[0])<fMaxY){
1208             if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypPion[theBin]->Fill(logdiff[0]);
1209             if(TMath::Abs(code)==11)fHistMCNegElHypPion[theBin]->Fill(logdiff[0]);
1210             if(TMath::Abs(code)==211)fHistMCNegPiHypPion[theBin]->Fill(logdiff[0]);
1211             if(TMath::Abs(code)==321)fHistMCNegKHypPion[theBin]->Fill(logdiff[0]);
1212             if(TMath::Abs(code)==2212)fHistMCNegPHypPion[theBin]->Fill(logdiff[0]);
1213           }
1214           if(TMath::Abs(y[1])<fMaxY){
1215             if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypKaon[theBin]->Fill(logdiff[1]);
1216             if(TMath::Abs(code)==11)fHistMCNegElHypKaon[theBin]->Fill(logdiff[1]);
1217             if(TMath::Abs(code)==211)fHistMCNegPiHypKaon[theBin]->Fill(logdiff[1]);
1218             if(TMath::Abs(code)==321)fHistMCNegKHypKaon[theBin]->Fill(logdiff[1]);
1219             if(TMath::Abs(code)==2212)fHistMCNegPHypKaon[theBin]->Fill(logdiff[1]);
1220           }
1221           if(TMath::Abs(y[2])<fMaxY){
1222             if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypProton[theBin]->Fill(logdiff[2]);
1223             if(TMath::Abs(code)==11)fHistMCNegElHypProton[theBin]->Fill(logdiff[2]);
1224             if(TMath::Abs(code)==211)fHistMCNegPiHypProton[theBin]->Fill(logdiff[2]);
1225             if(TMath::Abs(code)==321)fHistMCNegKHypProton[theBin]->Fill(logdiff[2]);
1226             if(TMath::Abs(code)==2212)fHistMCNegPHypProton[theBin]->Fill(logdiff[2]);
1227           }
1228         }
1229       }
1230     }                                                        
1231   
1232           
1233     //fill propaganda plot with dedx
1234     fHistDEDX->Fill(track->GetP(),dedx);
1235     fHistDEDXdouble->Fill(track->GetP()*track->GetSign(),dedx);
1236           
1237     //fill charge distribution histo to check the calibration
1238     for(Int_t j=0;j<4;j++){
1239       if(dedxLay[j]<5) continue;
1240       fHistCharge[j]->Fill(dedxLay[j]);
1241     }
1242   }
1243                 
1244   // Post output data.
1245   PostData(1,fOutput);
1246   PostData(2,fListCuts);
1247   //  Printf("............. end of Exec");
1248 }      
1249
1250 //________________________________________________________________________
1251 void AliAnalysisTaskSEITSsaSpectra::Terminate(Option_t *) {
1252   // Merge output
1253   // Called once at the end of the query
1254   
1255   fOutput = dynamic_cast<TList*>(GetOutputData(1));
1256   if (!fOutput) {
1257     printf("ERROR: fOutput not available\n");
1258     return;
1259   } 
1260   fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
1261   printf("Number of Analyzed Events = %f\n",fHistNEvents->GetBinContent(1));
1262   
1263   Printf("end of Terminate");
1264   return;
1265 }