]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/SPECTRA/PiKaPr/ITSsa/AliAnalysisTaskSEITSsaSpectra.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[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);
dc3531c1 343 fHistPrimMCposBefEvSelEta[j] = new TH1F(Form("fHistPrimMCposBefEvSelEta%d",j),Form("fHistPrimMCposBefEvSelEta%d",j),kNbins,fPtBinLimits);
344 fHistPrimMCposBefEvSelEtaY[j] = new TH1F(Form("fHistPrimMCposBefEvSelEtaY%d",j),Form("fHistPrimMCposBefEvSelEtaY%d",j),kNbins,fPtBinLimits);
10ba520e 345 fHistPrimMCnegBefEvSel[j] = new TH1F(Form("fHistPrimMCnegBefEvSel%d",j),Form("fHistPrimMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
dc3531c1 346 fHistPrimMCnegBefEvSelEta[j] = new TH1F(Form("fHistPrimMCnegBefEvSelEta%d",j),Form("fHistPrimMCnegBefEvSelEta%d",j),kNbins,fPtBinLimits);
347 fHistPrimMCnegBefEvSelEtaY[j] = new TH1F(Form("fHistPrimMCnegBefEvSelEtaY%d",j),Form("fHistPrimMCnegBefEvSelEtaY%d",j),kNbins,fPtBinLimits);
10ba520e 348 fOutput->Add(fHistPrimMCposBefEvSel[j]);
dc3531c1 349 fOutput->Add(fHistPrimMCposBefEvSelEta[j]);
350 fOutput->Add(fHistPrimMCposBefEvSelEtaY[j]);
351 fOutput->Add(fHistPrimMCnegBefEvSel[j]);
352 fOutput->Add(fHistPrimMCnegBefEvSelEta[j]);
353 fOutput->Add(fHistPrimMCnegBefEvSelEtaY[j]);
10ba520e 354 fHistSecStrMCposBefEvSel[j] = new TH1F(Form("fHistSecStrMCposBefEvSel%d",j),Form("fHistSecStrMCposBefEvSel%d",j),kNbins,fPtBinLimits);
355 fHistSecStrMCnegBefEvSel[j] = new TH1F(Form("fHistSecStrMCnegBefEvSel%d",j),Form("fHistSecStrMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
356 fOutput->Add(fHistSecStrMCnegBefEvSel[j]);
357 fOutput->Add(fHistSecStrMCposBefEvSel[j]);
358 fHistSecMatMCposBefEvSel[j] = new TH1F(Form("fHistSecMatMCposBefEvSel%d",j),Form("fHistSecMatMCposBefEvSel%d",j),kNbins,fPtBinLimits);
359 fHistSecMatMCnegBefEvSel[j] = new TH1F(Form("fHistSecMatMCnegBefEvSel%d",j),Form("fHistSecMatMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
360 fOutput->Add(fHistSecMatMCnegBefEvSel[j]);
361 fOutput->Add(fHistSecMatMCposBefEvSel[j]);
362 //
363 fHistPrimMCposReco[j] = new TH1F(Form("fHistPrimMCposReco%d",j),Form("fHistPrimMCposReco%d",j),kNbins,fPtBinLimits);
dc3531c1 364 fHistPrimMCposRecoEtaY[j] = new TH1F(Form("fHistPrimMCposRecoEtaY%d",j),Form("fHistPrimMCposRecoEtaY%d",j),kNbins,fPtBinLimits);
10ba520e 365 fHistPrimMCnegReco[j] = new TH1F(Form("fHistPrimMCnegReco%d",j),Form("fHistPrimMCnegReco%d",j),kNbins,fPtBinLimits);
dc3531c1 366 fHistPrimMCnegRecoEtaY[j] = new TH1F(Form("fHistPrimMCnegRecoEtaY%d",j),Form("fHistPrimMCnegRecoEtaY%d",j),kNbins,fPtBinLimits);
10ba520e 367 fOutput->Add(fHistPrimMCposReco[j]);
dc3531c1 368 fOutput->Add(fHistPrimMCposRecoEtaY[j]);
369 fOutput->Add(fHistPrimMCnegReco[j]);
370 fOutput->Add(fHistPrimMCnegRecoEtaY[j]);
10ba520e 371 fHistSecStrMCposReco[j] = new TH1F(Form("fHistSecStrMCposReco%d",j),Form("fHistSecStrMCposReco%d",j),kNbins,fPtBinLimits);
372 fHistSecStrMCnegReco[j] = new TH1F(Form("fHistSecStrMCnegReco%d",j),Form("fHistSecStrMCnegReco%d",j),kNbins,fPtBinLimits);
373 fOutput->Add(fHistSecStrMCnegReco[j]);
374 fOutput->Add(fHistSecStrMCposReco[j]);
375 fHistSecMatMCposReco[j] = new TH1F(Form("fHistSecMatMCposReco%d",j),Form("fHistSecMatMCposReco%d",j),kNbins,fPtBinLimits);
376 fHistSecMatMCnegReco[j] = new TH1F(Form("fHistSecMatMCnegReco%d",j),Form("fHistSecMatMCnegReco%d",j),kNbins,fPtBinLimits);
377 fOutput->Add(fHistSecMatMCnegReco[j]);
378 fOutput->Add(fHistSecMatMCposReco[j]);
36be14b3 379
36be14b3 380 }
381
382 for(Int_t i=0; i<4; i++){
383 fHistCharge[i] = new TH1F(Form("fHistChargeLay%d",i),Form("fHistChargeLay%d",i),100,0,300);
384 fOutput->Add(fHistCharge[i]);
385 }
386
387 for(Int_t i=0; i<kNbins; i++){
388 fHistPosPi[i] = new TH1F(Form("fHistPosPi%d",i),Form("fHistPosPi%d",i),175,-3.5,3.5);
389 fHistPosK[i] = new TH1F(Form("fHistPosK%d",i),Form("fHistPosK%d",i),175,-3.5,3.5);
390 fHistPosP[i] = new TH1F(Form("fHistPosP%d",i),Form("fHistPosP%d",i),175,-3.5,3.5);
391 fHistNegPi[i] = new TH1F(Form("fHistNegPi%d",i),Form("fHistNegPi%d",i),175,-3.5,3.5);
392 fHistNegK[i] = new TH1F(Form("fHistNegK%d",i),Form("fHistNegK%d",i),175,-3.5,3.5);
393 fHistNegP[i] = new TH1F(Form("fHistNegP%d",i),Form("fHistNegP%d",i),175,-3.5,3.5);
394
10ba520e 395 fHistDCAPosPi[i] = new TH1F(Form("fHistDCAPosPi%d",i),Form("fHistDCAPosPi%d",i),2000,-1,1); //DCA distr. with NSigma PID
36be14b3 396 fHistDCAPosK[i] = new TH1F(Form("fHistDCAPosK%d",i),Form("fHistDCAPosK%d",i),2000,-1,1);
397 fHistDCAPosP[i] = new TH1F(Form("fHistDCAPosP%d",i),Form("fHistDCAPosP%d",i),2000,-1,1);
398 fHistDCANegPi[i] = new TH1F(Form("fHistDCANegPi%d",i),Form("fHistDCANegPi%d",i),2000,-1,1);
399 fHistDCANegK[i] = new TH1F(Form("fHistDCANegK%d",i),Form("fHistDCANegK%d",i),2000,-1,1);
400 fHistDCANegP[i] = new TH1F(Form("fHistDCANegP%d",i),Form("fHistDCANegP%d",i),2000,-1,1);
401
10ba520e 402 fHistMCPrimDCAPosPi[i] = new TH1F(Form("fHistMCPrimDCAPosPi%d",i),Form("fHistMCPrimDCAPosPi%d",i),2000,-1,1); //DCA distr. with MC truth
403 fHistMCPrimDCAPosK[i] = new TH1F(Form("fHistMCPrimDCAPosK%d",i),Form("fHistMCPrimDCAPosK%d",i),2000,-1,1);
404 fHistMCPrimDCAPosP[i] = new TH1F(Form("fHistMCPrimDCAPosP%d",i),Form("fHistMCPrimDCAPosP%d",i),2000,-1,1);
405 fHistMCPrimDCANegPi[i] = new TH1F(Form("fHistMCPrimDCANegPi%d",i),Form("fHistMCPrimDCANegPi%d",i),2000,-1,1);
406 fHistMCPrimDCANegK[i] = new TH1F(Form("fHistMCPrimDCANegK%d",i),Form("fHistMCPrimDCANegK%d",i),2000,-1,1);
407 fHistMCPrimDCANegP[i] = new TH1F(Form("fHistMCPrimDCANegP%d",i),Form("fHistMCPrimDCANegP%d",i),2000,-1,1);
408
409 fHistMCSecStDCAPosPi[i] = new TH1F(Form("fHistMCSecStDCAPosPi%d",i),Form("fHistMCSecStDCAPosPi%d",i),2000,-1,1); //DCA distr. with MC truth
410 fHistMCSecStDCAPosK[i] = new TH1F(Form("fHistMCSecStDCAPosK%d",i),Form("fHistMCSecStDCAPosK%d",i),2000,-1,1);
411 fHistMCSecStDCAPosP[i] = new TH1F(Form("fHistMCSecStDCAPosP%d",i),Form("fHistMCSecStDCAPosP%d",i),2000,-1,1);
412 fHistMCSecStDCANegPi[i] = new TH1F(Form("fHistMCSecStDCANegPi%d",i),Form("fHistMCSecStDCANegPi%d",i),2000,-1,1);
413 fHistMCSecStDCANegK[i] = new TH1F(Form("fHistMCSecStDCANegK%d",i),Form("fHistMCSecStDCANegK%d",i),2000,-1,1);
414 fHistMCSecStDCANegP[i] = new TH1F(Form("fHistMCSecStDCANegP%d",i),Form("fHistMCSecStDCANegP%d",i),2000,-1,1);
415
416 fHistMCSecMatDCAPosPi[i] = new TH1F(Form("fHistMCSecMatDCAPosPi%d",i),Form("fHistMCSecMatDCAPosPi%d",i),2000,-1,1); //DCA distr. with MC truth
417 fHistMCSecMatDCAPosK[i] = new TH1F(Form("fHistMCSecMatDCAPosK%d",i),Form("fHistMCSecMatDCAPosK%d",i),2000,-1,1);
418 fHistMCSecMatDCAPosP[i] = new TH1F(Form("fHistMCSecMatDCAPosP%d",i),Form("fHistMCSecMatDCAPosP%d",i),2000,-1,1);
419 fHistMCSecMatDCANegPi[i] = new TH1F(Form("fHistMCSecMatDCANegPi%d",i),Form("fHistMCSecMatDCANegPi%d",i),2000,-1,1);
420 fHistMCSecMatDCANegK[i] = new TH1F(Form("fHistMCSecMatDCANegK%d",i),Form("fHistMCSecMatDCANegK%d",i),2000,-1,1);
421 fHistMCSecMatDCANegP[i] = new TH1F(Form("fHistMCSecMatDCANegP%d",i),Form("fHistMCSecMatDCANegP%d",i),2000,-1,1);
422
16c730d4 423 fHistMCPosOtherHypPion[i] = new TH1F(Form("fHistMCPosOtherHypPion%d",i),Form("fHistMCPosOtherHypPion%d",i),175,-3.5,3.5); //MC truth
424 fHistMCPosOtherHypKaon[i] = new TH1F(Form("fHistMCPosOtherHypKaon%d",i),Form("fHistMCPosOtherHypKaon%d",i),175,-3.5,3.5);
425 fHistMCPosOtherHypProton[i] = new TH1F(Form("fHistMCPosOtherHypProton%d",i),Form("fHistMCPosOtherHypProton%d",i),175,-3.5,3.5);
426 fHistMCPosElHypPion[i] = new TH1F(Form("fHistMCPosElHypPion%d",i),Form("fHistMCPosElHypPion%d",i),175,-3.5,3.5);
427 fHistMCPosElHypKaon[i] = new TH1F(Form("fHistMCPosElHypKaon%d",i),Form("fHistMCPosElHypKaon%d",i),175,-3.5,3.5);
428 fHistMCPosElHypProton[i] = new TH1F(Form("fHistMCPosElHypProton%d",i),Form("fHistMCPosElHypProton%d",i),175,-3.5,3.5);
429 fHistMCPosPiHypPion[i] = new TH1F(Form("fHistMCPosPiHypPion%d",i),Form("fHistMCPosPiHypPion%d",i),175,-3.5,3.5);
e6076bb8 430 fHistMCPosPiHypKaon[i] = new TH1F(Form("fHistMCPosPiHypKaon%d",i),Form("fHistMCPosPiHypKaon%d",i),175,-3.5,3.5);
431 fHistMCPosPiHypProton[i] = new TH1F(Form("fHistMCPosPiHypProton%d",i),Form("fHistMCPosPiHypProton%d",i),175,-3.5,3.5);
432 fHistMCPosKHypPion[i] = new TH1F(Form("fHistMCPosKHypPion%d",i),Form("fHistMCPosKHypPion%d",i),175,-3.5,3.5);
433 fHistMCPosKHypKaon[i] = new TH1F(Form("fHistMCPosKHypKaon%d",i),Form("fHistMCPosKHypKaon%d",i),175,-3.5,3.5);
434 fHistMCPosKHypProton[i] = new TH1F(Form("fHistMCPosKHypProton%d",i),Form("fHistMCPosKHypProton%d",i),175,-3.5,3.5);
435 fHistMCPosPHypPion[i] = new TH1F(Form("fHistMCPosPHypPion%d",i),Form("fHistMCPosPHypPion%d",i),175,-3.5,3.5);
436 fHistMCPosPHypKaon[i] = new TH1F(Form("fHistMCPosPHypKaon%d",i),Form("fHistMCPosPHypKaon%d",i),175,-3.5,3.5);
437 fHistMCPosPHypProton[i] = new TH1F(Form("fHistMCPosPHypProton%d",i),Form("fHistMCPosPHypProton%d",i),175,-3.5,3.5);
438
16c730d4 439 fHistMCNegOtherHypPion[i] = new TH1F(Form("fHistMCNegOtherHypPion%d",i),Form("fHistMCNegOtherHypPion%d",i),175,-3.5,3.5); //MC truth
440 fHistMCNegOtherHypKaon[i] = new TH1F(Form("fHistMCNegOtherHypKaon%d",i),Form("fHistMCNegOtherHypKaon%d",i),175,-3.5,3.5);
441 fHistMCNegOtherHypProton[i] = new TH1F(Form("fHistMCNegOtherHypProton%d",i),Form("fHistMCNegOtherHypProton%d",i),175,-3.5,3.5);
442 fHistMCNegElHypPion[i] = new TH1F(Form("fHistMCNegElHypPion%d",i),Form("fHistMCNegElHypPion%d",i),175,-3.5,3.5);
443 fHistMCNegElHypKaon[i] = new TH1F(Form("fHistMCNegElHypKaon%d",i),Form("fHistMCNegElHypKaon%d",i),175,-3.5,3.5);
444 fHistMCNegElHypProton[i] = new TH1F(Form("fHistMCNegElHypProton%d",i),Form("fHistMCNegElHypProton%d",i),175,-3.5,3.5);
445 fHistMCNegPiHypPion[i] = new TH1F(Form("fHistMCNegPiHypPion%d",i),Form("fHistMCNegPiHypPion%d",i),175,-3.5,3.5);
e6076bb8 446 fHistMCNegPiHypKaon[i] = new TH1F(Form("fHistMCNegPiHypKaon%d",i),Form("fHistMCNegPiHypKaon%d",i),175,-3.5,3.5);
447 fHistMCNegPiHypProton[i] = new TH1F(Form("fHistMCNegPiHypProton%d",i),Form("fHistMCNegPiHypProton%d",i),175,-3.5,3.5);
448 fHistMCNegKHypPion[i] = new TH1F(Form("fHistMCNegKHypPion%d",i),Form("fHistMCNegKHypPion%d",i),175,-3.5,3.5);
449 fHistMCNegKHypKaon[i] = new TH1F(Form("fHistMCNegKHypKaon%d",i),Form("fHistMCNegKHypKaon%d",i),175,-3.5,3.5);
450 fHistMCNegKHypProton[i] = new TH1F(Form("fHistMCNegKHypProton%d",i),Form("fHistMCNegKHypProton%d",i),175,-3.5,3.5);
451 fHistMCNegPHypPion[i] = new TH1F(Form("fHistMCNegPHypPion%d",i),Form("fHistMCNegPHypPion%d",i),175,-3.5,3.5);
452 fHistMCNegPHypKaon[i] = new TH1F(Form("fHistMCNegPHypKaon%d",i),Form("fHistMCNegPHypKaon%d",i),175,-3.5,3.5);
453 fHistMCNegPHypProton[i] = new TH1F(Form("fHistMCNegPHypProton%d",i),Form("fHistMCNegPHypProton%d",i),175,-3.5,3.5);
454
455
36be14b3 456 fOutput->Add(fHistPosPi[i]);
457 fOutput->Add(fHistPosK[i]);
458 fOutput->Add(fHistPosP[i]);
459 fOutput->Add(fHistNegPi[i]);
460 fOutput->Add(fHistNegK[i]);
461 fOutput->Add(fHistNegP[i]);
462
6b77d2c0 463 fOutput->Add(fHistDCAPosPi[i]); //DCA distr
36be14b3 464 fOutput->Add(fHistDCAPosK[i]);
465 fOutput->Add(fHistDCAPosP[i]);
466 fOutput->Add(fHistDCANegPi[i]);
467 fOutput->Add(fHistDCANegK[i]);
468 fOutput->Add(fHistDCANegP[i]);
469
10ba520e 470 fOutput->Add(fHistMCPrimDCAPosPi[i]);//DCA distr.
471 fOutput->Add(fHistMCPrimDCAPosK[i]);
472 fOutput->Add(fHistMCPrimDCAPosP[i]);
473 fOutput->Add(fHistMCPrimDCANegPi[i]);
474 fOutput->Add(fHistMCPrimDCANegK[i]);
475 fOutput->Add(fHistMCPrimDCANegP[i]);
476
477 fOutput->Add(fHistMCSecStDCAPosPi[i]);//DCA distr.
478 fOutput->Add(fHistMCSecStDCAPosK[i]);
479 fOutput->Add(fHistMCSecStDCAPosP[i]);
480 fOutput->Add(fHistMCSecStDCANegPi[i]);
481 fOutput->Add(fHistMCSecStDCANegK[i]);
482 fOutput->Add(fHistMCSecStDCANegP[i]);
483
484 fOutput->Add(fHistMCSecMatDCAPosPi[i]);//DCA distr.
485 fOutput->Add(fHistMCSecMatDCAPosK[i]);
486 fOutput->Add(fHistMCSecMatDCAPosP[i]);
487 fOutput->Add(fHistMCSecMatDCANegPi[i]);
488 fOutput->Add(fHistMCSecMatDCANegK[i]);
489 fOutput->Add(fHistMCSecMatDCANegP[i]);
490
16c730d4 491 fOutput->Add(fHistMCPosOtherHypPion[i]);//MC truth
492 fOutput->Add(fHistMCPosOtherHypKaon[i]);
493 fOutput->Add(fHistMCPosOtherHypProton[i]);
494 fOutput->Add(fHistMCPosElHypPion[i]);
495 fOutput->Add(fHistMCPosElHypKaon[i]);
496 fOutput->Add(fHistMCPosElHypProton[i]);
497 fOutput->Add(fHistMCPosPiHypPion[i]);
e6076bb8 498 fOutput->Add(fHistMCPosPiHypKaon[i]);
499 fOutput->Add(fHistMCPosPiHypProton[i]);
500 fOutput->Add(fHistMCPosKHypPion[i]);
501 fOutput->Add(fHistMCPosKHypKaon[i]);
502 fOutput->Add(fHistMCPosKHypProton[i]);
503 fOutput->Add(fHistMCPosPHypPion[i]);
504 fOutput->Add(fHistMCPosPHypKaon[i]);
505 fOutput->Add(fHistMCPosPHypProton[i]);
506
16c730d4 507 fOutput->Add(fHistMCNegOtherHypPion[i]);//MC truth
508 fOutput->Add(fHistMCNegOtherHypKaon[i]);
509 fOutput->Add(fHistMCNegOtherHypProton[i]);
510 fOutput->Add(fHistMCNegElHypPion[i]);
511 fOutput->Add(fHistMCNegElHypKaon[i]);
512 fOutput->Add(fHistMCNegElHypProton[i]);
513 fOutput->Add(fHistMCNegPiHypPion[i]);
e6076bb8 514 fOutput->Add(fHistMCNegPiHypKaon[i]);
515 fOutput->Add(fHistMCNegPiHypProton[i]);
516 fOutput->Add(fHistMCNegKHypPion[i]);
517 fOutput->Add(fHistMCNegKHypKaon[i]);
518 fOutput->Add(fHistMCNegKHypProton[i]);
519 fOutput->Add(fHistMCNegPHypPion[i]);
520 fOutput->Add(fHistMCNegPHypKaon[i]);
521 fOutput->Add(fHistMCNegPHypProton[i]);
522
36be14b3 523 }
524
525 //NSigma Histos
526 for(Int_t j=0;j<3;j++){
10ba520e 527
528 fHistPosNSigmaMean[j] = new TH1F(Form("hHistPosNSigmaMean%d",j),Form("hHistPosNSigmaMean%d",j),kNbins,fPtBinLimits);
529 fHistNegNSigmaMean[j] = new TH1F(Form("hHistNegNSigmaMean%d",j),Form("hHistNegNSigmaMean%d",j),kNbins,fPtBinLimits);
530 fHistPosNSigmaMCMean[j] = new TH1F(Form("hHistPosNSigmaMCMean%d",j),Form("hHistPosNSigmaMCMean%d",j),kNbins,fPtBinLimits);
531 fHistNegNSigmaMCMean[j] = new TH1F(Form("hHistNegNSigmaMCMean%d",j),Form("hHistNegNSigmaMCMean%d",j),kNbins,fPtBinLimits);
532 fHistPosNSigmaPrimMean[j] = new TH1F(Form("hHistPosNSigmaPrimMean%d",j),Form("hHistPosNSigmaPrimMean%d",j),kNbins,fPtBinLimits);
533 fHistNegNSigmaPrimMean[j] = new TH1F(Form("hHistNegNSigmaPrimMean%d",j),Form("hHistNegNSigmaPrimMean%d",j),kNbins,fPtBinLimits);
534 fHistPosNSigmaPrimMCMean[j] = new TH1F(Form("hHistPosNSigmaPrimMCMean%d",j),Form("hHistPosNSigmaPrimMCMean%d",j),kNbins,fPtBinLimits);
535 fHistNegNSigmaPrimMCMean[j] = new TH1F(Form("hHistNegNSigmaPrimMCMean%d",j),Form("hHistNegNSigmaPrimMCMean%d",j),kNbins,fPtBinLimits);
536 fOutput->Add(fHistPosNSigmaMean[j]);
537 fOutput->Add(fHistNegNSigmaMean[j]);
538 fOutput->Add(fHistPosNSigmaMCMean[j]);
539 fOutput->Add(fHistNegNSigmaMCMean[j]);
540 fOutput->Add(fHistPosNSigmaPrimMean[j]);
541 fOutput->Add(fHistNegNSigmaPrimMean[j]);
542 fOutput->Add(fHistPosNSigmaPrimMCMean[j]);
543 fOutput->Add(fHistNegNSigmaPrimMCMean[j]);
544
36be14b3 545 fHistPosNSigma[j] = new TH1F(Form("hHistPosNSigma%d",j),Form("hHistPosNSigma%d",j),kNbins,fPtBinLimits);
546 fHistNegNSigma[j] = new TH1F(Form("hHistNegNSigma%d",j),Form("hHistNegNSigma%d",j),kNbins,fPtBinLimits);
10ba520e 547 fHistPosNSigmaMC[j] = new TH1F(Form("hHistPosNSigmaMC%d",j),Form("hHistPosNSigmaMC%d",j),kNbins,fPtBinLimits);
548 fHistNegNSigmaMC[j] = new TH1F(Form("hHistNegNSigmaMC%d",j),Form("hHistNegNSigmaMC%d",j),kNbins,fPtBinLimits);
36be14b3 549 fHistPosNSigmaPrim[j] = new TH1F(Form("hHistPosNSigmaPrim%d",j),Form("hHistPosNSigmaPrim%d",j),kNbins,fPtBinLimits);
550 fHistNegNSigmaPrim[j] = new TH1F(Form("hHistNegNSigmaPrim%d",j),Form("hHistNegNSigmaPrim%d",j),kNbins,fPtBinLimits);
551 fHistPosNSigmaPrimMC[j] = new TH1F(Form("hHistPosNSigmaPrimMC%d",j),Form("hHistPosNSigmaPrimMC%d",j),kNbins,fPtBinLimits);
552 fHistNegNSigmaPrimMC[j] = new TH1F(Form("hHistNegNSigmaPrimMC%d",j),Form("hHistNegNSigmaPrimMC%d",j),kNbins,fPtBinLimits);
553 fOutput->Add(fHistPosNSigma[j]);
554 fOutput->Add(fHistNegNSigma[j]);
10ba520e 555 fOutput->Add(fHistPosNSigmaMC[j]);
556 fOutput->Add(fHistNegNSigmaMC[j]);
36be14b3 557 fOutput->Add(fHistPosNSigmaPrim[j]);
558 fOutput->Add(fHistNegNSigmaPrim[j]);
559 fOutput->Add(fHistPosNSigmaPrimMC[j]);
560 fOutput->Add(fHistNegNSigmaPrimMC[j]);
10ba520e 561
36be14b3 562 }
563
7c831837 564 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 565 fOutput->Add(fNtupleNSigma);
566 fNtupleMC = new TNtuple("fNtupleMC","fNtupleMC","ptMC:pdgcode:signMC:etaMC:yMC:isph:evSel:run");
567 fOutput->Add(fNtupleMC);
7f8f788e 568
a1c5e4a1 569
570
7f8f788e 571 // Post output data.
572 PostData(1,fOutput);
573
36be14b3 574 Printf("end of CreateOutputObjects");
575}
576
577//________________________________________________________________________
578void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
4d668cf7 579 // Main loop
580 // Called for each event
36be14b3 581
1f75ac7d 582 ///////////////////////////////////////
36be14b3 583 //variables
10ba520e 584 Float_t pdgmass[4]={0.13957,0.493677,0.938272,1.8756}; //mass for pi, K, P (Gev/c^2)
36be14b3 585 Int_t listcode[3]={211,321,2212};//code for pi, K, P (Gev/c^2)
146d18ee 586 Double_t dedxLay[4];
5c112001 587 Float_t ptMC=-999;
7c831837 588 Int_t code=-999, signMC=-999,isph=-999,mfl=-999,uniqueID=-999;
36be14b3 589 Float_t impactXY=-999, impactZ=-999;
590 Int_t evSel=1;
591 AliESDtrack* track;
592 UInt_t status;
593 AliStack* stack=0;
594 TParticle *part=0;
595 TParticlePDG *pdgPart=0;
596
36be14b3 597 /////////////////////
146d18ee 598
599 fESD=(AliESDEvent*)InputEvent();
600 if(!fESD) {
601 printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
602 return;
603 }
604 fHistNEvents->Fill(-1);
605
a1c5e4a1 606 UInt_t maskPhysSel = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
607 TString firedTriggerClasses=fESD->GetFiredTriggerClasses();
608 // if(!firedTriggerClasses.Contains("CINT1B")) return;
609 if((maskPhysSel & AliVEvent::kMB)==0) return;
146d18ee 610 fHistNEvents->Fill(0);
611
a1c5e4a1 612 if(fLowEnergypp && !fMC){ // remove events without SDD in pp 2.76 TeV
613 if(!firedTriggerClasses.Contains("ALL")) return;
614 }
615 fHistNEvents->Fill(1);
616
146d18ee 617
36be14b3 618 if(fMC){
619 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
620 if (!eventHandler) {
621 Printf("ERROR: Could not retrieve MC event handler");
622 return;
623 }
624 AliMCEvent* mcEvent = eventHandler->MCEvent();
625 if (!mcEvent) {
626 Printf("ERROR: Could not retrieve MC event");
627 return;
628 }
629 stack = mcEvent->Stack();
630 if (!stack) {
631 printf("ERROR: stack not available\n");
632 return;
633 }
634 }
98bb3219 635 if(!fITSPIDResponse){
636 fITSPIDResponse=new AliITSPIDResponse(fMC);
637 }
638
146d18ee 639 //flags for MC
36be14b3 640 Int_t nTrackMC=0;
641 if(stack) nTrackMC = stack->GetNtrack();
642 const AliESDVertex *vtx = fESD->GetPrimaryVertexSPD();
1f75ac7d 643
a1c5e4a1 644 ///////////selection of the centrality or multiplicity bin
1f75ac7d 645
646 //selection on the event centrality
647 if(fHImode){
648 if(!(fLowCentrality<0.0)&&fUpCentrality>0.0)
649 {
650 AliCentrality *centrality = fESD->GetCentrality();
651 if(!centrality->IsEventInCentralityClass(fLowCentrality,fUpCentrality,"V0M"))
652 return;
653 Printf("Centrality of the event: %.1f",centrality->GetCentralityPercentile("V0M"));
654 Printf("Centrality cut: %.1f to %.1f",fLowCentrality,fUpCentrality);
655 fHistCen->Fill(centrality->GetCentralityPercentile("V0M"));
656 }
657 }
658
659 //selection on the event multiplicity based on global tracks
a1c5e4a1 660 Int_t multiplicity = -1;
661 if(fMultEstimator==0){
662 // tracks+tracklets
663 multiplicity=AliESDtrackCuts::GetReferenceMultiplicity(fESD, AliESDtrackCuts::kTrackletsITSTPC, 0.8);
664 }else if(fMultEstimator==1){
665 // tracklets
666 multiplicity=AliESDtrackCuts::GetReferenceMultiplicity(fESD, AliESDtrackCuts::kTracklets, 0.8);
667 }else if(fMultEstimator==2){
668 // clusters in SPD1
1f75ac7d 669 const AliMultiplicity *mult = fESD->GetMultiplicity();
a1c5e4a1 670 Float_t nClu1 = (Float_t)mult->GetNumberOfITSClusters(1);
671 multiplicity =(Int_t)(AliESDUtils::GetCorrSPD2(nClu1,vtx->GetZ())+0.5);
1f75ac7d 672 }
a1c5e4a1 673 if(fLowMult>-1 && multiplicity<fLowMult) return;
674 if(fUpMult>-1 && multiplicity>fUpMult) return;
675 fHistMult->Fill(multiplicity);
676 fHistNEvents->Fill(2);
36be14b3 677
678 //event selection
36be14b3 679 if(!vtx)evSel=0;
680 else{
a1c5e4a1 681 fHistNEvents->Fill(3);
36be14b3 682 if(vtx->GetNContributors()<0) evSel=0;
683 else{
a1c5e4a1 684 fHistNEvents->Fill(4);
e690d4d0 685 if(TMath::Abs(vtx->GetZ())>10) evSel=0;
36be14b3 686 else{
a1c5e4a1 687 fHistNEvents->Fill(5);
36be14b3 688 if(vtx->GetZRes()>0.5) evSel=0;
689 else{
a1c5e4a1 690 fHistNEvents->Fill(6);
36be14b3 691 if(vtx->IsFromVertexerZ() && vtx->GetDispersion()>0.03) evSel=0;
a1c5e4a1 692 else fHistNEvents->Fill(7);
36be14b3 693 }
694 }
695 }
696 }
146d18ee 697
698
7c831837 699 /////first loop on stack, before event selection, filling MC ntuple
700
36be14b3 701 for(Int_t imc=0; imc<nTrackMC; imc++){
702 part = stack->Particle(imc);
10ba520e 703 isph=1;
704 if(!stack->IsPhysicalPrimary(imc)) isph=0;
36be14b3 705 pdgPart = part->GetPDG();
10ba520e 706 if(!pdgPart)continue;
36be14b3 707 if(pdgPart->Charge()==0) continue; //no neutral particles
5c112001 708 Float_t yMC=-999.;
36be14b3 709 if(part->Energy() != TMath::Abs(part->Pz())) yMC = 0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
36be14b3 710 if(pdgPart->Charge()>0) signMC=1;
711 else signMC=-1;
712 ptMC=part->Pt();
713 code=pdgPart->PdgCode();
10ba520e 714 Int_t jpart=-1;
715 for(Int_t j=0; j<3; j++){
716 if(TMath::Abs(code)==listcode[j]){
717 jpart=j;
718 break;
719 }
720 }
dc3531c1 721
722 if(jpart>=0 && stack->IsPhysicalPrimary(imc) && TMath::Abs(part->Eta())< fEtaRange){
723 if(signMC>0) fHistPrimMCposBefEvSelEta[jpart]->Fill(ptMC);
724 else fHistPrimMCnegBefEvSelEta[jpart]->Fill(ptMC);
725 }
726
727 if(TMath::Abs(yMC) > fMaxY) continue; //rapidity cut
728
10ba520e 729 Int_t indexMoth=part->GetFirstMother();
730 if(indexMoth>=0){
731 TParticle* moth = stack->Particle(indexMoth);
732 Float_t codemoth = TMath::Abs(moth->GetPdgCode());
733 mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
734 }
7c831837 735 uniqueID = part->GetUniqueID();
36be14b3 736
36be14b3 737 //filling MC ntuple
738 if(TMath::Abs(code)==211 || TMath::Abs(code)==321 || TMath::Abs(code)==2212){
739 Float_t xntMC[8];
740 Int_t indexMC=0;
741 xntMC[indexMC++]=(Float_t)ptMC;
742 xntMC[indexMC++]=(Float_t)code;
743 xntMC[indexMC++]=(Float_t)signMC;
744 xntMC[indexMC++]=(Float_t)part->Eta();
745 xntMC[indexMC++]=(Float_t)yMC;
746 xntMC[indexMC++]=(Float_t)isph;
747 xntMC[indexMC++]=(Float_t)evSel;
748 xntMC[indexMC++]=(Float_t)fESD->GetRunNumber();
10ba520e 749
36be14b3 750 if(fFillNtuple) fNtupleMC->Fill(xntMC);
751 }
752
10ba520e 753 if(jpart>=0){
754 if(stack->IsPhysicalPrimary(imc)){
dc3531c1 755 if(signMC>0){
756 fHistPrimMCposBefEvSel[jpart]->Fill(ptMC);
757 if(TMath::Abs(part->Eta())<fEtaRange) fHistPrimMCposBefEvSelEtaY[jpart]->Fill(ptMC);
758 }
759 else{
760 fHistPrimMCnegBefEvSel[jpart]->Fill(ptMC);
761 if(TMath::Abs(part->Eta())<fEtaRange) fHistPrimMCnegBefEvSelEtaY[jpart]->Fill(ptMC);
762 }
10ba520e 763 if(evSel==1){
a1c5e4a1 764 if(signMC>0) fHistPrimMCpos[jpart]->Fill(ptMC);
765 else fHistPrimMCneg[jpart]->Fill(ptMC);
10ba520e 766 }
767 }else{
7c831837 768 if(mfl==3 && uniqueID == kPDecay){ // If a particle is not a physical primary, check if it comes from weak decay
a1c5e4a1 769 if(signMC>0) fHistSecStrMCposBefEvSel[jpart]->Fill(ptMC);
770 else fHistSecStrMCnegBefEvSel[jpart]->Fill(ptMC);
10ba520e 771 if(evSel==1){
a1c5e4a1 772 if(signMC>0) fHistSecStrMCpos[jpart]->Fill(ptMC);
773 else fHistSecStrMCneg[jpart]->Fill(ptMC);
10ba520e 774 }
775 }else{
a1c5e4a1 776 if(signMC>0) fHistSecMatMCposBefEvSel[jpart]->Fill(ptMC);
777 else fHistSecMatMCnegBefEvSel[jpart]->Fill(ptMC);
10ba520e 778 if(evSel==1){
a1c5e4a1 779 if(signMC>0) fHistSecMatMCpos[jpart]->Fill(ptMC);
780 else fHistSecMatMCneg[jpart]->Fill(ptMC);
10ba520e 781 }
36be14b3 782 }
783 }
10ba520e 784 }
785 }
786
36be14b3 787
10ba520e 788 if(evSel==0)return; //event selection
1f75ac7d 789
790
36be14b3 791 //loop on tracks
792 for (Int_t iTrack=0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
10ba520e 793 isph=-999;
36be14b3 794 code=-999;
795 mfl=-999;
7c831837 796 uniqueID=-999;
10ba520e 797
36be14b3 798 track = (AliESDtrack*)fESD->GetTrack(iTrack);
799 if (!track) continue;
800
801 //track selection
c841b89f 802 Int_t countBinTrk=1;
803 TString label;
804
805 label="no selection";
806 fHistNTracks->Fill(countBinTrk);
807 if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
808 if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
809 fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
810 fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
811 fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
812 countBinTrk++;
813
36be14b3 814 status=track->GetStatus();
815 if((status&AliESDtrack::kITSpureSA)==0) continue; //its standalone
c841b89f 816
817 label="ITSsa";
818 fHistNTracks->Fill(countBinTrk);
819 if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
820 if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
821 fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
822 fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
823 fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
824 countBinTrk++;
825
36be14b3 826 if((status&AliESDtrack::kITSrefit)==0) continue; //its refit
c841b89f 827
828 label="ITSrefit";
829 fHistNTracks->Fill(countBinTrk);
830 if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
831 if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
832 fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
833 fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
834 fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
835 countBinTrk++;
836
029e1796 837 if(TMath::Abs(track->GetSign())<0.0001) continue; //no neutral particles
c841b89f 838
839 label="neutral particle";
840 fHistNTracks->Fill(countBinTrk);
841 if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
842 if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
843 fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
844 fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
845 fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
846 countBinTrk++;
10ba520e 847
36be14b3 848 //cluster in ITS
849 UInt_t clumap = track->GetITSClusterMap();
850 Int_t nSPD=0;
851 for(Int_t il=0; il<2; il++) if(TESTBIT(clumap,il)) nSPD++;
852 if(nSPD<fMinSPDPts) continue;
c841b89f 853
854 label="SPDcls";
855 fHistNTracks->Fill(countBinTrk);
856 if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
857 if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
858 fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
859 fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
860 fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
861 countBinTrk++;
862
146d18ee 863 Int_t nPtsForPid=0;
864 for(Int_t j=2;j<6;j++) if(TESTBIT(clumap,j)) nPtsForPid++;
865 if(nPtsForPid<fMinNdEdxSamples) continue; //at least 3 points on SSD/SDD
c841b89f 866
867 label="SDD+SSD cls";
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 //chisquare/nclusters
146d18ee 877 Int_t nclu=nSPD+nPtsForPid;
36be14b3 878 if(track->GetITSchi2()/nclu > fMaxChi2Clu) continue;
c841b89f 879
880 label="chi2/ncls";
881 fHistNTracks->Fill(countBinTrk);
882 if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
883 if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
884 fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
885 fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
886 fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
887 countBinTrk++;
888
36be14b3 889 //pseudorapidity and rapidity
dac69cb0 890 if(TMath::Abs(track->Eta()) > fEtaRange) continue;
c841b89f 891
892 label="eta";
893 fHistNTracks->Fill(countBinTrk);
894 if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
895 if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
896 fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
897 fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
898 fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
899 countBinTrk++;
900
36be14b3 901 //truncated mean
146d18ee 902 track->GetITSdEdxSamples(dedxLay);
903 Double_t dedx = CookdEdx(dedxLay);
36be14b3 904 if(dedx<0) continue;
36be14b3 905
c841b89f 906 label="de/dx<0";
907 fHistNTracks->Fill(countBinTrk);
908 if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
909 if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
910 fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
911 fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
912 fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
913 countBinTrk++;
914
36be14b3 915 Float_t pt = track->Pt();
916 Int_t theBin=-1;
917 for(Int_t m=0; m<kNbins; m++){
918 if(TMath::Abs(pt) > fPtBinLimits[m] && TMath::Abs(pt) < fPtBinLimits[m+1]){
919 theBin=m;
920 break;
921 }
922 }
923 track->GetImpactParameters(impactXY, impactZ);
10ba520e 924
36be14b3 925 //Filling Ntuple
926 //information from the MC kinematics
927 if(fMC){
10ba520e 928 if(track->GetLabel()<0)isph=-1;
36be14b3 929 if(track->GetLabel()>=0){
930 part = (TParticle*)stack->Particle(track->GetLabel());
931 pdgPart = part->GetPDG();
932 code = pdgPart->PdgCode();
10ba520e 933 if(stack->IsPhysicalPrimary(track->GetLabel())) isph=1;
36be14b3 934 else{
10ba520e 935 isph=0;
936 Int_t indexMoth=part->GetFirstMother();
937 if(indexMoth>=0){
938 TParticle* moth = stack->Particle(indexMoth);
939 Float_t codemoth = TMath::Abs(moth->GetPdgCode());
940 mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
941 }
7c831837 942 uniqueID = part->GetUniqueID();
10ba520e 943 }
944
945 //Filling DCA distribution with MC truth
946
947 if(theBin>=0 && theBin<kNbins){
948 if(isph==1){//primaries in MC
949 if(track->GetSign()>0){
950 if(TMath::Abs(code)==listcode[0]) fHistMCPrimDCAPosPi[theBin]->Fill(impactXY);
951 if(TMath::Abs(code)==listcode[1]) fHistMCPrimDCAPosK[theBin]->Fill(impactXY);
952 if(TMath::Abs(code)==listcode[2]) fHistMCPrimDCAPosP[theBin]->Fill(impactXY);
953 }else{
954 if(TMath::Abs(code)==listcode[0]) fHistMCPrimDCANegPi[theBin]->Fill(impactXY);
955 if(TMath::Abs(code)==listcode[1]) fHistMCPrimDCANegK[theBin]->Fill(impactXY);
956 if(TMath::Abs(code)==listcode[2]) fHistMCPrimDCANegP[theBin]->Fill(impactXY);
957 }
958 }
959
960 if(isph==0){//primaries in MC
7c831837 961 if(mfl==3 && uniqueID == kPDecay){
10ba520e 962 if(track->GetSign()>0){
963 if(TMath::Abs(code)==listcode[0]) fHistMCSecStDCAPosPi[theBin]->Fill(impactXY);
964 if(TMath::Abs(code)==listcode[1]) fHistMCSecStDCAPosK[theBin]->Fill(impactXY);
965 if(TMath::Abs(code)==listcode[2]) fHistMCSecStDCAPosP[theBin]->Fill(impactXY);
966 }else{
967 if(TMath::Abs(code)==listcode[0]) fHistMCSecStDCANegPi[theBin]->Fill(impactXY);
968 if(TMath::Abs(code)==listcode[1]) fHistMCSecStDCANegK[theBin]->Fill(impactXY);
969 if(TMath::Abs(code)==listcode[2]) fHistMCSecStDCANegP[theBin]->Fill(impactXY);
970 }
971 }else{
972 if(track->GetSign()>0){
973 if(TMath::Abs(code)==listcode[0]) fHistMCSecMatDCAPosPi[theBin]->Fill(impactXY);
974 if(TMath::Abs(code)==listcode[1]) fHistMCSecMatDCAPosK[theBin]->Fill(impactXY);
975 if(TMath::Abs(code)==listcode[2]) fHistMCSecMatDCAPosP[theBin]->Fill(impactXY);
976 }else{
977 if(TMath::Abs(code)==listcode[0]) fHistMCSecMatDCANegPi[theBin]->Fill(impactXY);
978 if(TMath::Abs(code)==listcode[1]) fHistMCSecMatDCANegK[theBin]->Fill(impactXY);
979 if(TMath::Abs(code)==listcode[2]) fHistMCSecMatDCANegP[theBin]->Fill(impactXY);
980 }
981 }
982 }
36be14b3 983 }
984 }
985 }
7c831837 986 Float_t xnt[19];
36be14b3 987 Int_t index=0;
988 xnt[index++]=(Float_t)track->GetP();
989 xnt[index++]=(Float_t)track->Pt();
990 xnt[index++]=(Float_t)dedx;
146d18ee 991 xnt[index++]=(Float_t)dedxLay[0];
992 xnt[index++]=(Float_t)dedxLay[1];
993 xnt[index++]=(Float_t)dedxLay[2];
994 xnt[index++]=(Float_t)dedxLay[3];
995 xnt[index++]=(Float_t)nclu;
996 xnt[index++]=(Float_t)nPtsForPid;
36be14b3 997 xnt[index++]=(Float_t)track->GetSign();
998 xnt[index++]=(Float_t)fESD->GetRunNumber();
999 xnt[index++]=(Float_t)track->Eta();
1000 xnt[index++]=(Float_t)impactXY;
1001 xnt[index++]=(Float_t)impactZ;
1002 xnt[index++]=(Float_t)isph;
1003 xnt[index++]=(Float_t)code;
c841b89f 1004 xnt[index++]=(Float_t)mfl;
7c831837 1005 xnt[index++]=(Float_t)uniqueID;
c841b89f 1006 xnt[index]=(Float_t)track->GetITSchi2()/nclu;
36be14b3 1007
1008 if(fFillNtuple) fNtupleNSigma->Fill(xnt);
1009
1010
1011
1012 //Compute y and bb
10ba520e 1013 Double_t y[4],bbtheo[4],logdiff[4];
36be14b3 1014 Float_t p=track->GetP();
1015 if(fMC && fSmearMC){
1016 dedx=fRandGener->Gaus(dedx,fSmeardEdx*dedx);
1017 p=fRandGener->Gaus(p,fSmearP*p);
1018 }
98bb3219 1019
1020 //Nsigma Method
1021 Float_t resodedx[4];
1022 for(Int_t ires=0;ires<4;ires++){
1023 resodedx[ires]=fITSPIDResponse->GetResolution(1,ires+1,kTRUE);
1024 }
1025
10ba520e 1026 for(Int_t i=0;i<4;i++){
36be14b3 1027 y[i] = Eta2y(pt,pdgmass[i],track->Eta());
7fcf3a38 1028 //bbtheo[i]=fITSPIDResponse->Bethe(p,pdgmass[i],kTRUE); //Pure PHOBOS BB
1029 bbtheo[i]=fITSPIDResponse->BetheITSsaHybrid(p,pdgmass[i]); //PHOBOS + polinomial at low pt (below beta*gamma = 0.76)
36be14b3 1030 logdiff[i]=TMath::Log(dedx) - TMath::Log(bbtheo[i]);
1031 }
98bb3219 1032
146d18ee 1033 Int_t resocls=(Int_t)nPtsForPid-1;
10ba520e 1034
1035 //NSigma Method, with asymmetric bands
10ba520e 1036 Int_t minPosMean=-1;
1037 for(Int_t isp=0; isp<3; isp++){
1038 if(dedx<bbtheo[0])continue;
1039 Double_t bb=bbtheo[isp];
1040 if(dedx<bb){
1041 Double_t bbdistance=TMath::Abs((bbtheo[isp]-bbtheo[isp-1])/2);
1042 Double_t nsigma=TMath::Abs((dedx-bb)/bbdistance);
1043 if(nsigma<1.)minPosMean=isp;
1044 }
1045 else{
1046 Double_t bbdistance=TMath::Abs((bbtheo[isp]-bbtheo[isp+1])/2);
1047 Double_t nsigma=TMath::Abs((dedx-bb)/bbdistance);
1048 if(nsigma<1.)minPosMean=isp;
1049 }
1050 }
1051 if(dedx<bbtheo[0] && TMath::Abs((dedx-bbtheo[0])/(resodedx[resocls]*bbtheo[0]))<2)minPosMean=0;
1052
1053 //NSigma method with simmetric bands
98bb3219 1054
36be14b3 1055 Double_t nsigmas[3];
1056 Double_t min=999999.;
1057 Int_t minPos=-1;
1058 for(Int_t isp=0; isp<3; isp++){
1059 Double_t bb=bbtheo[isp];
1060 nsigmas[isp]=TMath::Abs((dedx-bb)/(resodedx[resocls]*bb));
1061 if(nsigmas[isp]<min){
1062 min=nsigmas[isp];
1063 minPos=isp;
1064 }
98bb3219 1065 //Filling histos with nsigma separation
1066 if(track->GetSign()>0)fHistPosNSigmaSep[isp]->Fill(track->GetP(),((dedx-bb)/(resodedx[resocls]*bb)));
1067 else fHistNegNSigmaSep[isp]->Fill(track->GetP(),((dedx-bb)/(resodedx[resocls]*bb)));
36be14b3 1068 }
10ba520e 1069
1070 // y calculation
c92ab4a5 1071 Double_t yPartMean=-999;
1072 Double_t yPart=-999;
1073 if(minPosMean>-1) yPartMean=y[minPosMean];
1074 if(minPos>-1) yPart=y[minPos];
36be14b3 1075
34f309eb 1076 if(TMath::Abs(yPartMean)<fMaxY){
10ba520e 1077 //DCA distributions, before the DCA cuts, based on asymmetrinc nsigma approach
36be14b3 1078 if(theBin>=0 && theBin<kNbins){
1079 if(track->GetSign()>0){
10ba520e 1080 if(minPosMean==0) fHistDCAPosPi[theBin]->Fill(impactXY);
1081 else if(minPosMean==1) fHistDCAPosK[theBin]->Fill(impactXY);
1082 else if(minPosMean==2) fHistDCAPosP[theBin]->Fill(impactXY);
36be14b3 1083 }else{
10ba520e 1084 if(minPosMean==0) fHistDCANegPi[theBin]->Fill(impactXY);
1085 else if(minPosMean==1) fHistDCANegK[theBin]->Fill(impactXY);
1086 else if(minPosMean==2) fHistDCANegP[theBin]->Fill(impactXY);
36be14b3 1087 }
1088 }
1089 }
1090
1091 //DCA cut on xy and z
a1c5e4a1 1092 if(!DCAcut(impactXY,impactZ,pt)) continue;
c841b89f 1093
1094 label="DCA";
1095 fHistNTracks->Fill(countBinTrk);
1096 if(track->GetSign()>0)fHistNTracksPos->Fill(countBinTrk);
1097 if(track->GetSign()<0)fHistNTracksNeg->Fill(countBinTrk);
1098 fHistNTracks->GetXaxis()->SetBinLabel(fHistNTracks->FindBin(countBinTrk),label.Data());
1099 fHistNTracksPos->GetXaxis()->SetBinLabel(fHistNTracksPos->FindBin(countBinTrk),label.Data());
1100 fHistNTracksNeg->GetXaxis()->SetBinLabel(fHistNTracksNeg->FindBin(countBinTrk),label.Data());
1101 countBinTrk++;
16c730d4 1102
1103
1104 Int_t jpart=-1;
10ba520e 1105
1106 //Filling Histos for Reco Efficiency
1107 //information from the MC kinematics
1108 if(fMC){
1109 if(track->GetLabel()<0)isph=-1;
1110 if(track->GetLabel()>=0){
1111 part = (TParticle*)stack->Particle(track->GetLabel());
1112 pdgPart = part->GetPDG();
1113 code = pdgPart->PdgCode();
10ba520e 1114 for(Int_t j=0; j<3; j++){
1115 if(TMath::Abs(code)==listcode[j]){
1116 jpart=j;
1117 break;
1118 }
1119 }
16c730d4 1120 if(jpart>=0){
1121 if(pdgPart->Charge()>0) signMC=1;
1122 else signMC=-1;
1123 ptMC=part->Pt();
1124 if(stack->IsPhysicalPrimary(track->GetLabel())){
dc3531c1 1125 if(signMC>0){
1126 fHistPrimMCposReco[jpart]->Fill(ptMC);
1127 if(TMath::Abs(y[jpart])<fMaxY) fHistPrimMCposRecoEtaY[jpart]->Fill(ptMC);
1128 }
1129 else{
1130 fHistPrimMCnegReco[jpart]->Fill(ptMC);
1131 if(TMath::Abs(y[jpart])<fMaxY) fHistPrimMCnegRecoEtaY[jpart]->Fill(ptMC);
1132 }
16c730d4 1133 }else{
1134 Int_t indexMoth=part->GetFirstMother();
1135 if(indexMoth>=0){
1136 TParticle* moth = stack->Particle(indexMoth);
1137 Float_t codemoth = TMath::Abs(moth->GetPdgCode());
1138 mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
1139 }
7c831837 1140 uniqueID = part->GetUniqueID();
1141 if(mfl==3 && uniqueID == kPDecay){ // strangeness
a1c5e4a1 1142 if(signMC>0) fHistSecStrMCposReco[jpart]->Fill(ptMC);
1143 else fHistSecStrMCnegReco[jpart]->Fill(ptMC);
16c730d4 1144 }else{
a1c5e4a1 1145 if(signMC>0) fHistSecMatMCposReco[jpart]->Fill(ptMC);
1146 else fHistSecMatMCnegReco[jpart]->Fill(ptMC);
16c730d4 1147 }
10ba520e 1148 }
1149 }
1150 }
1151 }
1152
1153 //Nsigma histos with MC truth
1154
1155 //asymmetric approach
c92ab4a5 1156 if(TMath::Abs(yPartMean)<fMaxY && minPosMean>-1){
10ba520e 1157 //nsigma histos
1158 if(track->GetSign()>0) fHistPosNSigmaMean[minPosMean]->Fill(pt);
1159 else fHistNegNSigmaMean[minPosMean]->Fill(pt);
1160 if(fMC){
1161 //nsigma histos with MC truth on PID
1162 if(TMath::Abs(code)==listcode[minPosMean]){
1163 if(track->GetSign()>0) fHistPosNSigmaMCMean[minPosMean]->Fill(pt);
1164 else fHistNegNSigmaMCMean[minPosMean]->Fill(pt);
1165 }
1166 //nsigma histos with MC truth on IsPhysicalPrimary
1167 if(isph==1){
1168 if(track->GetSign()>0) fHistPosNSigmaPrimMean[minPosMean]->Fill(pt);
1169 else fHistNegNSigmaPrimMean[minPosMean]->Fill(pt);
1170 //nsigma histos with MC truth on IsPhysicalPrimary and PID
1171 if(TMath::Abs(code)==listcode[minPosMean]){
1172 if(track->GetSign()>0) fHistPosNSigmaPrimMCMean[minPosMean]->Fill(pt);
1173 else fHistNegNSigmaPrimMCMean[minPosMean]->Fill(pt);
1174 }
1175 }
1176 }
1177 }
36be14b3 1178
10ba520e 1179 //symmetric bands
34f309eb 1180 if(min<fMinNSigma && TMath::Abs(yPart)<fMaxY){
10ba520e 1181 //nsigma histos
36be14b3 1182 if(track->GetSign()>0) fHistPosNSigma[minPos]->Fill(pt);
1183 else fHistNegNSigma[minPos]->Fill(pt);
1184 if(fMC){
10ba520e 1185 //nsigma histos with MC truth on PID
1186 if(TMath::Abs(code)==listcode[minPos]){
1187 if(track->GetSign()>0) fHistPosNSigmaMC[minPos]->Fill(pt);
1188 else fHistNegNSigmaMC[minPos]->Fill(pt);
1189 }
1190 //nsigma histos with MC truth on IsPhysicalPrimary
36be14b3 1191 if(isph==1){
1192 if(track->GetSign()>0) fHistPosNSigmaPrim[minPos]->Fill(pt);
1193 else fHistNegNSigmaPrim[minPos]->Fill(pt);
16c730d4 1194 //nsigma histos with MC truth on IsPhysicalPrimary and PID
36be14b3 1195 if(TMath::Abs(code)==listcode[minPos]){
1196 if(track->GetSign()>0) fHistPosNSigmaPrimMC[minPos]->Fill(pt);
1197 else fHistNegNSigmaPrimMC[minPos]->Fill(pt);
1198 }
1199 }
1200 }
1201 }
1202
10ba520e 1203
1204 //integral approach histograms
36be14b3 1205 if(theBin>=0 && theBin<kNbins){
1206 if(track->GetSign()>0){
1207 if(TMath::Abs(y[0]) < fMaxY)fHistPosPi[theBin]->Fill(logdiff[0]);
1208 if(TMath::Abs(y[1]) < fMaxY)fHistPosK[theBin]->Fill(logdiff[1]);
1209 if(TMath::Abs(y[2]) < fMaxY)fHistPosP[theBin]->Fill(logdiff[2]);
1210 if(fMC){
e6076bb8 1211 if(TMath::Abs(y[0])<fMaxY){
16c730d4 1212 if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypPion[theBin]->Fill(logdiff[0]);
1213 if(TMath::Abs(code)==11)fHistMCPosElHypPion[theBin]->Fill(logdiff[0]);
e6076bb8 1214 if(TMath::Abs(code)==211)fHistMCPosPiHypPion[theBin]->Fill(logdiff[0]);
98bb3219 1215 if(TMath::Abs(code)==321)fHistMCPosKHypPion[theBin]->Fill(logdiff[0]);
e6076bb8 1216 if(TMath::Abs(code)==2212)fHistMCPosPHypPion[theBin]->Fill(logdiff[0]);
1217 }
1218 if(TMath::Abs(y[1])<fMaxY){
16c730d4 1219 if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypKaon[theBin]->Fill(logdiff[1]);
1220 if(TMath::Abs(code)==11)fHistMCPosElHypKaon[theBin]->Fill(logdiff[1]);
e6076bb8 1221 if(TMath::Abs(code)==211)fHistMCPosPiHypKaon[theBin]->Fill(logdiff[1]);
98bb3219 1222 if(TMath::Abs(code)==321)fHistMCPosKHypKaon[theBin]->Fill(logdiff[1]);
e6076bb8 1223 if(TMath::Abs(code)==2212)fHistMCPosPHypKaon[theBin]->Fill(logdiff[1]);
1224 }
1225 if(TMath::Abs(y[2])<fMaxY){
16c730d4 1226 if(TMath::Abs(code)!=11 && jpart<0)fHistMCPosOtherHypProton[theBin]->Fill(logdiff[2]);
1227 if(TMath::Abs(code)==11)fHistMCPosElHypProton[theBin]->Fill(logdiff[2]);
1228 if(TMath::Abs(code)==211)fHistMCPosPiHypProton[theBin]->Fill(logdiff[2]);
1229 if(TMath::Abs(code)==321)fHistMCPosKHypProton[theBin]->Fill(logdiff[2]);
1230 if(TMath::Abs(code)==2212)fHistMCPosPHypProton[theBin]->Fill(logdiff[2]);
e6076bb8 1231 }
36be14b3 1232 }
1233 }else{
1234 if(TMath::Abs(y[0]) < fMaxY)fHistNegPi[theBin]->Fill(logdiff[0]);
1235 if(TMath::Abs(y[1]) < fMaxY)fHistNegK[theBin]->Fill(logdiff[1]);
1236 if(TMath::Abs(y[2]) < fMaxY)fHistNegP[theBin]->Fill(logdiff[2]);
1237 if(fMC){
e6076bb8 1238 if(TMath::Abs(y[0])<fMaxY){
16c730d4 1239 if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypPion[theBin]->Fill(logdiff[0]);
1240 if(TMath::Abs(code)==11)fHistMCNegElHypPion[theBin]->Fill(logdiff[0]);
98bb3219 1241 if(TMath::Abs(code)==211)fHistMCNegPiHypPion[theBin]->Fill(logdiff[0]);
1242 if(TMath::Abs(code)==321)fHistMCNegKHypPion[theBin]->Fill(logdiff[0]);
1243 if(TMath::Abs(code)==2212)fHistMCNegPHypPion[theBin]->Fill(logdiff[0]);
e6076bb8 1244 }
1245 if(TMath::Abs(y[1])<fMaxY){
16c730d4 1246 if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypKaon[theBin]->Fill(logdiff[1]);
1247 if(TMath::Abs(code)==11)fHistMCNegElHypKaon[theBin]->Fill(logdiff[1]);
98bb3219 1248 if(TMath::Abs(code)==211)fHistMCNegPiHypKaon[theBin]->Fill(logdiff[1]);
1249 if(TMath::Abs(code)==321)fHistMCNegKHypKaon[theBin]->Fill(logdiff[1]);
1250 if(TMath::Abs(code)==2212)fHistMCNegPHypKaon[theBin]->Fill(logdiff[1]);
e6076bb8 1251 }
1252 if(TMath::Abs(y[2])<fMaxY){
16c730d4 1253 if(TMath::Abs(code)!=11 && jpart<0)fHistMCNegOtherHypProton[theBin]->Fill(logdiff[2]);
1254 if(TMath::Abs(code)==11)fHistMCNegElHypProton[theBin]->Fill(logdiff[2]);
98bb3219 1255 if(TMath::Abs(code)==211)fHistMCNegPiHypProton[theBin]->Fill(logdiff[2]);
1256 if(TMath::Abs(code)==321)fHistMCNegKHypProton[theBin]->Fill(logdiff[2]);
1257 if(TMath::Abs(code)==2212)fHistMCNegPHypProton[theBin]->Fill(logdiff[2]);
e6076bb8 1258 }
36be14b3 1259 }
1260 }
1261 }
1262
1263
1264 //fill propaganda plot with dedx
1265 fHistDEDX->Fill(track->GetP(),dedx);
1266 fHistDEDXdouble->Fill(track->GetP()*track->GetSign(),dedx);
1267
1268 //fill charge distribution histo to check the calibration
1269 for(Int_t j=0;j<4;j++){
146d18ee 1270 if(dedxLay[j]<5) continue;
1271 fHistCharge[j]->Fill(dedxLay[j]);
36be14b3 1272 }
1273 }
1274
1275 // Post output data.
1276 PostData(1,fOutput);
a1c5e4a1 1277 PostData(2,fListCuts);
1278 // Printf("............. end of Exec");
36be14b3 1279}
1280
1281//________________________________________________________________________
1282void AliAnalysisTaskSEITSsaSpectra::Terminate(Option_t *) {
4d668cf7 1283 // Merge output
1284 // Called once at the end of the query
36be14b3 1285
1286 fOutput = dynamic_cast<TList*>(GetOutputData(1));
1287 if (!fOutput) {
1288 printf("ERROR: fOutput not available\n");
1289 return;
1290 }
1291 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
a1c5e4a1 1292 printf("Number of Analyzed Events = %f\n",fHistNEvents->GetBinContent(1));
36be14b3 1293
1294 Printf("end of Terminate");
1295 return;
0763e995 1296}