]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/EBYE/Fluctuations/AliEbyEHigherMomentsTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Major changes in the task. Shift from TH2D to THnSparse to avoid high memory consumption.
[u/mrichter/AliRoot.git] / PWGCF / EBYE / Fluctuations / AliEbyEHigherMomentsTask.cxx
CommitLineData
8768ec6a 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: Satyajit Jena. *
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//=========================================================================//
18// //
19// Analysis Task for Net-Charge Higher Moment Analysis //
20// Author: Satyajit Jena || Nirbhay K. Behera //
21// sjena@cern.ch || nbehera@cern.ch //
22// //
23//=========================================================================//
24
25#include "TChain.h"
26#include "TList.h"
27#include "TFile.h"
28#include "TTree.h"
29#include "TH1D.h"
e405901b 30#include "TH2D.h"
31#include "TH3D.h"
32#include "THnSparse.h"
8768ec6a 33#include "TCanvas.h"
34
35#include "AliAnalysisTaskSE.h"
36#include "AliAnalysisManager.h"
37
38#include "AliESDVertex.h"
39#include "AliESDEvent.h"
40#include "AliESDInputHandler.h"
41#include "AliAODEvent.h"
42#include "AliAODTrack.h"
43#include "AliAODInputHandler.h"
44#include "AliESD.h"
45#include "AliESDEvent.h"
46#include "AliAODEvent.h"
47#include "AliStack.h"
48#include "AliESDtrackCuts.h"
49#include "AliAODMCHeader.h"
50#include "AliAODMCParticle.h"
51#include "AliMCEventHandler.h"
52#include "AliMCEvent.h"
53#include "AliStack.h"
54#include "AliGenHijingEventHeader.h"
55#include "AliGenEventHeader.h"
56#include "AliPID.h"
57#include "AliAODPid.h"
58#include "AliPIDResponse.h"
59#include "AliAODpidUtil.h"
60#include "AliPIDCombined.h"
61
62#include "AliEbyEHigherMomentsTask.h"
63
64using std::cout;
65using std::endl;
66
67ClassImp(AliEbyEHigherMomentsTask)
68
69
70//-----------------------------------------------------------------------
71AliEbyEHigherMomentsTask::AliEbyEHigherMomentsTask( const char *name )
72: AliAnalysisTaskSE( name ),
e405901b 73 fListOfHistosQA(0),
8768ec6a 74 fListOfHistos(0),
e405901b 75 fAOD(0),
76 fArrayMC(0),
77 fPIDResponse(0),
78 fParticleSpecies(AliPID::kProton),
beb13b6c 79 fAnalysisType(0),
8768ec6a 80 fCentralityEstimator("V0M"),
e405901b 81 fCentrality(0),
8768ec6a 82 fVxMax(3.),
83 fVyMax(3.),
e405901b 84 fVzMax(10.),
8768ec6a 85 fDCAxy(2.4),
86 fDCAz(3.),
87 fPtLowerLimit(0.3),
88 fPtHigherLimit(1.5),
89 fEtaLowerLimit(-0.8),
90 fEtaHigherLimit(0.8),
e405901b 91 fRapidityCut(0.5),
92 fNSigmaCut(3.),
8768ec6a 93 fTPCNClus(80),
abe73236 94 fChi2perNDF(4.),
e405901b 95 fAODtrackCutBit(128),
96 fLabel(NULL),
97 fUsePid(kFALSE),
98 fCheckEff(kFALSE),
99 fEventCounter(0),
100 fHistDCA(0),
101 fTPCSig(0),
102 fTPCSigA(0),
103 fTHnCentNplusNminusCh(0),
104 fTHnCentNplusNminusChTruth(0),
105 fTHnCentNplusNminus(0),
106 fTHnEfficiencyHisto(0)
8768ec6a 107{
108
abe73236 109 for ( Int_t i = 0; i < 13; i++) {
8768ec6a 110 fHistQA[i] = NULL;
111 }
abe73236 112
e405901b 113 for ( Int_t i = 0; i < 5; i++ ){
114 fTHnCentNplusNminusPid[i] = NULL;
115 fTHnCentNplusNminusPidTruth[i] = NULL;
116 }
117
8768ec6a 118 DefineOutput(1, TList::Class()); // Outpput....
e405901b 119 DefineOutput(2, TList::Class());
8768ec6a 120}
121
122AliEbyEHigherMomentsTask::~AliEbyEHigherMomentsTask() {
e405901b 123 if(fListOfHistosQA) delete fListOfHistosQA;
8768ec6a 124 if(fListOfHistos) delete fListOfHistos;
e405901b 125 if(fLabel[0] ) delete [] (fLabel[0]);
126 if(fLabel[1] ) delete [] (fLabel[1]);
8768ec6a 127}
128
129//---------------------------------------------------------------------------------
130void AliEbyEHigherMomentsTask::UserCreateOutputObjects() {
e405901b 131
132 fListOfHistosQA = new TList();
133 fListOfHistosQA->SetOwner(kTRUE);
8768ec6a 134 fListOfHistos = new TList();
135 fListOfHistos->SetOwner(kTRUE);
8768ec6a 136
e405901b 137 fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
138 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
139 fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-90% centrality");
140 fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
141 fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
142 fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
143 fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
144 fListOfHistosQA->Add(fEventCounter);
145
146 //For QA-Histograms
8768ec6a 147 fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.);
148 fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.);
149 fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
150 fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.);
151 fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.);
152 fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.);
153 fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.);
154 fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.);
155 fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
156 fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
157 fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
abe73236 158 fHistQA[11] = new TH1D("fHistQANCls","Number of TPC cluster",200,0,200);
159 fHistQA[12] = new TH1D("fHistQAChi2","Chi2 per NDF",100,0,10);
abe73236 160 for(Int_t i = 0; i < 13; i++)
8768ec6a 161 {
e405901b 162 fListOfHistosQA->Add(fHistQA[i]);
8768ec6a 163 }
abe73236 164
e405901b 165 fHistDCA = new TH2D("fHistDCA","DCAxy Vs DCAz", 500, -5., 5., 500, -5., 5.);
166 fTPCSig = new TH2D("fTPCSig","TPC signal",200, 0.0, 10. ,1000,0.,1000);
167 fTPCSig->SetMarkerColor(kRed);
168 fTPCSigA = new TH2D("fTPCSigA","TPC signal all ",200, 0.0, 10. ,1000,0.,1000);
169 fListOfHistosQA->Add(fHistDCA);
170 fListOfHistosQA->Add(fTPCSig);
171 fListOfHistosQA->Add(fTPCSigA);
172
173 const Int_t nDim = 3;
174 const Int_t nPid = 5;
175 TString Species[nPid] = {"Electron","Muon","Pion","Kaon","Proton"};
8768ec6a 176
e405901b 177 Int_t fBins[nPid][nDim];
178 Double_t fMin[nPid][nDim];
179 Double_t fMax[nPid][nDim];
8768ec6a 180
e405901b 181 for( Int_t i = 0; i < nPid; i++ ){
182 for( Int_t j = 0; j < nDim; j++ ){
183 fBins[i][j] = 0;
184 fMin[i][j] = 0.;
185 fMax[i][j] = 0.;
186 }
187 }
8768ec6a 188
abe73236 189
e405901b 190 Int_t fBinsCh[nDim] = {100, 1500, 1500};
191 Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
192 Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
193
194 fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
195 fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
196 fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
197 fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
198 fListOfHistos->Add(fTHnCentNplusNminusCh);
199
200 if( fAnalysisType == "MCAOD" ){
201
202 fTHnCentNplusNminusChTruth = new THnSparseD("fTHnCentNplusNminusChTruth","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
203 fTHnCentNplusNminusChTruth->GetAxis(0)->SetTitle("Centrality");
204 fTHnCentNplusNminusChTruth->GetAxis(1)->SetTitle("Nplus");
205 fTHnCentNplusNminusChTruth->GetAxis(2)->SetTitle("Nminus");
206 fListOfHistos->Add(fTHnCentNplusNminusChTruth);
207
208 if(fCheckEff){
209
210 Int_t effBin[10] = { 100, 3, 2, 2, 49, 180, 100, 49, 180, 100};
211 Double_t effBinL[10] = { -0.5, -1.5, -0.5, -0.5, 0.1, -0.9, 0., 0.1, -0.9, 0.};
212 Double_t effBinH[10] = { 99.5, 1.5, 1.5, 1.5, 5.0, 0.9, 6.8, 5.0, 0.9, 6.8 };
213
214 fTHnEfficiencyHisto = new THnSparseD("fTHnEfficiencyHisto","Cent-charge-pid-pt-eta-phi", 10, effBin, effBinL, effBinH);
215 fTHnEfficiencyHisto->GetAxis(0)->SetTitle("Centrality");
216 fTHnEfficiencyHisto->GetAxis(1)->SetTitle("Charge");
217 fTHnEfficiencyHisto->GetAxis(2)->SetTitle("RecStatus");
218 fTHnEfficiencyHisto->GetAxis(3)->SetTitle("PidRecStatus");
219 fTHnEfficiencyHisto->GetAxis(4)->SetTitle("P_{T}MC");
220 fTHnEfficiencyHisto->GetAxis(5)->SetTitle("#eta_{MC}");
221 fTHnEfficiencyHisto->GetAxis(6)->SetTitle("#phi_{MC}");
222 fTHnEfficiencyHisto->GetAxis(7)->SetTitle("P_{T}Rec");
223 fTHnEfficiencyHisto->GetAxis(8)->SetTitle("#eta_{Rec}");
224 fTHnEfficiencyHisto->GetAxis(9)->SetTitle("#phi_{Rec}");
225 fListOfHistos->Add(fTHnEfficiencyHisto);
226 }//Efficiency Sparese---
227
228 }//MCAOD---
229
230 TString hname1, hname11;
231 TString htitle1, axisTitle1, axisTitle2;
232
233
234 if( fUsePid ){
235
236 Int_t gPid = (Int_t)fParticleSpecies;
237
238 if( gPid > 1 && gPid < 5 ){
239 //Pion----
240 fBins[2][0] = 100; fBins[2][1] = 1000; fBins[2][2] = 1000;
241 fMin[2][0] = -0.5; fMin[2][1] = -0.5; fMin[2][2] = -0.5;
242 fMax[2][0] = 99.5; fMax[2][1] = 999.5; fMax[2][2] = 999.5;
243 //Kaon------
244 fBins[3][0] = 100; fBins[3][1] = 600; fBins[3][2] = 600;
245 fMin[3][0] = -0.5; fMin[3][1] = -0.5; fMin[3][2] = -0.5;
246 fMax[3][0] = 99.5; fMax[3][1] = 599.5; fMax[3][2] = 599.5;
247 //Proton-----
248 fBins[4][0] = 100; fBins[4][1] = 400; fBins[4][2] = 400;
249 fMin[4][0] = -0.5; fMin[4][1] = -0.5; fMin[4][2] = -0.5;
250 fMax[4][0] = 99.5; fMax[4][1] = 399.5; fMax[4][2] = 399.5;
251
252 hname1 = "fCentNplusNminusPid"; hname1 +=gPid;
253 htitle1 = Species[gPid]; htitle1 +=" And Neg-"; htitle1 +=Species[gPid];
254 axisTitle1 ="Pos-"; axisTitle1 += Species[gPid];
255 axisTitle2 ="Neg-"; axisTitle2 += Species[gPid];
256
257 fTHnCentNplusNminusPid[gPid] = new THnSparseD(hname1.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
258
259 fTHnCentNplusNminusPid[gPid]->GetAxis(0)->SetTitle("Centrality");
260 fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
261 fTHnCentNplusNminusPid[gPid]->GetAxis(2)->SetTitle(axisTitle2.Data());
262
263 fListOfHistos->Add(fTHnCentNplusNminusPid[gPid]);
264
265 if( fAnalysisType == "MCAOD" ){
266
267 hname11 = "fCentNplusNminusPidTruth"; hname11 +=gPid;
268 fTHnCentNplusNminusPidTruth[gPid] = new THnSparseD(hname11.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
269
270 fTHnCentNplusNminusPidTruth[gPid]->GetAxis(0)->SetTitle("Centrality");
271 fTHnCentNplusNminusPidTruth[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
272 fTHnCentNplusNminusPidTruth[gPid]->GetAxis(2)->SetTitle(axisTitle2.Data());
273 fListOfHistos->Add(fTHnCentNplusNminusPidTruth[gPid]);
274
275 }//MCAOD-----
276 }//Pion, Koan and Proton-------
277 else{
278
279 Int_t fBinsX[nDim] = {100, 1500, 1500};
280 Double_t fMinX[nDim] = { -0.5, -0.5, -0.5 };
281 Double_t fMaxX[nDim] = { 99.5, 1499.5, 1499.5};
282
283 fTHnCentNplusNminus = new THnSparseD("fTHnCentNplusNminus","Cent-NplusChrg-NminusChrg", nDim, fBinsX, fMinX, fMaxX);
284 fTHnCentNplusNminus->GetAxis(0)->SetTitle("Centrality");
285 fTHnCentNplusNminus->GetAxis(1)->SetTitle("UnwantedPlus");
286 fTHnCentNplusNminus->GetAxis(2)->SetTitle("UnwantedMinus");
287 fListOfHistos->Add(fTHnCentNplusNminus);
288
289 }//Unwanted particle -(other than pi, K or proton)
290
291 }//fUsePid-------
292
293 PostData(1, fListOfHistosQA);
294 PostData(2, fListOfHistos);
295
abe73236 296
8768ec6a 297}
298
299
300//----------------------------------------------------------------------------------
301void AliEbyEHigherMomentsTask::UserExec( Option_t * ){
abe73236 302
8768ec6a 303 fEventCounter->Fill(1);
e405901b 304
305 if(fAnalysisType == "AOD") {
306
307 doAODEvent();
308
309 }//AOD--analysis-----
310
311 else if(fAnalysisType == "MCAOD") {
312
313 doMCAODEvent();
314
315 }
316
317 else return;
318
319 fEventCounter->Fill(8);
320
321 PostData(1, fListOfHistosQA);
322 PostData(2, fListOfHistos);
323
324}
325
326//--------------------------------------------------------------------------------------
327void AliEbyEHigherMomentsTask::doAODEvent(){
328
8768ec6a 329 Double_t nPlusCharge = 0.;
330 Double_t nMinusCharge = 0.;
e405901b 331 Double_t nParticle = 0.;
332 Double_t nAntiParticle = 0.;
333 Int_t gPid = 0;
8768ec6a 334
e405901b 335 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
336 if (!manager) {
337 cout<<"ERROR: Analysis manager not found."<<endl;
8768ec6a 338 return;
339 }
e405901b 340 //coneect to the inputHandler------------
341 AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
342 if (!inputHandler) {
343 cout<<"ERROR: Input handler not found."<<endl;
344 return;
345 }
346
347 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
348 if (!fAOD)
349 {
350 cout<< "ERROR 01: AOD not found " <<endl;
351 return;
352 }
8768ec6a 353
e405901b 354 fPIDResponse =inputHandler->GetPIDResponse();
8768ec6a 355
e405901b 356
357 if (!fPIDResponse){
358 AliFatal("This Task needs the PID response attached to the inputHandler");
359 return;
360 }
361
362 AliAODHeader *aodHeader = fAOD->GetHeader();
363
364 fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
365
366 if(fCentrality < 0 || fCentrality >= 91) return;
367
368
369
370 fEventCounter->Fill(2);
371
372 if(!ProperVertex()) return;
373
374 Int_t nTracks = fAOD->GetNumberOfTracks();
375
376 TExMap *trackMap = new TExMap();//Mapping matrix----
377
378 for(Int_t i = 0; i < nTracks; i++) {
8768ec6a 379
e405901b 380 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
abe73236 381
e405901b 382 if(!aodTrack) {
383 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
384 continue;
385 }
abe73236 386
e405901b 387 Double_t tpcSignalAll = aodTrack->GetTPCsignal();
388 fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
abe73236 389
e405901b 390 Int_t gID = aodTrack->GetID();
abe73236 391
e405901b 392 if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
393 }//1st track loop----
394
395 AliAODTrack* newAodTrack;
396
397 for( Int_t j = 0; j < nTracks; j++ )
398 {
399
400 AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
401
402 if(!aodTrack1) {
403 AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
404 continue;
405 }
406
407
408 if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
409
410 Int_t gID = aodTrack1->GetID();
411
412 //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
413 newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
414
415
416 Double_t pt = aodTrack1->Pt();
417 Double_t eta = aodTrack1->Eta();
418 Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
419 Double_t chi2ndf = aodTrack1->Chi2perNDF();
420
421 if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
422 if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
423 if( nclus < fTPCNClus ) continue;
424 if( chi2ndf > fChi2perNDF ) continue;
425
426 if( fAODtrackCutBit == 128 ){
427 //TPC only tracks----
428 Float_t dxy = 0., dz = 0.;
429 dxy = aodTrack1->DCA();
430 dz = aodTrack1->ZAtDCA();
431
432 if( fabs(dxy) > fDCAxy ) continue;
433 if( fabs(dz) > fDCAz ) continue;
434 //Extra cut on DCA---( Similar to BF Task.. )
435 if( fDCAxy !=0 && fDCAz !=0 ){
436 if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
437 }
438
439 fHistQA[6]->Fill(dxy);
440 fHistQA[7]->Fill(dz);
441 fHistDCA->Fill(dxy,dz);
442
443 }
444
445 fHistQA[8]->Fill(pt);
446 fHistQA[9]->Fill(eta);
447 fHistQA[10]->Fill(aodTrack1->Phi());
448 fHistQA[11]->Fill(nclus);
449 fHistQA[12]->Fill(chi2ndf);
450
451 Short_t gCharge = aodTrack1->Charge();
452
453 if(gCharge > 0) nPlusCharge += 1.;
454 if(gCharge < 0) nMinusCharge += 1.;
455
456 if( fUsePid ) {
457
458 gPid = (Int_t)fParticleSpecies;
459
460 Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
461 Double_t tpcSignal = newAodTrack->GetTPCsignal();
462 //Double_t rap = aodTrack1->Y(AliPID::ParticleMass(fParticleSpecies));
463 //Double_t tpcSignal = aodTrack1->GetTPCsignal();
464
465 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
466
467 fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
468
469 Float_t nsigmaTPCPID = -999.;
470 Float_t nsigmaTOFPID = -999.;
471 //Float_t nsigmaTPCTOFPID = -999.;
472
473 nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
474 nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
475
476 if ( nsigmaTPCPID < fNSigmaCut ){
abe73236 477
e405901b 478 if (gCharge > 0) nParticle +=1.;
479 if( gCharge < 0 ) nAntiParticle +=1.;
abe73236 480
e405901b 481 }
482 }//fUsepid-----
483
484 }//--------- Track Loop to select with filterbit
8768ec6a 485
e405901b 486
487 //cout << fCentrality <<" "<< nPlusCharge << " " << nMinusCharge << endl;
488 //cout << fCentrality <<" "<< nParticle << " " << nAntiParticle << endl;
489
490 Double_t fContainerCh[3] = { fCentrality, nPlusCharge, nMinusCharge};
491 Double_t fContainerPid[3] = { fCentrality, nParticle, nAntiParticle};
492
493
494 fTHnCentNplusNminusCh->Fill(fContainerCh);
495
496 if( fUsePid ){
497 gPid = (Int_t)fParticleSpecies;
498 fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
499 }
500
501 fEventCounter->Fill(7);
502
503 return;
504
505}
506//--------------------------------------------------------------------------------------
507//--------------------------------------------------------------------------------------
508void AliEbyEHigherMomentsTask::doMCAODEvent(){
509
510
511 //---------
512 Double_t nPlusCharge = 0.;
513 Double_t nMinusCharge = 0.;
514
515 Double_t nPlusChargeTruth = 0.;
516 Double_t nMinusChargeTruth = 0.;
517
518
519 Double_t nParticle = 0.;
520 Double_t nAntiParticle = 0.;
521 Double_t nParticleTruth = 0.;
522 Double_t nAntiParticleTruth = 0.;
523
524
525 Int_t gPid = 0;
526 Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
527
528 //Connect to Anlaysis manager------
529 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
530 if (!manager) {
531 cout<<"ERROR: Analysis manager not found."<<endl;
532 return;
533 }
534
535 AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
536 if (!inputHandler) {
537 cout<<"ERROR: Input handler not found."<<endl;
538 return;
539 }
540
541 //AOD event------
542 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
543 if (!fAOD)
544 {
545 cout<< "ERROR 01: AOD not found " <<endl;
beb13b6c 546 return;
8768ec6a 547 }
e405901b 548
549 fPIDResponse =inputHandler->GetPIDResponse();
550
551
552 if (!fPIDResponse){
553 AliFatal("This Task needs the PID response attached to the inputHandler");
554 return;
555 }
556
557 // -- Get MC header
558 // ------------------------------------------------------------------
559
560 fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
561 if (!fArrayMC) {
562 AliFatal("Error: MC particles branch not found!\n");
563 return;
564 }
565
566 AliAODMCHeader *mcHdr=NULL;
567 mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
568 if(!mcHdr) {
569 printf("MC header branch not found!\n");
570 return;
571 }
572
573 AliAODHeader *aodHeader = fAOD->GetHeader();
574
575 fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
576
577
578
579 if( fCentrality < 0 || fCentrality >= 91) return;
580
581 fEventCounter->Fill(2);
582
583 if(!ProperVertex()) return;
584
585 Int_t nTracks = fAOD->GetNumberOfTracks();
586
587 fLabel = new Int_t*[2];
588 fLabel[0] = new Int_t[nTracks]; //All charged hadrons----
589 fLabel[1] = new Int_t[nTracks]; //For Pid-----------
590 //Initialize labels----
591
592 if(!fLabel[0]){
593 AliError("Can't Get fLabel[0] ");
594 return;
595 }
596 if(!fLabel[1]){
597 AliError("Can't Get fLabel[1] ");
598 return;
599 }
600
601 for(Int_t i=0; i < 2; i++){
602 for(Int_t j=0; j < nTracks; j++){
603 fLabel[i][j] = 0;
8768ec6a 604 }
e405901b 605 }
606
607
608 TExMap *trackMap = new TExMap();//Mapping matrix----
609
610 for(Int_t i = 0; i < nTracks; i++) {
8768ec6a 611
e405901b 612 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
8768ec6a 613
e405901b 614 if(!aodTrack) {
615 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
616 continue;
617 }
8768ec6a 618
e405901b 619 Double_t tpcSignalAll = aodTrack->GetTPCsignal();
620 fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
621
622 Int_t gID = aodTrack->GetID();
623
624 if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global(primary) tracks-----
625
626 }//1st track loop----
627
628 AliAODTrack* newAodTrack;
629
630 for( Int_t j = 0; j < nTracks; j++ ){
631
632 AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
633
634 if(!aodTrack1) {
635 AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
636 continue;
abe73236 637 }
8768ec6a 638
e405901b 639
640 if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
641
642 Int_t gID = aodTrack1->GetID();
643
644 //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
645
646 newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
647
648
649 //cout << dxy << endl;
650 Double_t pt = aodTrack1->Pt();
651 Double_t eta = aodTrack1->Eta();
652 Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
653 Double_t chi2ndf = aodTrack1->Chi2perNDF();
654
655 if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
656 if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
657 if( nclus < fTPCNClus ) continue;
658 if( chi2ndf > fChi2perNDF ) continue;
659
660 if( fAODtrackCutBit == 128 ){
661 Float_t dxy = 0., dz = 0.;
662 dxy = aodTrack1->DCA();
663 dz = aodTrack1->ZAtDCA();
664
665 if( fabs(dxy) > fDCAxy ) continue;
666 if( fabs(dz) > fDCAz ) continue;
667 //Extra cut on DCA---( Similar to BF Task.. )
668 if( fDCAxy !=0 && fDCAz !=0 ){
669 if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
beb13b6c 670 }
8768ec6a 671
e405901b 672 fHistQA[6]->Fill(dxy);
673 fHistQA[7]->Fill(dz);
674 fHistDCA->Fill(dxy,dz);
abe73236 675
e405901b 676 }
677
678
679
680 fHistQA[8]->Fill(pt);
681 fHistQA[9]->Fill(eta);
682 fHistQA[10]->Fill(aodTrack1->Phi());
683 fHistQA[11]->Fill(nclus);
684 fHistQA[12]->Fill(chi2ndf);
685
686
687 Short_t gCharge = aodTrack1->Charge();
688
689 if( gCharge == 0 ) continue;
abe73236 690
e405901b 691 //Check the labels----------
692 Int_t label = TMath::Abs(aodTrack1->GetLabel());
693 //fill the labels--------
694 fLabel[0][j] = label;//charged particle---
695
696 AliAODMCParticle* particle = (AliAODMCParticle*)fArrayMC->At(label);
697 if (!particle) return;
698
699 if(gCharge > 0) nPlusCharge += 1.;
700 if(gCharge < 0) nMinusCharge += 1.;
701
702 if( fUsePid ) {
703
704 gPid = (Int_t)fParticleSpecies;
705 Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
706 Double_t tpcSignal = newAodTrack->GetTPCsignal();
707
708 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
709
710 fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
711
712 Float_t nsigmaTPCPID = -999.;
713 Float_t nsigmaTOFPID = -999.;
714 //Float_t nsigmaTPCTOFPID = -999.;
715
716 nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
717 nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
718
719 if( nsigmaTPCPID < fNSigmaCut ){
720 fLabel[1][j] = label;//pid labels----
721 if (gCharge > 0) nParticle +=1;
722 if( gCharge < 0 ) nAntiParticle +=1.;
723 }//nSigmaCut-----
724 }//fUsepid-----
725
726 }//--------- Track Loop to select with filterbit
abe73236 727
e405901b 728
729 Double_t fContainerCh[3] = { fCentrality, nPlusCharge, nMinusCharge};//Reco. values ch. hadrons
730 Double_t fContainerPid[3] = { fCentrality, nParticle, nAntiParticle};//Reco. values pid.
731
732 fTHnCentNplusNminusCh->Fill(fContainerCh);//Fill the ch. particles---
733 if( fUsePid ){
734 gPid = (Int_t)fParticleSpecies;
735 fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);//Fill the pid
736 }
737
738
739 //cout << fCentrality << " "<< nPlusCharge << " "<< nMinusCharge << endl;
740 //cout << nParticle <<" And " << nAntiParticle << endl;
741 //===========================================
742 //--------MC Truth---------------------------
743 //===========================================
744
745 Int_t nMCTrack = fArrayMC->GetEntriesFast();
746
747 for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
beb13b6c 748
e405901b 749 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
750
751 if(!partMC){
752 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
753 continue;
754 }
755
756 if(partMC->Charge() == 0) continue;
757 if(!partMC->IsPhysicalPrimary()) continue;
758
759 if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
760 if (partMC->Pt() < fPtLowerLimit || partMC->Pt() > fPtHigherLimit) continue;
761
762 Short_t gCharge = partMC->Charge();
763 Int_t chargeState = ( gCharge < 0)? -1 : 1 ;
764 if(gCharge > 0) nPlusChargeTruth += 1.;
765 if(gCharge < 0) nMinusChargeTruth += 1.;
766
767 if(fUsePid){
768
769 Double_t rap = partMC->Y();
770 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
771 if(TMath::Abs(partMC->GetPdgCode()) != gPdgCode) continue;
772
773 if(gCharge > 0) nParticleTruth += 1.;
774 if(gCharge < 0) nAntiParticleTruth += 1.;
775
776 }//if(fUsePid) ----
777
778 if( fCheckEff ){
779
780 Int_t chrgRecStatus = 0;
781 Int_t pidRecStatus = 0;
782
783 Double_t ptRec = 0.;
784 Double_t etaRec = 0.;
785 Double_t phiRec = 0.;
786
787 for( Int_t iRec = 0; iRec < nTracks; iRec++ ){
788
789 if( iMC == fLabel[0][iRec] ){
790
791 AliAODTrack* aodTrack = NULL;
792 if(fAOD)
793 aodTrack = fAOD->GetTrack(iRec);
794 if(aodTrack){
795
796 chrgRecStatus = 1;
797
798 if(fUsePid){
799
800 if( iMC == fLabel[1][iRec] ){
801 pidRecStatus = 1;
802 ptRec = aodTrack->Pt();
803 etaRec = aodTrack->Eta();
804 phiRec = aodTrack->Phi();
805 }
806 }//fUsePid--
807
808 }//aodTrack
809
810 break;
811
812 }
813
814 }//Reconstd track--
815
816 Double_t effContainer[10] = { fCentrality, chargeState, chrgRecStatus, pidRecStatus, partMC->Pt(), partMC->Eta(), partMC->Phi(), ptRec, etaRec, phiRec };
817
818 fTHnEfficiencyHisto->Fill(effContainer);
819 }
820
821 }//MC-Truth Track loop--
822
823 Double_t fContainerChTruth[3] = { fCentrality, nPlusChargeTruth, nMinusChargeTruth };
824 Double_t fContainerPidTruth[3] = { fCentrality, nParticleTruth, nAntiParticleTruth };
825
826 fTHnCentNplusNminusChTruth->Fill(fContainerChTruth);//MC -Truth ch. particles
827
828 if( fUsePid ){
829 gPid = (Int_t)fParticleSpecies;
830 fTHnCentNplusNminusPidTruth[gPid]->Fill(fContainerPidTruth);//MC-Truth pid
831 }
832
833 fEventCounter->Fill(7);
834
835 return;
8768ec6a 836
837}
838
e405901b 839//---------------------------------------------------------------------------------------
840Bool_t AliEbyEHigherMomentsTask::ProperVertex(){
841
842 Bool_t IsVtx = kFALSE;
843
844 const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
845
846 if(vertex) {
847 Double32_t fCov[6];
848 vertex->GetCovarianceMatrix(fCov);
849 if(vertex->GetNContributors() > 0) {
850 if(fCov[5] != 0) {
851
852 Double_t lvx = vertex->GetX();
853 Double_t lvy = vertex->GetY();
854 Double_t lvz = vertex->GetZ();
855
856 fEventCounter->Fill(5);
857
858 fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz);
859
860 if(TMath::Abs(lvx) < fVxMax) {
861 if(TMath::Abs(lvy) < fVyMax) {
862 if(TMath::Abs(lvz) < fVzMax) {
863
864 fEventCounter->Fill(6);
865 fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz);
866
867 IsVtx = kTRUE;
868
869 }//Z-Vertex cut---
870 }//Y-vertex cut--
871 }//X-vertex cut---
872 }//Covariance------
873 }//Contributors check---
874 }//If vertex-----------
875
876 return IsVtx;
877}
878//---------------------------------------------------------------------------------------
8768ec6a 879
880//------------------------------------------------------------------------
881void AliEbyEHigherMomentsTask::Terminate( Option_t * ){
882
883 Info("AliEbyEHigherMomentTask"," Task Successfully finished");
884
885}
886
887//------------------------------------------------------------------------