]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/SPECTRA/AliAnalysisTaskSEITSsaSpectra.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
All formulas now work for any case of the overlap between POI1, POI2 and RP in the...
[u/mrichter/AliRoot.git] / PWG2 / SPECTRA / 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>
28#include <TRandom3.h>
29#include <TH2F.h>
30#include <TChain.h>
31#include <TNtuple.h>
32#include <TParticle.h>
33#include <Rtypes.h>
34#include "AliAnalysisTaskSE.h"
35#include "AliAnalysisManager.h"
36#include "AliAnalysisDataContainer.h"
37#include "AliESDEvent.h"
38#include "AliESDInputHandler.h"
39#include "AliESDtrack.h"
40#include "AliStack.h"
41#include "AliMCEventHandler.h"
42#include "AliMCEvent.h"
43#include "AliPhysicsSelection.h"
44#include "AliAnalysisTaskSEITSsaSpectra.h"
45
46ClassImp(AliAnalysisTaskSEITSsaSpectra)
47
48//________________________________________________________________________
49AliAnalysisTaskSEITSsaSpectra::AliAnalysisTaskSEITSsaSpectra():
50AliAnalysisTaskSE("Task CFits"),
51 fESD(0),
52 fOutput(0),
53 fHistNEvents(0),
54 fHistNTracks(0),
55 fHistDEDX(0),
56 fHistDEDXdouble(0),
57 fHistBeforeEvSel(0),
58 fHistAfterEvSel(0),
59 fMinSPDPts(1),
60 fMinNdEdxSamples(3),
61 fMindEdx(0.),
62 fMinNSigma(3.),
63 fMaxY(0.5),
64 fMaxChi2Clu(1.),
65 fNSigmaDCAxy(7.),
66 fNSigmaDCAz(7.),
67 fMC(kFALSE),
68 fSmearMC(kFALSE),
69 fSmearP(0.),
70 fSmeardEdx(0.),
71 fRandGener(0),
72 fFillNtuple(kFALSE),
73 fNtupleNSigma(0),
74 fNtupleMC(0)
75{
76 // Constructor
77 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};
78 for(Int_t iBin=0; iBin<kNbins+1; iBin++) fPtBinLimits[iBin]=xbins[iBin];
79 fRandGener=new TRandom3(0);
80
81 DefineInput(0, TChain::Class());
82 DefineOutput(1, TList::Class());
83 Printf("end of AliAnalysisTaskSEITSsaSpectra");
84}
85
86//___________________________________________________________________________
87AliAnalysisTaskSEITSsaSpectra::~AliAnalysisTaskSEITSsaSpectra(){
88 // Destructor
89 if (fOutput) {
90 delete fOutput;
91 fOutput = 0;
92 }
93 if(fRandGener) delete fRandGener;
94}
95
96//________________________________________________________________________
97Double_t AliAnalysisTaskSEITSsaSpectra::CookdEdx(Double_t *s){
98 // truncated mean for the dEdx
99 Int_t nc=0;
100 Double_t dedx[4]={0.,0.,0.,0.};
101 for (Int_t il=0; il<4; il++) { // count good (>0) dE/dx values
102 if(s[il]>fMindEdx){
103 dedx[nc]= s[il];
104 nc++;
105 }
106 }
107 if(nc<fMinNdEdxSamples) return -1.;
108
109 Double_t tmp;
110 Int_t swap; // sort in ascending order
111 do {
112 swap=0;
113 for (Int_t i=0; i<nc-1; i++) {
114 if (dedx[i]<=dedx[i+1]) continue;
115 tmp=dedx[i];
116 dedx[i]=dedx[i+1];
117 dedx[i+1]=tmp;
118 swap++;
119 }
120 } while (swap);
121
122 Double_t sumamp=0,sumweight=0;
123 Double_t weight[4]={1.,1.,0.,0.};
124 if(nc==3) weight[1]=0.5;
125 else if(nc<3) weight[1]=0.;
126 for (Int_t i=0; i<nc; i++) {
127 sumamp+= dedx[i]*weight[i];
128 sumweight+=weight[i];
129 }
130 return sumamp/sumweight;
131}
132
133
134//________________________________________________________________________
135Bool_t AliAnalysisTaskSEITSsaSpectra::DCAcut(Double_t impactXY, Double_t impactZ, Double_t pt, Bool_t optMC){
136 // cut on transverse impact parameter updaated on 20-5-2010
137 // from the study of L. Milano, F. Prino on the ITS standalone tracks
138 // using the common binning of the TPC tracks
139 Double_t xyP[3];
140 Double_t zP[3];
141 if(optMC){
142 xyP[0]=88.63;//SIMpass4a12
143 xyP[1]=19.57;
144 xyP[2]=1.65;
145 zP[0]=140.98;
146 zP[1]=62.33;
147 zP[2]=1.15;
148 }
149 else{
150 xyP[0]=85.28;//DATApass6
151 xyP[1]=25.78;
152 xyP[2]=1.55;
153 zP[0]=146.80;
154 zP[1]=70.07;
155 zP[2]=1.11;
156 }
157 Double_t xySigma = xyP[0] + xyP[1]/TMath::Power(TMath::Abs(pt),xyP[2]);
158 Double_t xyMax = fNSigmaDCAxy*xySigma; //in micron
159 if((TMath::Abs(impactXY)*10000)>xyMax) return kFALSE;
160
161 Double_t zSigma = zP[0] + zP[1]/TMath::Power(TMath::Abs(pt),zP[2]);
162 Double_t zMax = fNSigmaDCAz*zSigma; //in micron
163 if((TMath::Abs(impactZ)*10000)>zMax) return kFALSE;
164
165 return kTRUE;
166}
167
168//________________________________________________________________________
169Double_t AliAnalysisTaskSEITSsaSpectra::Eta2y(Double_t pt, Double_t m, Double_t eta){
170 Double_t mt = TMath::Sqrt(m*m + pt*pt);
171 return TMath::ASinH(pt/mt*TMath::SinH(eta));
172}
173
174
175//________________________________________________________________________
176Double_t AliAnalysisTaskSEITSsaSpectra::BetheBloch(Double_t bg,Bool_t optMC) {
177 // BB PHOBOS parametrization tuned by G. Ortona on 900 GeV pp data of 2009
178 Double_t par[5];
179 if(optMC){//Double_t par[5]={139.1,23.36,0.06052,0.2043,-0.0004999};
180 par[0]=139.1;
181 par[1]=23.36;
182 par[2]=0.06052;
183 par[3]=0.2043;
184 par[4]=-0.0004999;
185 }else {
186 //Double_t par[5]={5.33458e+04,1.65303e+01,2.60065e-03,3.59533e-04,7.51168e-05};}
187 par[0]=5.33458e+04;
188 par[1]=1.65303e+01;
189 par[2]=2.60065e-03;
190 par[3]=3.59533e-04;
191 par[4]=7.51168e-05;
192 }
193 Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
194 Double_t gamma=bg/beta;
195 Double_t eff=1.0;
196 if(bg<par[2]) eff=(bg-par[3])*(bg-par[3])+par[4];
197 else eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
198 return (par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
199}
200
201
202//________________________________________________________________________
203void AliAnalysisTaskSEITSsaSpectra::UserCreateOutputObjects(){
029e1796 204 // Create a TList with histograms and a TNtuple
205 // Called once
206
36be14b3 207 fOutput = new TList();
208 fOutput->SetOwner();
209 fOutput->SetName("Spiderman");
210
211 fHistNEvents = new TH1F("fHistNEvents", "Number of processed events",6,0.5,6.5);
212 fHistNEvents->Sumw2();
213 fHistNEvents->SetMinimum(0);
214 fOutput->Add(fHistNEvents);
215
216 fHistNTracks = new TH1F("fHistNTracks", "Number of ITSsa tracks",20,0.5,20.5);
217 fHistNTracks->Sumw2();
218 fHistNTracks->SetMinimum(0);
219 fOutput->Add(fHistNTracks);
220
221 //binning for the histogram
222 const Int_t hnbins=400;
223 Double_t hxmin = 0.01;
224 Double_t hxmax = 10;
225 Double_t hlogxmin = TMath::Log10(hxmin);
226 Double_t hlogxmax = TMath::Log10(hxmax);
227 Double_t hbinwidth = (hlogxmax-hlogxmin)/hnbins;
228 Double_t hxbins[hnbins+1];
229 hxbins[0] = 0.01;
230 for (Int_t i=1;i<=hnbins;i++) {
231 hxbins[i] = hxmin + TMath::Power(10,hlogxmin+i*hbinwidth);
232 }
233
234 fHistDEDX = new TH2F("fHistDEDX","",hnbins,hxbins,900,0,1000);
235 fOutput->Add(fHistDEDX);
236
237 fHistDEDXdouble = new TH2F("fHistDEDXdouble","",500,-5,5,900,0,1000);
238 fOutput->Add(fHistDEDXdouble);
239
240
241 fHistBeforeEvSel = new TH1F("fHistBeforeEvSel","fHistBeforeEvSel",kNbins,fPtBinLimits);
242 fHistAfterEvSel = new TH1F("fHistAfterEvSel","fHistAfterEvSel",kNbins,fPtBinLimits);
243 fOutput->Add(fHistBeforeEvSel);
244 fOutput->Add(fHistAfterEvSel);
245
246
247
248 for(Int_t j=0;j<3;j++){
249 fHistMCpos[j] = new TH1F(Form("fHistMCpos%d",j),Form("fHistMCpos%d",j),kNbins,fPtBinLimits);
250 fHistMCneg[j] = new TH1F(Form("fHistMCneg%d",j),Form("fHistMCneg%d",j),kNbins,fPtBinLimits);
251 fOutput->Add(fHistMCneg[j]);
252 fOutput->Add(fHistMCpos[j]);
253 }
254
255
256 for(Int_t j=0;j<3;j++){
257 fHistMCposBefEvSel[j] = new TH1F(Form("fHistMCposBefEvSel%d",j),Form("fHistMCposBefEvSel%d",j),kNbins,fPtBinLimits);
258 fHistMCnegBefEvSel[j] = new TH1F(Form("fHistMCnegBefEvSel%d",j),Form("fHistMCnegBefEvSel%d",j),kNbins,fPtBinLimits);
259 fOutput->Add(fHistMCnegBefEvSel[j]);
260 fOutput->Add(fHistMCposBefEvSel[j]);
261 }
262
263 for(Int_t i=0; i<4; i++){
264 fHistCharge[i] = new TH1F(Form("fHistChargeLay%d",i),Form("fHistChargeLay%d",i),100,0,300);
265 fOutput->Add(fHistCharge[i]);
266 }
267
268 for(Int_t i=0; i<kNbins; i++){
269 fHistPosPi[i] = new TH1F(Form("fHistPosPi%d",i),Form("fHistPosPi%d",i),175,-3.5,3.5);
270 fHistPosK[i] = new TH1F(Form("fHistPosK%d",i),Form("fHistPosK%d",i),175,-3.5,3.5);
271 fHistPosP[i] = new TH1F(Form("fHistPosP%d",i),Form("fHistPosP%d",i),175,-3.5,3.5);
272 fHistNegPi[i] = new TH1F(Form("fHistNegPi%d",i),Form("fHistNegPi%d",i),175,-3.5,3.5);
273 fHistNegK[i] = new TH1F(Form("fHistNegK%d",i),Form("fHistNegK%d",i),175,-3.5,3.5);
274 fHistNegP[i] = new TH1F(Form("fHistNegP%d",i),Form("fHistNegP%d",i),175,-3.5,3.5);
275
276 fHistDCAPosPi[i] = new TH1F(Form("fHistDCAPosPi%d",i),Form("fHistDCAPosPi%d",i),2000,-1,1); //DCA distr.
277 fHistDCAPosK[i] = new TH1F(Form("fHistDCAPosK%d",i),Form("fHistDCAPosK%d",i),2000,-1,1);
278 fHistDCAPosP[i] = new TH1F(Form("fHistDCAPosP%d",i),Form("fHistDCAPosP%d",i),2000,-1,1);
279 fHistDCANegPi[i] = new TH1F(Form("fHistDCANegPi%d",i),Form("fHistDCANegPi%d",i),2000,-1,1);
280 fHistDCANegK[i] = new TH1F(Form("fHistDCANegK%d",i),Form("fHistDCANegK%d",i),2000,-1,1);
281 fHistDCANegP[i] = new TH1F(Form("fHistDCANegP%d",i),Form("fHistDCANegP%d",i),2000,-1,1);
282
283 fHistMCPosPi[i] = new TH1F(Form("fHistMCPosPi%d",i),Form("fHistMCPosPi%d",i),175,-3.5,3.5); //MC truth
284 fHistMCPosK[i] = new TH1F(Form("fHistMCPosK%d",i),Form("fHistMCPosK%d",i),175,-3.5,3.5);
285 fHistMCPosP[i] = new TH1F(Form("fHistMCPosP%d",i),Form("fHistMCPosP%d",i),175,-3.5,3.5);
286 fHistMCNegPi[i] = new TH1F(Form("fHistMCNegPi%d",i),Form("fHistMCNegPi%d",i),175,-3.5,3.5);
287 fHistMCNegK[i] = new TH1F(Form("fHistMCNegK%d",i),Form("fHistMCNegK%d",i),175,-3.5,3.5);
288 fHistMCNegP[i] = new TH1F(Form("fHistMCNegP%d",i),Form("fHistMCNegP%d",i),175,-3.5,3.5);
289 fOutput->Add(fHistPosPi[i]);
290 fOutput->Add(fHistPosK[i]);
291 fOutput->Add(fHistPosP[i]);
292 fOutput->Add(fHistNegPi[i]);
293 fOutput->Add(fHistNegK[i]);
294 fOutput->Add(fHistNegP[i]);
295
296 fOutput->Add(fHistDCAPosPi[i]);//DCA distr.
297 fOutput->Add(fHistDCAPosK[i]);
298 fOutput->Add(fHistDCAPosP[i]);
299 fOutput->Add(fHistDCANegPi[i]);
300 fOutput->Add(fHistDCANegK[i]);
301 fOutput->Add(fHistDCANegP[i]);
302
303 fOutput->Add(fHistMCPosPi[i]);//MC truth
304 fOutput->Add(fHistMCPosK[i]);
305 fOutput->Add(fHistMCPosP[i]);
306 fOutput->Add(fHistMCNegPi[i]);
307 fOutput->Add(fHistMCNegK[i]);
308 fOutput->Add(fHistMCNegP[i]);
309 }
310
311 //NSigma Histos
312 for(Int_t j=0;j<3;j++){
313 fHistPosNSigma[j] = new TH1F(Form("hHistPosNSigma%d",j),Form("hHistPosNSigma%d",j),kNbins,fPtBinLimits);
314 fHistNegNSigma[j] = new TH1F(Form("hHistNegNSigma%d",j),Form("hHistNegNSigma%d",j),kNbins,fPtBinLimits);
315 fHistPosNSigmaPrim[j] = new TH1F(Form("hHistPosNSigmaPrim%d",j),Form("hHistPosNSigmaPrim%d",j),kNbins,fPtBinLimits);
316 fHistNegNSigmaPrim[j] = new TH1F(Form("hHistNegNSigmaPrim%d",j),Form("hHistNegNSigmaPrim%d",j),kNbins,fPtBinLimits);
317 fHistPosNSigmaPrimMC[j] = new TH1F(Form("hHistPosNSigmaPrimMC%d",j),Form("hHistPosNSigmaPrimMC%d",j),kNbins,fPtBinLimits);
318 fHistNegNSigmaPrimMC[j] = new TH1F(Form("hHistNegNSigmaPrimMC%d",j),Form("hHistNegNSigmaPrimMC%d",j),kNbins,fPtBinLimits);
319 fOutput->Add(fHistPosNSigma[j]);
320 fOutput->Add(fHistNegNSigma[j]);
321 fOutput->Add(fHistPosNSigmaPrim[j]);
322 fOutput->Add(fHistNegNSigmaPrim[j]);
323 fOutput->Add(fHistPosNSigmaPrimMC[j]);
324 fOutput->Add(fHistNegNSigmaPrimMC[j]);
325 }
326
327 fNtupleNSigma = new TNtuple("fNtupleNSigma","fNtupleNSigma","p:pt:dedx:ncls:sign:run:eta:impactXY:impactZ:isph:pdgcode:mfl");
328 fOutput->Add(fNtupleNSigma);
329 fNtupleMC = new TNtuple("fNtupleMC","fNtupleMC","ptMC:pdgcode:signMC:etaMC:yMC:isph:evSel:run");
330 fOutput->Add(fNtupleMC);
331
332 Printf("end of CreateOutputObjects");
333}
334
335//________________________________________________________________________
336void AliAnalysisTaskSEITSsaSpectra::UserExec(Option_t *){
029e1796 337 // Main loop
338 // Called for each event
36be14b3 339
340 fESD=(AliESDEvent*)InputEvent();
341 if(!fESD) {
342 printf("AliAnalysisTaskSDDRP::Exec(): bad ESD\n");
343 return;
344 }
345
346 //binning for the dEdx distributions
347
348 //variables
349 Float_t pdgmass[3]={0.13957,0.493677,0.938272}; //mass for pi, K, P (Gev/c^2)
350 Int_t listcode[3]={211,321,2212};//code for pi, K, P (Gev/c^2)
351 Double_t s[4];
352 Float_t ptMC=-999,yMC=-999,code=-999, signMC=-999,isph=-999,mfl=-999.;
353 Float_t impactXY=-999, impactZ=-999;
354 Int_t evSel=1;
355 AliESDtrack* track;
356 UInt_t status;
357 AliStack* stack=0;
358 TParticle *part=0;
359 TParticlePDG *pdgPart=0;
360
361 //Nsigma Method
362 Float_t resodedx[4];
363 if(fMC){
364 resodedx[0]=0.13;
365 resodedx[1]=0.13;
366 resodedx[2]=0.134;
367 resodedx[3]=0.127;
368 }else{
369 resodedx[0]=0.23;
370 resodedx[1]=0.18;
371 resodedx[2]=0.16;
372 resodedx[3]=0.14;
373 }
374 /////////////////////
375 if(fMC){
376 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
377 if (!eventHandler) {
378 Printf("ERROR: Could not retrieve MC event handler");
379 return;
380 }
381 AliMCEvent* mcEvent = eventHandler->MCEvent();
382 if (!mcEvent) {
383 Printf("ERROR: Could not retrieve MC event");
384 return;
385 }
386 stack = mcEvent->Stack();
387 if (!stack) {
388 printf("ERROR: stack not available\n");
389 return;
390 }
391 }
392 //flags for MC
393 Int_t nTrackMC=0;
394 if(stack) nTrackMC = stack->GetNtrack();
395 const AliESDVertex *vtx = fESD->GetPrimaryVertexSPD();
396
397 //event selection
398 fHistNEvents->Fill(1);
399 if(!vtx)evSel=0;
400 else{
401 fHistNEvents->Fill(2);
402 if(vtx->GetNContributors()<0) evSel=0;
403 else{
404 fHistNEvents->Fill(3);
405 if(TMath::Abs(vtx->GetZv())>10) evSel=0;
406 else{
407 fHistNEvents->Fill(4);
408 if(vtx->GetZRes()>0.5) evSel=0;
409 else{
410 fHistNEvents->Fill(5);
411 if(vtx->IsFromVertexerZ() && vtx->GetDispersion()>0.03) evSel=0;
412 else fHistNEvents->Fill(6);
413 }
414 }
415 }
416 }
417
418 /////first loop on stack, before event selection, filling MC ntuple
419
420 for(Int_t imc=0; imc<nTrackMC; imc++){
421 part = stack->Particle(imc);
422 if(!stack->IsPhysicalPrimary(imc))continue;//no secondary in the MC sample
423 isph=1.;
424 pdgPart = part->GetPDG();
425 if(pdgPart->Charge()==0) continue; //no neutral particles
426 if(TMath::Abs(part->Eta()) > 0.9) continue; //pseudorapidity-acceptance cut
427 if(part->Energy() != TMath::Abs(part->Pz())) yMC = 0.5*TMath::Log((part->Energy()+part->Pz())/(part->Energy()-part->Pz()));
428 if(TMath::Abs(yMC) > fMaxY) continue; //rapidity cut
429
430 if(pdgPart->Charge()>0) signMC=1;
431 else signMC=-1;
432 ptMC=part->Pt();
433 code=pdgPart->PdgCode();
434
435
436 //filling MC ntuple
437 if(TMath::Abs(code)==211 || TMath::Abs(code)==321 || TMath::Abs(code)==2212){
438 Float_t xntMC[8];
439 Int_t indexMC=0;
440 xntMC[indexMC++]=(Float_t)ptMC;
441 xntMC[indexMC++]=(Float_t)code;
442 xntMC[indexMC++]=(Float_t)signMC;
443 xntMC[indexMC++]=(Float_t)part->Eta();
444 xntMC[indexMC++]=(Float_t)yMC;
445 xntMC[indexMC++]=(Float_t)isph;
446 xntMC[indexMC++]=(Float_t)evSel;
447 xntMC[indexMC++]=(Float_t)fESD->GetRunNumber();
448
449 if(fFillNtuple) fNtupleMC->Fill(xntMC);
450 }
451
452 for(Int_t j=0; j<3; j++){
453 if(TMath::Abs(code)==listcode[j]){
454 if(signMC>0) fHistMCposBefEvSel[j]->Fill(TMath::Abs(ptMC));
455 else fHistMCnegBefEvSel[j]->Fill(TMath::Abs(ptMC));
456 }
457 }
458 if(evSel==1){
459 for(Int_t j=0; j<3; j++){
460 if(TMath::Abs(code)==listcode[j]){
461 if(signMC>0) fHistMCpos[j]->Fill(TMath::Abs(ptMC));
462 else fHistMCneg[j]->Fill(TMath::Abs(ptMC));
463 }
464 }
465 }
466 }
467
468 if(evSel==0)return;
469
470 //loop on tracks
471 for (Int_t iTrack=0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
472 isph=-999.;
473 code=-999;
474 mfl=-999;
475
476 track = (AliESDtrack*)fESD->GetTrack(iTrack);
477 if (!track) continue;
478
479 //track selection
480 fHistNTracks->Fill(1);
481 status=track->GetStatus();
482 if((status&AliESDtrack::kITSpureSA)==0) continue; //its standalone
483 fHistNTracks->Fill(2);
484 if((status&AliESDtrack::kITSrefit)==0) continue; //its refit
485 fHistNTracks->Fill(3);
029e1796 486 if(TMath::Abs(track->GetSign())<0.0001) continue; //no neutral particles
36be14b3 487 fHistNTracks->Fill(4);
488
489
490 //cluster in ITS
491 UInt_t clumap = track->GetITSClusterMap();
492 Int_t nSPD=0;
493 for(Int_t il=0; il<2; il++) if(TESTBIT(clumap,il)) nSPD++;
494 if(nSPD<fMinSPDPts) continue;
495 fHistNTracks->Fill(5);
496 Int_t count=0;
497 for(Int_t j=2;j<6;j++) if(TESTBIT(clumap,j)) count++;
498 if(count<fMinNdEdxSamples) continue; //at least 3 points on SSD/SDD
499 fHistNTracks->Fill(6);
500 //chisquare/nclusters
501 Int_t nclu=nSPD+count;
502 if(track->GetITSchi2()/nclu > fMaxChi2Clu) continue;
503 fHistNTracks->Fill(7);
504 //pseudorapidity and rapidity
505 if(TMath::Abs(track->Eta()) > 0.9) continue;
506 fHistNTracks->Fill(8);
507 //truncated mean
508 //if(fMC) for(Int_t j=0;j<2;j++) s[j]*=3.34/5.43;//correction for SDD miscalibration of the MCpass4
509 track->GetITSdEdxSamples(s);
510 Double_t dedx = CookdEdx(s);
511 if(dedx<0) continue;
512 fHistNTracks->Fill(9);
513
514
515 Float_t pt = track->Pt();
516 Int_t theBin=-1;
517 for(Int_t m=0; m<kNbins; m++){
518 if(TMath::Abs(pt) > fPtBinLimits[m] && TMath::Abs(pt) < fPtBinLimits[m+1]){
519 theBin=m;
520 break;
521 }
522 }
523 track->GetImpactParameters(impactXY, impactZ);
524
525 //Filling Ntuple
526 //information from the MC kinematics
527 if(fMC){
528 if(track->GetLabel()<0)isph=-1.;
529 if(track->GetLabel()>=0){
530 part = (TParticle*)stack->Particle(track->GetLabel());
531 pdgPart = part->GetPDG();
532 code = pdgPart->PdgCode();
533 if(stack->IsPhysicalPrimary(track->GetLabel()))isph=1.;
534 else{
535 isph=0.;
536 TParticle* moth = stack->Particle(part->GetFirstMother());
537 Float_t codemoth = TMath::Abs(moth->GetPdgCode());
538 mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));
539 }
540 }
541 }
542 Float_t xnt[12];
543 Int_t index=0;
544 xnt[index++]=(Float_t)track->GetP();
545 xnt[index++]=(Float_t)track->Pt();
546 xnt[index++]=(Float_t)dedx;
547 xnt[index++]=(Float_t)count;
548 xnt[index++]=(Float_t)track->GetSign();
549 xnt[index++]=(Float_t)fESD->GetRunNumber();
550 xnt[index++]=(Float_t)track->Eta();
551 xnt[index++]=(Float_t)impactXY;
552 xnt[index++]=(Float_t)impactZ;
553 xnt[index++]=(Float_t)isph;
554 xnt[index++]=(Float_t)code;
555 xnt[index]=(Float_t)mfl;
556
557 if(fFillNtuple) fNtupleNSigma->Fill(xnt);
558
559
560
561 //Compute y and bb
562 Double_t y[3],bbtheo[3],logdiff[3];
563 Float_t p=track->GetP();
564 if(fMC && fSmearMC){
565 dedx=fRandGener->Gaus(dedx,fSmeardEdx*dedx);
566 p=fRandGener->Gaus(p,fSmearP*p);
567 }
568
569 for(Int_t i=0;i<3;i++){
570 y[i] = Eta2y(pt,pdgmass[i],track->Eta());
571 bbtheo[i]=BetheBloch(p/pdgmass[i],fMC);
572 logdiff[i]=TMath::Log(dedx) - TMath::Log(bbtheo[i]);
573 }
574
575 //NSigma Method
576 Int_t resocls=(Int_t)count-1;
577
578 Double_t nsigmas[3];
579 Double_t min=999999.;
580 Int_t minPos=-1;
581 for(Int_t isp=0; isp<3; isp++){
582 Double_t bb=bbtheo[isp];
583 nsigmas[isp]=TMath::Abs((dedx-bb)/(resodedx[resocls]*bb));
584 if(nsigmas[isp]<min){
585 min=nsigmas[isp];
586 minPos=isp;
587 }
588 }
589 Double_t yPart=y[minPos];
590
591 if(min<fMinNSigma && yPart<fMaxY){
592 //DCA distributions, before the DCA cuts
593 if(theBin>=0 && theBin<kNbins){
594 if(track->GetSign()>0){
595 if(minPos==0) fHistDCAPosPi[theBin]->Fill(impactXY);
596 else if(minPos==1) fHistDCAPosK[theBin]->Fill(impactXY);
597 else if(minPos==2) fHistDCAPosP[theBin]->Fill(impactXY);
598 }else{
599 if(minPos==0) fHistDCANegPi[theBin]->Fill(impactXY);
600 else if(minPos==1) fHistDCANegK[theBin]->Fill(impactXY);
601 else if(minPos==2) fHistDCANegP[theBin]->Fill(impactXY);
602 }
603 }
604 }
605
606 //DCA cut on xy and z
607 if(!DCAcut(impactXY,impactZ,pt,fMC)) continue;
608 fHistNTracks->Fill(10);
609
610 if(min<fMinNSigma && yPart<fMaxY){
611 if(track->GetSign()>0) fHistPosNSigma[minPos]->Fill(pt);
612 else fHistNegNSigma[minPos]->Fill(pt);
613 if(fMC){
614 if(isph==1){
615 if(track->GetSign()>0) fHistPosNSigmaPrim[minPos]->Fill(pt);
616 else fHistNegNSigmaPrim[minPos]->Fill(pt);
617 if(TMath::Abs(code)==listcode[minPos]){
618 if(track->GetSign()>0) fHistPosNSigmaPrimMC[minPos]->Fill(pt);
619 else fHistNegNSigmaPrimMC[minPos]->Fill(pt);
620 }
621 }
622 }
623 }
624
625 if(theBin>=0 && theBin<kNbins){
626 if(track->GetSign()>0){
627 if(TMath::Abs(y[0]) < fMaxY)fHistPosPi[theBin]->Fill(logdiff[0]);
628 if(TMath::Abs(y[1]) < fMaxY)fHistPosK[theBin]->Fill(logdiff[1]);
629 if(TMath::Abs(y[2]) < fMaxY)fHistPosP[theBin]->Fill(logdiff[2]);
630 if(fMC){
631 if((TMath::Abs(y[0])<fMaxY) && (TMath::Abs(code)==211)) fHistMCPosPi[theBin]->Fill(logdiff[0]);
632 if((TMath::Abs(y[1])<fMaxY) && (TMath::Abs(code)==321)) fHistMCPosK[theBin]->Fill(logdiff[1]);
633 if((TMath::Abs(y[2])<fMaxY) && (TMath::Abs(code)==2212)) fHistMCPosP[theBin]->Fill(logdiff[2]);
634 }
635 }else{
636 if(TMath::Abs(y[0]) < fMaxY)fHistNegPi[theBin]->Fill(logdiff[0]);
637 if(TMath::Abs(y[1]) < fMaxY)fHistNegK[theBin]->Fill(logdiff[1]);
638 if(TMath::Abs(y[2]) < fMaxY)fHistNegP[theBin]->Fill(logdiff[2]);
639 if(fMC){
640 if((TMath::Abs(y[0])<fMaxY) && (TMath::Abs(code)==211)) fHistMCNegPi[theBin]->Fill(logdiff[0]);
641 if((TMath::Abs(y[1])<fMaxY) && (TMath::Abs(code)==321)) fHistMCNegK[theBin]->Fill(logdiff[1]);
642 if((TMath::Abs(y[2])<fMaxY) && (TMath::Abs(code)==2212)) fHistMCNegP[theBin]->Fill(logdiff[2]);
643 }
644 }
645 }
646
647
648 //fill propaganda plot with dedx
649 fHistDEDX->Fill(track->GetP(),dedx);
650 fHistDEDXdouble->Fill(track->GetP()*track->GetSign(),dedx);
651
652 //fill charge distribution histo to check the calibration
653 for(Int_t j=0;j<4;j++){
654 if(s[j]<5) continue;
655 fHistCharge[j]->Fill(s[j]);
656 }
657 }
658
659 // Post output data.
660 PostData(1,fOutput);
661 Printf("............. end of Exec");
662}
663
664//________________________________________________________________________
665void AliAnalysisTaskSEITSsaSpectra::Terminate(Option_t *) {
029e1796 666 // Merge output
667 // Called once at the end of the query
36be14b3 668
669 fOutput = dynamic_cast<TList*>(GetOutputData(1));
670 if (!fOutput) {
671 printf("ERROR: fOutput not available\n");
672 return;
673 }
674 fHistNEvents = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNEvents"));
675 fHistNTracks = dynamic_cast<TH1F*>(fOutput->FindObject("fHistNTracks"));
676 fHistDEDX = dynamic_cast<TH2F*>(fOutput->FindObject("fHistDEDX"));
677 fHistDEDXdouble = dynamic_cast<TH2F*>(fOutput->FindObject("fHistDEDXdouble"));
678
679 fHistBeforeEvSel = dynamic_cast<TH1F*>(fOutput->FindObject("fHistBeforeEvSel"));
680 fHistAfterEvSel = dynamic_cast<TH1F*>(fOutput->FindObject("fHistAfterEvSel"));
681
682
683 for(Int_t j=0;j<3;j++){
684 fHistMCpos[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCpos%d",j)));
685 fHistMCneg[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCneg%d",j)));
686 }
687
688
689 for(Int_t j=0;j<3;j++){
690 fHistMCposBefEvSel[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCposBefEvSel%d",j)));
691 fHistMCnegBefEvSel[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCnegBefEvSel%d",j)));
692 }
693
694 for(Int_t i=0; i<4; i++){
695 fHistCharge[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistChargeLay%d",i)));
696 }
697
698 for(Int_t i=0; i<kNbins; i++){
699 fHistPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPosPi%d",i)));
700 fHistPosK[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPosK%d",i)));
701 fHistPosP[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistPosP%d",i)));
702 fHistNegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistNegPi%d",i)));
703 fHistNegK[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistNegK%d",i)));
704 fHistNegP[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistNegP%d",i)));
705
706 fHistDCAPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCAPosPi%d",i)));
707 fHistDCAPosK[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCAPosK%d",i)));
708 fHistDCAPosP[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCAPosP%d",i)));
709 fHistDCANegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCANegPi%d",i)));
710 fHistDCANegK[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCANegK%d",i)));
711 fHistDCANegP[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistDCANegP%d",i)));
712
713
714 if(fMC){
715 fHistMCPosPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosPi%d",i)));
716 fHistMCPosK[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosK%d",i)));
717 fHistMCPosP[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCPosP%d",i)));
718 fHistMCNegPi[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegPi%d",i)));
719 fHistMCNegK[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegK%d",i)));
720 fHistMCNegP[i] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("fHistMCNegP%d",i)));
721 }
722 }
723
724 for(Int_t j=0;j<3;j++){
725 fHistPosNSigma[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigma%d",j)));
726 fHistNegNSigma[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigma%d",j)));
727 fHistPosNSigmaPrim[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaPrim%d",j)));
728 fHistNegNSigmaPrim[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaPrim%d",j)));
729 fHistPosNSigmaPrimMC[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistPosNSigmaPrimMC%d",j)));
730 fHistNegNSigmaPrimMC[j] = dynamic_cast<TH1F*>(fOutput->FindObject(Form("hHistNegNSigmaPrimMC%d",j)));
731 }
732
733 fNtupleNSigma = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleNSigma"));
734 fNtupleMC = dynamic_cast<TNtuple*>(fOutput->FindObject("fNtupleMC"));
735
736 Printf("end of Terminate");
737 return;
738}