Added protection for unidentified tracks (Nicolo' Jacazio)
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / ITSsa / AliAnalysisTaskSEITSsaSpectra.cxx
CommitLineData
36be14b3 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
36be14b3 27#include <TH1F.h>
a1c5e4a1 28#include <TF1.h>
36be14b3 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"
c841b89f 46#include "AliESDtrackCuts.h"
1f75ac7d 47#include "AliCentrality.h"
48#include "AliMultiplicity.h"
49#include "AliESDUtils.h"
98bb3219 50#include "AliITSPIDResponse.h"
36be14b3 51
52ClassImp(AliAnalysisTaskSEITSsaSpectra)
6b77d2c0 53/* $Id$ */
36be14b3 54
55//________________________________________________________________________
56AliAnalysisTaskSEITSsaSpectra::AliAnalysisTaskSEITSsaSpectra():
a1c5e4a1 57AliAnalysisTaskSE("taskITSsaSpectra"),
36be14b3 58 fESD(0),
59 fOutput(0),
a1c5e4a1 60 fListCuts(0),
36be14b3 61 fHistNEvents(0),
c841b89f 62 fHistMult(0),
1f75ac7d 63 fHistCen(0),
36be14b3 64 fHistNTracks(0),
c841b89f 65 fHistNTracksPos(0),
66 fHistNTracksNeg(0),
36be14b3 67 fHistDEDX(0),
68 fHistDEDXdouble(0),
69 fHistBeforeEvSel(0),
70 fHistAfterEvSel(0),
a1c5e4a1 71 fDCAxyCutFunc(0x0),
72 fDCAzCutFunc(0x0),
98bb3219 73 fITSPIDResponse(0),
36be14b3 74 fMinSPDPts(1),
75 fMinNdEdxSamples(3),
76 fMindEdx(0.),
c841b89f 77 fMinNSigma(1.5),
36be14b3 78 fMaxY(0.5),
c841b89f 79 fMaxChi2Clu(2.5),
36be14b3 80 fNSigmaDCAxy(7.),
81 fNSigmaDCAz(7.),
dac69cb0 82 fEtaRange(0.8),
1f75ac7d 83 fLowMult(-1),
84 fUpMult(-1),
85 fLowCentrality(-1.0),
86 fUpCentrality(-1.0),
a1c5e4a1 87 fMultEstimator(0),
1f75ac7d 88 fHImode(0),
6b77d2c0 89 fYear(2010),
36be14b3 90 fMC(kFALSE),
91 fSmearMC(kFALSE),
92 fSmearP(0.),
93 fSmeardEdx(0.),
94 fRandGener(0),
95 fFillNtuple(kFALSE),
146d18ee 96 fLowEnergypp(kFALSE),
36be14b3 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);
36be14b3 104 DefineInput(0, TChain::Class());
105 DefineOutput(1, TList::Class());
a1c5e4a1 106 DefineOutput(2,TList::Class());
36be14b3 107 Printf("end of AliAnalysisTaskSEITSsaSpectra");
108}
109
110//___________________________________________________________________________
111AliAnalysisTaskSEITSsaSpectra::~AliAnalysisTaskSEITSsaSpectra(){
112 // Destructor
113 if (fOutput) {
114 delete fOutput;
115 fOutput = 0;
116 }
a1c5e4a1 117 if (fListCuts) {
118 delete fListCuts;
119 fListCuts = 0;
120 }
36be14b3 121 if(fRandGener) delete fRandGener;
98bb3219 122 if(fITSPIDResponse) delete fITSPIDResponse;
a1c5e4a1 123 delete fDCAxyCutFunc;
124 delete fDCAzCutFunc;
36be14b3 125}
126
127//________________________________________________________________________
4d668cf7 128Double_t AliAnalysisTaskSEITSsaSpectra::CookdEdx(Double_t *s) const {
36be14b3 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
7f8f788e 165//________________________________________________________________________
a1c5e4a1 166Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt) const {
36be14b3 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
a1c5e4a1 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//________________________________________________________________________
180Double_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//________________________________________________________________________
188void AliAnalysisTaskSEITSsaSpectra::Init(){
189 //
190 // Initialization
191 //
192 fListCuts=new TList();
193 fListCuts->SetOwner();
36be14b3 194 Double_t xyP[3];
195 Double_t zP[3];
a1c5e4a1 196 if(fMC){
6b77d2c0 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 }
a1c5e4a1 212 }else{
6b77d2c0 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 }
36be14b3 228 }
a1c5e4a1 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);
36be14b3 232
36be14b3 233
a1c5e4a1 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}
36be14b3 244
36be14b3 245//________________________________________________________________________
246void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){
4d668cf7 247 // Create a TList with histograms and a TNtuple
248 // Called once
029e1796 249
36be14b3 250 fOutput = new TList();
251 fOutput->SetOwner();
252 fOutput->SetName("Spiderman");
253
a1c5e4a1 254 fHistNEvents = new TH1F("fHistNEvents", "Number of processed events",9,-1.5,7.5);
36be14b3 255 fHistNEvents->Sumw2();
256 fHistNEvents->SetMinimum(0);
a1c5e4a1 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");
36be14b3 266 fOutput->Add(fHistNEvents);
267
1f75ac7d 268 fHistMult = new TH1F("fHistMult", "Event Multiplicity",3000,-0.5,2999.5);
c841b89f 269 fHistMult->Sumw2();
270 fHistMult->SetMinimum(0);
a1c5e4a1 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");
c841b89f 274 fOutput->Add(fHistMult);
275
1f75ac7d 276 fHistCen = new TH1F("fHistCen", "Event Centrality",101,-0.5,100.5);
277 fHistCen->Sumw2();
278 fHistCen->SetMinimum(0);
279 fOutput->Add(fHistCen);
280
36be14b3 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
c841b89f 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
36be14b3 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
36be14b3 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++){
98bb3219 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
10ba520e 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]);
36be14b3 367
36be14b3 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
10ba520e 383 fHistDCAPosPi[i] = new TH1F(Form("fHistDCAPosPi%d",i),Form("fHistDCAPosPi%d",i),2000,-1,1); //DCA distr. with NSigma PID
36be14b3 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
10ba520e 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
16c730d4 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);
e6076bb8 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
16c730d4 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);
e6076bb8 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
36be14b3 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
6b77d2c0 451 fOutput->Add(fHistDCAPosPi[i]); //DCA distr
36be14b3 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
10ba520e 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
16c730d4 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]);
e6076bb8 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
16c730d4 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]);
e6076bb8 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
36be14b3 511 }
512
513 //NSigma Histos
514 for(Int_t j=0;j<3;j++){
10ba520e 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
36be14b3 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);
10ba520e 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);
36be14b3 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]);
10ba520e 543 fOutput->Add(fHistPosNSigmaMC[j]);
544 fOutput->Add(fHistNegNSigmaMC[j]);
36be14b3 545 fOutput->Add(fHistPosNSigmaPrim[j]);
546 fOutput->Add(fHistNegNSigmaPrim[j]);
547 fOutput->Add(fHistPosNSigmaPrimMC[j]);
548 fOutput->Add(fHistNegNSigmaPrimMC[j]);
10ba520e 549
36be14b3 550 }
551
7c831837 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");
36be14b3 553 fOutput->Add(fNtupleNSigma);
554 fNtupleMC = new TNtuple("fNtupleMC","fNtupleMC","ptMC:pdgcode:signMC:etaMC:yMC:isph:evSel:run");
555 fOutput->Add(fNtupleMC);
7f8f788e 556
a1c5e4a1 557
558
7f8f788e 559 // Post output data.
560 PostData(1,fOutput);
561
36be14b3 562 Printf("end of CreateOutputObjects");
563}
564
565//________________________________________________________________________
566void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
4d668cf7 567 // Main loop
568 // Called for each event
36be14b3 569
1f75ac7d 570 ///////////////////////////////////////
36be14b3 571 //variables
10ba520e 572 Float_t pdgmass[4]={0.13957,0.493677,0.938272,1.8756}; //mass for pi, K, P (Gev/c^2)
36be14b3 573 Int_t listcode[3]={211,321,2212};//code for pi, K, P (Gev/c^2)
146d18ee 574 Double_t dedxLay[4];
5c112001 575 Float_t ptMC=-999;
7c831837 576 Int_t code=-999, signMC=-999,isph=-999,mfl=-999,uniqueID=-999;
36be14b3 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
36be14b3 585 /////////////////////
146d18ee 586
587 fESD=(AliESDEvent*)InputEvent();
588 if(!fESD) {
589 printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
590 return;
591 }
592 fHistNEvents->Fill(-1);
593
a1c5e4a1 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;
146d18ee 598 fHistNEvents->Fill(0);
599
a1c5e4a1 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
146d18ee 605
36be14b3 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 }
98bb3219 623 if(!fITSPIDResponse){
624 fITSPIDResponse=new AliITSPIDResponse(fMC);
625 }
626
146d18ee 627 //flags for MC
36be14b3 628 Int_t nTrackMC=0;
629 if(stack) nTrackMC = stack->GetNtrack();
630 const AliESDVertex *vtx = fESD->GetPrimaryVertexSPD();
1f75ac7d 631
a1c5e4a1 632 ///////////selection of the centrality or multiplicity bin
1f75ac7d 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
a1c5e4a1 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
1f75ac7d 657 const AliMultiplicity *mult = fESD->GetMultiplicity();
a1c5e4a1 658 Float_t nClu1 = (Float_t)mult->GetNumberOfITSClusters(1);
659 multiplicity =(Int_t)(AliESDUtils::GetCorrSPD2(nClu1,vtx->GetZ())+0.5);
1f75ac7d 660 }
a1c5e4a1 661 if(fLowMult>-1 && multiplicity<fLowMult) return;
662 if(fUpMult>-1 && multiplicity>fUpMult) return;
663 fHistMult->Fill(multiplicity);
664 fHistNEvents->Fill(2);
36be14b3 665
666 //event selection
36be14b3 667 if(!vtx)evSel=0;
668 else{
a1c5e4a1 669 fHistNEvents->Fill(3);
36be14b3 670 if(vtx->GetNContributors()<0) evSel=0;
671 else{
a1c5e4a1 672 fHistNEvents->Fill(4);
36be14b3 673 if(TMath::Abs(vtx->GetZv())>10) evSel=0;
674 else{
a1c5e4a1 675 fHistNEvents->Fill(5);
36be14b3 676 if(vtx->GetZRes()>0.5) evSel=0;
677 else{
a1c5e4a1 678 fHistNEvents->Fill(6);
36be14b3 679 if(vtx->IsFromVertexerZ() && vtx->GetDispersion()>0.03) evSel=0;
a1c5e4a1 680 else fHistNEvents->Fill(7);
36be14b3 681 }
682 }
683 }
684 }
146d18ee 685
686
7c831837 687 /////first loop on stack, before event selection, filling MC ntuple
688
36be14b3 689 for(Int_t imc=0; imc<nTrackMC; imc++){
690 part = stack->Particle(imc);
10ba520e 691 isph=1;
692 if(!stack->IsPhysicalPrimary(imc)) isph=0;
36be14b3 693 pdgPart = part->GetPDG();
10ba520e 694 if(!pdgPart)continue;
36be14b3 695 if(pdgPart->Charge()==0) continue; //no neutral particles
5c112001 696 Float_t yMC=-999.;
36be14b3 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
36be14b3 699 if(pdgPart->Charge()>0) signMC=1;
700 else signMC=-1;
701 ptMC=part->Pt();
702 code=pdgPart->PdgCode();
10ba520e 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 }
7c831837 716 uniqueID = part->GetUniqueID();
36be14b3 717
36be14b3 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();
10ba520e 730
36be14b3 731 if(fFillNtuple) fNtupleMC->Fill(xntMC);
732 }
733
10ba520e 734 if(jpart>=0){
735 if(stack->IsPhysicalPrimary(imc)){
a1c5e4a1 736 if(signMC>0) fHistPrimMCposBefEvSel[jpart]->Fill(ptMC);
737 else fHistPrimMCnegBefEvSel[jpart]->Fill(ptMC);
10ba520e 738 if(evSel==1){
a1c5e4a1 739 if(signMC>0) fHistPrimMCpos[jpart]->Fill(ptMC);
740 else fHistPrimMCneg[jpart]->Fill(ptMC);
10ba520e 741 }
742 }else{
7c831837 743 if(mfl==3 && uniqueID == kPDecay){ // If a particle is not a physical primary, check if it comes from weak decay
a1c5e4a1 744 if(signMC>0) fHistSecStrMCposBefEvSel[jpart]->Fill(ptMC);
745 else fHistSecStrMCnegBefEvSel[jpart]->Fill(ptMC);
10ba520e 746 if(evSel==1){
a1c5e4a1 747 if(signMC>0) fHistSecStrMCpos[jpart]->Fill(ptMC);
748 else fHistSecStrMCneg[jpart]->Fill(ptMC);
10ba520e 749 }
750 }else{
a1c5e4a1 751 if(signMC>0) fHistSecMatMCposBefEvSel[jpart]->Fill(ptMC);
752 else fHistSecMatMCnegBefEvSel[jpart]->Fill(ptMC);
10ba520e 753 if(evSel==1){
a1c5e4a1 754 if(signMC>0) fHistSecMatMCpos[jpart]->Fill(ptMC);
755 else fHistSecMatMCneg[jpart]->Fill(ptMC);
10ba520e 756 }
36be14b3 757 }
758 }
10ba520e 759 }
760 }
761
36be14b3 762
10ba520e 763 if(evSel==0)return; //event selection
1f75ac7d 764
765
36be14b3 766 //loop on tracks
767 for (Int_t iTrack=0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
10ba520e 768 isph=-999;
36be14b3 769 code=-999;
770 mfl=-999;
7c831837 771 uniqueID=-999;
10ba520e 772
36be14b3 773 track = (AliESDtrack*)fESD->GetTrack(iTrack);
774 if (!track) continue;
775
776 //track selection
c841b89f 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
36be14b3 789 status=track->GetStatus();
790 if((status&AliESDtrack::kITSpureSA)==0) continue; //its standalone
c841b89f 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
36be14b3 801 if((status&AliESDtrack::kITSrefit)==0) continue; //its refit
c841b89f 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
029e1796 812 if(TMath::Abs(track->GetSign())<0.0001) continue; //no neutral particles
c841b89f 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++;
10ba520e 822
36be14b3 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;
c841b89f 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
146d18ee 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
c841b89f 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
36be14b3 851 //chisquare/nclusters
146d18ee 852 Int_t nclu=nSPD+nPtsForPid;
36be14b3 853 if(track->GetITSchi2()/nclu > fMaxChi2Clu) continue;
c841b89f 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
36be14b3 864 //pseudorapidity and rapidity
dac69cb0 865 if(TMath::Abs(track->Eta()) > fEtaRange) continue;
c841b89f 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
36be14b3 876 //truncated mean
146d18ee 877 track->GetITSdEdxSamples(dedxLay);
878 Double_t dedx = CookdEdx(dedxLay);
36be14b3 879 if(dedx<0) continue;
36be14b3 880
c841b89f 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
36be14b3 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);
10ba520e 899
36be14b3 900 //Filling Ntuple
901 //information from the MC kinematics
902 if(fMC){
10ba520e 903 if(track->GetLabel()<0)isph=-1;
36be14b3 904 if(track->GetLabel()>=0){
905 part = (TParticle*)stack->Particle(track->GetLabel());
906 pdgPart = part->GetPDG();
907 code = pdgPart->PdgCode();
10ba520e 908 if(stack->IsPhysicalPrimary(track->GetLabel())) isph=1;
36be14b3 909 else{
10ba520e 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 }
7c831837 917 uniqueID = part->GetUniqueID();
10ba520e 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
7c831837 936 if(mfl==3 && uniqueID == kPDecay){
10ba520e 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 }
36be14b3 958 }
959 }
960 }
7c831837 961 Float_t xnt[19];
36be14b3 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;
146d18ee 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;
36be14b3 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;
c841b89f 979 xnt[index++]=(Float_t)mfl;
7c831837 980 xnt[index++]=(Float_t)uniqueID;
c841b89f 981 xnt[index]=(Float_t)track->GetITSchi2()/nclu;
36be14b3 982
983 if(fFillNtuple) fNtupleNSigma->Fill(xnt);
984
985
986
987 //Compute y and bb
10ba520e 988 Double_t y[4],bbtheo[4],logdiff[4];
36be14b3 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 }
98bb3219 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
10ba520e 1001 for(Int_t i=0;i<4;i++){
36be14b3 1002 y[i] = Eta2y(pt,pdgmass[i],track->Eta());
7fcf3a38 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)
36be14b3 1005 logdiff[i]=TMath::Log(dedx) - TMath::Log(bbtheo[i]);
1006 }
98bb3219 1007
146d18ee 1008 Int_t resocls=(Int_t)nPtsForPid-1;
10ba520e 1009
1010 //NSigma Method, with asymmetric bands
10ba520e 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
98bb3219 1029
36be14b3 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 }
98bb3219 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)));
36be14b3 1043 }
10ba520e 1044
1045 // y calculation
c92ab4a5 1046 Double_t yPartMean=-999;
1047 Double_t yPart=-999;
1048 if(minPosMean>-1) yPartMean=y[minPosMean];
1049 if(minPos>-1) yPart=y[minPos];
36be14b3 1050
34f309eb 1051 if(TMath::Abs(yPartMean)<fMaxY){
10ba520e 1052 //DCA distributions, before the DCA cuts, based on asymmetrinc nsigma approach
36be14b3 1053 if(theBin>=0 && theBin<kNbins){
1054 if(track->GetSign()>0){
10ba520e 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);
36be14b3 1058 }else{
10ba520e 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);
36be14b3 1062 }
1063 }
1064 }
1065
1066 //DCA cut on xy and z
a1c5e4a1 1067 if(!DCAcut(impactXY,impactZ,pt)) continue;
c841b89f 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++;
16c730d4 1077
1078
1079 Int_t jpart=-1;
10ba520e 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();
10ba520e 1089 for(Int_t j=0; j<3; j++){
1090 if(TMath::Abs(code)==listcode[j]){
1091 jpart=j;
1092 break;
1093 }
1094 }
16c730d4 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())){
a1c5e4a1 1100 if(signMC>0) fHistPrimMCposReco[jpart]->Fill(ptMC);
1101 else fHistPrimMCnegReco[jpart]->Fill(ptMC);
16c730d4 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 }
7c831837 1109 uniqueID = part->GetUniqueID();
1110 if(mfl==3 && uniqueID == kPDecay){ // strangeness
a1c5e4a1 1111 if(signMC>0) fHistSecStrMCposReco[jpart]->Fill(ptMC);
1112 else fHistSecStrMCnegReco[jpart]->Fill(ptMC);
16c730d4 1113 }else{
a1c5e4a1 1114 if(signMC>0) fHistSecMatMCposReco[jpart]->Fill(ptMC);
1115 else fHistSecMatMCnegReco[jpart]->Fill(ptMC);
16c730d4 1116 }
10ba520e 1117 }
1118 }
1119 }
1120 }
1121
1122 //Nsigma histos with MC truth
1123
1124 //asymmetric approach
c92ab4a5 1125 if(TMath::Abs(yPartMean)<fMaxY && minPosMean>-1){
10ba520e 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 }
36be14b3 1147
10ba520e 1148 //symmetric bands
34f309eb 1149 if(min<fMinNSigma && TMath::Abs(yPart)<fMaxY){
10ba520e 1150 //nsigma histos
36be14b3 1151 if(track->GetSign()>0) fHistPosNSigma[minPos]->Fill(pt);
1152 else fHistNegNSigma[minPos]->Fill(pt);
1153 if(fMC){
10ba520e 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
36be14b3 1160 if(isph==1){
1161 if(track->GetSign()>0) fHistPosNSigmaPrim[minPos]->Fill(pt);
1162 else fHistNegNSigmaPrim[minPos]->Fill(pt);
16c730d4 1163 //nsigma histos with MC truth on IsPhysicalPrimary and PID
36be14b3 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
10ba520e 1172
1173 //integral approach histograms
36be14b3 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){
e6076bb8 1180 if(TMath::Abs(y[0])<fMaxY){
16c730d4 1181 if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypPion[theBin]->Fill(logdiff[0]);
1182 if(TMath::Abs(code)==11)fHistMCPosElHypPion[theBin]->Fill(logdiff[0]);
e6076bb8 1183 if(TMath::Abs(code)==211)fHistMCPosPiHypPion[theBin]->Fill(logdiff[0]);
98bb3219 1184 if(TMath::Abs(code)==321)fHistMCPosKHypPion[theBin]->Fill(logdiff[0]);
e6076bb8 1185 if(TMath::Abs(code)==2212)fHistMCPosPHypPion[theBin]->Fill(logdiff[0]);
1186 }
1187 if(TMath::Abs(y[1])<fMaxY){
16c730d4 1188 if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypKaon[theBin]->Fill(logdiff[1]);
1189 if(TMath::Abs(code)==11)fHistMCPosElHypKaon[theBin]->Fill(logdiff[1]);
e6076bb8 1190 if(TMath::Abs(code)==211)fHistMCPosPiHypKaon[theBin]->Fill(logdiff[1]);
98bb3219 1191 if(TMath::Abs(code)==321)fHistMCPosKHypKaon[theBin]->Fill(logdiff[1]);
e6076bb8 1192 if(TMath::Abs(code)==2212)fHistMCPosPHypKaon[theBin]->Fill(logdiff[1]);
1193 }
1194 if(TMath::Abs(y[2])<fMaxY){
16c730d4 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]);
e6076bb8 1200 }
36be14b3 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){
e6076bb8 1207 if(TMath::Abs(y[0])<fMaxY){
16c730d4 1208 if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypPion[theBin]->Fill(logdiff[0]);
1209 if(TMath::Abs(code)==11)fHistMCNegElHypPion[theBin]->Fill(logdiff[0]);
98bb3219 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]);
e6076bb8 1213 }
1214 if(TMath::Abs(y[1])<fMaxY){
16c730d4 1215 if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypKaon[theBin]->Fill(logdiff[1]);
1216 if(TMath::Abs(code)==11)fHistMCNegElHypKaon[theBin]->Fill(logdiff[1]);
98bb3219 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]);
e6076bb8 1220 }
1221 if(TMath::Abs(y[2])<fMaxY){
16c730d4 1222 if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypProton[theBin]->Fill(logdiff[2]);
1223 if(TMath::Abs(code)==11)fHistMCNegElHypProton[theBin]->Fill(logdiff[2]);
98bb3219 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]);
e6076bb8 1227 }
36be14b3 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++){
146d18ee 1239 if(dedxLay[j]<5) continue;
1240 fHistCharge[j]->Fill(dedxLay[j]);
36be14b3 1241 }
1242 }
1243
1244 // Post output data.
1245 PostData(1,fOutput);
a1c5e4a1 1246 PostData(2,fListCuts);
1247 // Printf("............. end of Exec");
36be14b3 1248}
1249
1250//________________________________________________________________________
1251void AliAnalysisTaskSEITSsaSpectra::Terminate(Option_t *) {
4d668cf7 1252 // Merge output
1253 // Called once at the end of the query
36be14b3 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"));
a1c5e4a1 1261 printf("Number of Analyzed Events = %f\n",fHistNEvents->GetBinContent(1));
36be14b3 1262
1263 Printf("end of Terminate");
1264 return;
1265}