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