]> 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
new analysis task: multiplicity fluctuations (Maitreyee Mukherjee <maitreyee.mukherje...
[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 ),
73 fListOfHistos(0),
e405901b 74 fAOD(0),
75 fArrayMC(0),
76 fPIDResponse(0),
77 fParticleSpecies(AliPID::kProton),
beb13b6c 78 fAnalysisType(0),
8768ec6a 79 fCentralityEstimator("V0M"),
e405901b 80 fCentrality(0),
8768ec6a 81 fVxMax(3.),
82 fVyMax(3.),
e405901b 83 fVzMax(10.),
8768ec6a 84 fDCAxy(2.4),
85 fDCAz(3.),
86 fPtLowerLimit(0.3),
87 fPtHigherLimit(1.5),
88 fEtaLowerLimit(-0.8),
89 fEtaHigherLimit(0.8),
e405901b 90 fRapidityCut(0.5),
91 fNSigmaCut(3.),
8768ec6a 92 fTPCNClus(80),
abe73236 93 fChi2perNDF(4.),
e405901b 94 fAODtrackCutBit(128),
e405901b 95 fUsePid(kFALSE),
e405901b 96 fEventCounter(0),
97 fHistDCA(0),
98 fTPCSig(0),
99 fTPCSigA(0),
100 fTHnCentNplusNminusCh(0),
101 fTHnCentNplusNminusChTruth(0),
bfa00cce 102 fTHnCentNplusNminus(0)
8768ec6a 103{
104
abe73236 105 for ( Int_t i = 0; i < 13; i++) {
8768ec6a 106 fHistQA[i] = NULL;
107 }
abe73236 108
e405901b 109 for ( Int_t i = 0; i < 5; i++ ){
110 fTHnCentNplusNminusPid[i] = NULL;
111 fTHnCentNplusNminusPidTruth[i] = NULL;
112 }
113
8768ec6a 114 DefineOutput(1, TList::Class()); // Outpput....
bfa00cce 115 //DefineOutput(2, TList::Class());
8768ec6a 116}
117
118AliEbyEHigherMomentsTask::~AliEbyEHigherMomentsTask() {
bfa00cce 119 //if(fListOfHistosQA) delete fListOfHistosQA;
8768ec6a 120 if(fListOfHistos) delete fListOfHistos;
bfa00cce 121
8768ec6a 122}
123
124//---------------------------------------------------------------------------------
125void AliEbyEHigherMomentsTask::UserCreateOutputObjects() {
e405901b 126
bfa00cce 127 // fListOfHistosQA = new TList();
128 //fListOfHistosQA->SetOwner(kTRUE);
8768ec6a 129 fListOfHistos = new TList();
130 fListOfHistos->SetOwner(kTRUE);
8768ec6a 131
e405901b 132 fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
133 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
134 fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-90% centrality");
135 fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
136 fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
137 fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
138 fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
bfa00cce 139 fListOfHistos->Add(fEventCounter);
e405901b 140
141 //For QA-Histograms
8768ec6a 142 fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.);
143 fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.);
144 fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
145 fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.);
146 fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.);
147 fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.);
148 fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.);
149 fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.);
150 fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
151 fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
152 fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
abe73236 153 fHistQA[11] = new TH1D("fHistQANCls","Number of TPC cluster",200,0,200);
154 fHistQA[12] = new TH1D("fHistQAChi2","Chi2 per NDF",100,0,10);
abe73236 155 for(Int_t i = 0; i < 13; i++)
8768ec6a 156 {
bfa00cce 157 fListOfHistos->Add(fHistQA[i]);
8768ec6a 158 }
abe73236 159
e405901b 160 fHistDCA = new TH2D("fHistDCA","DCAxy Vs DCAz", 500, -5., 5., 500, -5., 5.);
161 fTPCSig = new TH2D("fTPCSig","TPC signal",200, 0.0, 10. ,1000,0.,1000);
162 fTPCSig->SetMarkerColor(kRed);
163 fTPCSigA = new TH2D("fTPCSigA","TPC signal all ",200, 0.0, 10. ,1000,0.,1000);
bfa00cce 164 fListOfHistos->Add(fHistDCA);
165 fListOfHistos->Add(fTPCSig);
166 fListOfHistos->Add(fTPCSigA);
e405901b 167
168 const Int_t nDim = 3;
169 const Int_t nPid = 5;
170 TString Species[nPid] = {"Electron","Muon","Pion","Kaon","Proton"};
8768ec6a 171
e405901b 172 Int_t fBins[nPid][nDim];
173 Double_t fMin[nPid][nDim];
174 Double_t fMax[nPid][nDim];
8768ec6a 175
e405901b 176 for( Int_t i = 0; i < nPid; i++ ){
177 for( Int_t j = 0; j < nDim; j++ ){
178 fBins[i][j] = 0;
179 fMin[i][j] = 0.;
180 fMax[i][j] = 0.;
181 }
182 }
8768ec6a 183
abe73236 184
e405901b 185 Int_t fBinsCh[nDim] = {100, 1500, 1500};
186 Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
187 Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
188
189 fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
190 fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
191 fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
192 fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
193 fListOfHistos->Add(fTHnCentNplusNminusCh);
194
195 if( fAnalysisType == "MCAOD" ){
196
197 fTHnCentNplusNminusChTruth = new THnSparseD("fTHnCentNplusNminusChTruth","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
198 fTHnCentNplusNminusChTruth->GetAxis(0)->SetTitle("Centrality");
199 fTHnCentNplusNminusChTruth->GetAxis(1)->SetTitle("Nplus");
200 fTHnCentNplusNminusChTruth->GetAxis(2)->SetTitle("Nminus");
201 fListOfHistos->Add(fTHnCentNplusNminusChTruth);
202
bfa00cce 203
e405901b 204 }//MCAOD---
205
206 TString hname1, hname11;
207 TString htitle1, axisTitle1, axisTitle2;
208
209
210 if( fUsePid ){
211
212 Int_t gPid = (Int_t)fParticleSpecies;
213
214 if( gPid > 1 && gPid < 5 ){
215 //Pion----
216 fBins[2][0] = 100; fBins[2][1] = 1000; fBins[2][2] = 1000;
217 fMin[2][0] = -0.5; fMin[2][1] = -0.5; fMin[2][2] = -0.5;
218 fMax[2][0] = 99.5; fMax[2][1] = 999.5; fMax[2][2] = 999.5;
219 //Kaon------
220 fBins[3][0] = 100; fBins[3][1] = 600; fBins[3][2] = 600;
221 fMin[3][0] = -0.5; fMin[3][1] = -0.5; fMin[3][2] = -0.5;
222 fMax[3][0] = 99.5; fMax[3][1] = 599.5; fMax[3][2] = 599.5;
223 //Proton-----
224 fBins[4][0] = 100; fBins[4][1] = 400; fBins[4][2] = 400;
225 fMin[4][0] = -0.5; fMin[4][1] = -0.5; fMin[4][2] = -0.5;
226 fMax[4][0] = 99.5; fMax[4][1] = 399.5; fMax[4][2] = 399.5;
227
228 hname1 = "fCentNplusNminusPid"; hname1 +=gPid;
229 htitle1 = Species[gPid]; htitle1 +=" And Neg-"; htitle1 +=Species[gPid];
230 axisTitle1 ="Pos-"; axisTitle1 += Species[gPid];
231 axisTitle2 ="Neg-"; axisTitle2 += Species[gPid];
232
233 fTHnCentNplusNminusPid[gPid] = new THnSparseD(hname1.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
234
235 fTHnCentNplusNminusPid[gPid]->GetAxis(0)->SetTitle("Centrality");
236 fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
237 fTHnCentNplusNminusPid[gPid]->GetAxis(2)->SetTitle(axisTitle2.Data());
238
239 fListOfHistos->Add(fTHnCentNplusNminusPid[gPid]);
240
241 if( fAnalysisType == "MCAOD" ){
242
243 hname11 = "fCentNplusNminusPidTruth"; hname11 +=gPid;
244 fTHnCentNplusNminusPidTruth[gPid] = new THnSparseD(hname11.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
245
246 fTHnCentNplusNminusPidTruth[gPid]->GetAxis(0)->SetTitle("Centrality");
247 fTHnCentNplusNminusPidTruth[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
248 fTHnCentNplusNminusPidTruth[gPid]->GetAxis(2)->SetTitle(axisTitle2.Data());
249 fListOfHistos->Add(fTHnCentNplusNminusPidTruth[gPid]);
250
251 }//MCAOD-----
252 }//Pion, Koan and Proton-------
253 else{
254
255 Int_t fBinsX[nDim] = {100, 1500, 1500};
256 Double_t fMinX[nDim] = { -0.5, -0.5, -0.5 };
257 Double_t fMaxX[nDim] = { 99.5, 1499.5, 1499.5};
258
259 fTHnCentNplusNminus = new THnSparseD("fTHnCentNplusNminus","Cent-NplusChrg-NminusChrg", nDim, fBinsX, fMinX, fMaxX);
260 fTHnCentNplusNminus->GetAxis(0)->SetTitle("Centrality");
261 fTHnCentNplusNminus->GetAxis(1)->SetTitle("UnwantedPlus");
262 fTHnCentNplusNminus->GetAxis(2)->SetTitle("UnwantedMinus");
263 fListOfHistos->Add(fTHnCentNplusNminus);
264
265 }//Unwanted particle -(other than pi, K or proton)
266
267 }//fUsePid-------
268
bfa00cce 269 //PostData(1, fListOfHistosQA);
270 PostData(1, fListOfHistos);
e405901b 271
abe73236 272
8768ec6a 273}
274
275
276//----------------------------------------------------------------------------------
277void AliEbyEHigherMomentsTask::UserExec( Option_t * ){
abe73236 278
8768ec6a 279 fEventCounter->Fill(1);
e405901b 280
281 if(fAnalysisType == "AOD") {
282
283 doAODEvent();
284
285 }//AOD--analysis-----
286
287 else if(fAnalysisType == "MCAOD") {
288
289 doMCAODEvent();
290
291 }
292
293 else return;
294
295 fEventCounter->Fill(8);
296
e405901b 297
298}
299
300//--------------------------------------------------------------------------------------
301void AliEbyEHigherMomentsTask::doAODEvent(){
302
bfa00cce 303 Double_t positiveSum = 0.;
304 Double_t negativeSum = 0.;
305 Double_t posPidSum = 0.;
306 Double_t negPidSum = 0.;
e405901b 307 Int_t gPid = 0;
8768ec6a 308
e405901b 309 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
310 if (!manager) {
311 cout<<"ERROR: Analysis manager not found."<<endl;
8768ec6a 312 return;
313 }
e405901b 314 //coneect to the inputHandler------------
315 AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
316 if (!inputHandler) {
317 cout<<"ERROR: Input handler not found."<<endl;
318 return;
319 }
320
321 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
322 if (!fAOD)
323 {
324 cout<< "ERROR 01: AOD not found " <<endl;
325 return;
326 }
8768ec6a 327
e405901b 328 fPIDResponse =inputHandler->GetPIDResponse();
8768ec6a 329
e405901b 330
331 if (!fPIDResponse){
332 AliFatal("This Task needs the PID response attached to the inputHandler");
333 return;
334 }
335
336 AliAODHeader *aodHeader = fAOD->GetHeader();
337
338 fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
339
340 if(fCentrality < 0 || fCentrality >= 91) return;
341
342
343
344 fEventCounter->Fill(2);
345
346 if(!ProperVertex()) return;
347
348 Int_t nTracks = fAOD->GetNumberOfTracks();
349
350 TExMap *trackMap = new TExMap();//Mapping matrix----
351
352 for(Int_t i = 0; i < nTracks; i++) {
8768ec6a 353
e405901b 354 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
abe73236 355
e405901b 356 if(!aodTrack) {
357 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
358 continue;
359 }
abe73236 360
e405901b 361 Double_t tpcSignalAll = aodTrack->GetTPCsignal();
362 fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
abe73236 363
e405901b 364 Int_t gID = aodTrack->GetID();
abe73236 365
e405901b 366 if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
367 }//1st track loop----
368
369 AliAODTrack* newAodTrack;
370
371 for( Int_t j = 0; j < nTracks; j++ )
372 {
373
374 AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
375
376 if(!aodTrack1) {
377 AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
378 continue;
379 }
380
381
382 if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
383
384 Int_t gID = aodTrack1->GetID();
385
386 //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
387 newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
388
389
390 Double_t pt = aodTrack1->Pt();
391 Double_t eta = aodTrack1->Eta();
392 Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
393 Double_t chi2ndf = aodTrack1->Chi2perNDF();
394
395 if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
396 if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
397 if( nclus < fTPCNClus ) continue;
398 if( chi2ndf > fChi2perNDF ) continue;
399
400 if( fAODtrackCutBit == 128 ){
401 //TPC only tracks----
402 Float_t dxy = 0., dz = 0.;
403 dxy = aodTrack1->DCA();
404 dz = aodTrack1->ZAtDCA();
405
406 if( fabs(dxy) > fDCAxy ) continue;
407 if( fabs(dz) > fDCAz ) continue;
408 //Extra cut on DCA---( Similar to BF Task.. )
409 if( fDCAxy !=0 && fDCAz !=0 ){
410 if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
411 }
412
413 fHistQA[6]->Fill(dxy);
414 fHistQA[7]->Fill(dz);
415 fHistDCA->Fill(dxy,dz);
416
417 }
418
419 fHistQA[8]->Fill(pt);
420 fHistQA[9]->Fill(eta);
421 fHistQA[10]->Fill(aodTrack1->Phi());
422 fHistQA[11]->Fill(nclus);
423 fHistQA[12]->Fill(chi2ndf);
424
425 Short_t gCharge = aodTrack1->Charge();
426
bfa00cce 427 if(gCharge > 0) positiveSum += 1.;
428 if(gCharge < 0) negativeSum += 1.;
e405901b 429
430 if( fUsePid ) {
431
432 gPid = (Int_t)fParticleSpecies;
433
434 Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
435 Double_t tpcSignal = newAodTrack->GetTPCsignal();
436 //Double_t rap = aodTrack1->Y(AliPID::ParticleMass(fParticleSpecies));
437 //Double_t tpcSignal = aodTrack1->GetTPCsignal();
438
439 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
440
441 fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
442
443 Float_t nsigmaTPCPID = -999.;
444 Float_t nsigmaTOFPID = -999.;
445 //Float_t nsigmaTPCTOFPID = -999.;
446
447 nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
448 nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
449
450 if ( nsigmaTPCPID < fNSigmaCut ){
abe73236 451
bfa00cce 452 if (gCharge > 0) posPidSum +=1.;
453 if( gCharge < 0 ) negPidSum +=1.;
abe73236 454
e405901b 455 }
456 }//fUsepid-----
457
458 }//--------- Track Loop to select with filterbit
8768ec6a 459
e405901b 460
461 //cout << fCentrality <<" "<< nPlusCharge << " " << nMinusCharge << endl;
462 //cout << fCentrality <<" "<< nParticle << " " << nAntiParticle << endl;
463
bfa00cce 464 Double_t fContainerCh[3] = { fCentrality, positiveSum, negativeSum};
465
e405901b 466
467
468 fTHnCentNplusNminusCh->Fill(fContainerCh);
469
470 if( fUsePid ){
471 gPid = (Int_t)fParticleSpecies;
bfa00cce 472 Double_t fContainerPid[3] = { fCentrality, posPidSum, negPidSum};
e405901b 473 fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
474 }
475
476 fEventCounter->Fill(7);
bfa00cce 477 PostData(1, fListOfHistos);
e405901b 478 return;
479
480}
481//--------------------------------------------------------------------------------------
482//--------------------------------------------------------------------------------------
483void AliEbyEHigherMomentsTask::doMCAODEvent(){
484
485
486 //---------
bfa00cce 487 Double_t positiveSumMCRec = 0.;
488 Double_t negativeSumMCRec = 0.;
489 Double_t posPidSumMCRec = 0.;
490 Double_t negPidSumMCRec = 0.;
491
492 Double_t positiveSumMCTruth = 0.;
493 Double_t negativeSumMCTruth = 0.;
494 Double_t posPidSumMCTruth = 0.;
495 Double_t negPidSumMCTruth = 0.;
e405901b 496
497 Int_t gPid = 0;
498 Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
499
500 //Connect to Anlaysis manager------
501 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
502 if (!manager) {
503 cout<<"ERROR: Analysis manager not found."<<endl;
504 return;
505 }
506
507 AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
508 if (!inputHandler) {
509 cout<<"ERROR: Input handler not found."<<endl;
510 return;
511 }
512
513 //AOD event------
514 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
515 if (!fAOD)
516 {
517 cout<< "ERROR 01: AOD not found " <<endl;
beb13b6c 518 return;
8768ec6a 519 }
e405901b 520
521 fPIDResponse =inputHandler->GetPIDResponse();
522
523
524 if (!fPIDResponse){
525 AliFatal("This Task needs the PID response attached to the inputHandler");
526 return;
527 }
528
529 // -- Get MC header
530 // ------------------------------------------------------------------
531
532 fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
533 if (!fArrayMC) {
534 AliFatal("Error: MC particles branch not found!\n");
535 return;
536 }
537
538 AliAODMCHeader *mcHdr=NULL;
539 mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
540 if(!mcHdr) {
541 printf("MC header branch not found!\n");
542 return;
543 }
544
545 AliAODHeader *aodHeader = fAOD->GetHeader();
546
547 fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
548
549
550
551 if( fCentrality < 0 || fCentrality >= 91) return;
552
553 fEventCounter->Fill(2);
554
555 if(!ProperVertex()) return;
556
557 Int_t nTracks = fAOD->GetNumberOfTracks();
558
e405901b 559 TExMap *trackMap = new TExMap();//Mapping matrix----
560
561 for(Int_t i = 0; i < nTracks; i++) {
8768ec6a 562
e405901b 563 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
8768ec6a 564
e405901b 565 if(!aodTrack) {
566 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
567 continue;
568 }
8768ec6a 569
e405901b 570 Double_t tpcSignalAll = aodTrack->GetTPCsignal();
571 fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
572
573 Int_t gID = aodTrack->GetID();
574
575 if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global(primary) tracks-----
576
577 }//1st track loop----
578
579 AliAODTrack* newAodTrack;
580
581 for( Int_t j = 0; j < nTracks; j++ ){
582
583 AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
584
585 if(!aodTrack1) {
586 AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
587 continue;
abe73236 588 }
8768ec6a 589
e405901b 590
591 if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
592
593 Int_t gID = aodTrack1->GetID();
594
595 //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
596
597 newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
598
599
600 //cout << dxy << endl;
601 Double_t pt = aodTrack1->Pt();
602 Double_t eta = aodTrack1->Eta();
603 Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
604 Double_t chi2ndf = aodTrack1->Chi2perNDF();
605
606 if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
607 if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
608 if( nclus < fTPCNClus ) continue;
609 if( chi2ndf > fChi2perNDF ) continue;
610
611 if( fAODtrackCutBit == 128 ){
612 Float_t dxy = 0., dz = 0.;
613 dxy = aodTrack1->DCA();
614 dz = aodTrack1->ZAtDCA();
615
616 if( fabs(dxy) > fDCAxy ) continue;
617 if( fabs(dz) > fDCAz ) continue;
618 //Extra cut on DCA---( Similar to BF Task.. )
619 if( fDCAxy !=0 && fDCAz !=0 ){
620 if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
beb13b6c 621 }
8768ec6a 622
e405901b 623 fHistQA[6]->Fill(dxy);
624 fHistQA[7]->Fill(dz);
625 fHistDCA->Fill(dxy,dz);
abe73236 626
e405901b 627 }
628
629
630
631 fHistQA[8]->Fill(pt);
632 fHistQA[9]->Fill(eta);
633 fHistQA[10]->Fill(aodTrack1->Phi());
634 fHistQA[11]->Fill(nclus);
635 fHistQA[12]->Fill(chi2ndf);
636
637
638 Short_t gCharge = aodTrack1->Charge();
639
640 if( gCharge == 0 ) continue;
bfa00cce 641
642 if(gCharge > 0) positiveSumMCRec += 1.;
643 if(gCharge < 0) negativeSumMCRec += 1.;
e405901b 644
645 if( fUsePid ) {
646
647 gPid = (Int_t)fParticleSpecies;
648 Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
649 Double_t tpcSignal = newAodTrack->GetTPCsignal();
650
651 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
652
653 fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
654
655 Float_t nsigmaTPCPID = -999.;
656 Float_t nsigmaTOFPID = -999.;
657 //Float_t nsigmaTPCTOFPID = -999.;
658
659 nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
660 nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
661
662 if( nsigmaTPCPID < fNSigmaCut ){
bfa00cce 663
664 if (gCharge > 0) posPidSumMCRec +=1;
665 if( gCharge < 0 ) negPidSumMCRec +=1.;
e405901b 666 }//nSigmaCut-----
667 }//fUsepid-----
668
669 }//--------- Track Loop to select with filterbit
abe73236 670
e405901b 671
bfa00cce 672 Double_t fContainerCh[3] = { fCentrality, positiveSumMCRec, negativeSumMCRec};//Reco. values ch. hadrons
bed5186b 673 fTHnCentNplusNminusCh->Fill(fContainerCh);//Fill the rec. ch. particles---
bfa00cce 674
e405901b 675 if( fUsePid ){
676 gPid = (Int_t)fParticleSpecies;
bfa00cce 677 Double_t fContainerPid[3] = { fCentrality, posPidSumMCRec, negPidSumMCRec};//Reco. values pid.
bed5186b 678 fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);//Fill the rec. pid tracks
e405901b 679 }
680
681
682 //cout << fCentrality << " "<< nPlusCharge << " "<< nMinusCharge << endl;
683 //cout << nParticle <<" And " << nAntiParticle << endl;
684 //===========================================
685 //--------MC Truth---------------------------
686 //===========================================
687
688 Int_t nMCTrack = fArrayMC->GetEntriesFast();
689
690 for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
beb13b6c 691
e405901b 692 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
693
694 if(!partMC){
695 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
696 continue;
697 }
698
699 if(partMC->Charge() == 0) continue;
700 if(!partMC->IsPhysicalPrimary()) continue;
701
702 if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
703 if (partMC->Pt() < fPtLowerLimit || partMC->Pt() > fPtHigherLimit) continue;
704
705 Short_t gCharge = partMC->Charge();
0660fa85 706
bfa00cce 707 if(gCharge > 0) positiveSumMCTruth += 1.;
708 if(gCharge < 0) negativeSumMCTruth += 1.;
e405901b 709
710 if(fUsePid){
711
712 Double_t rap = partMC->Y();
713 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
714 if(TMath::Abs(partMC->GetPdgCode()) != gPdgCode) continue;
715
bfa00cce 716 if(gCharge > 0) posPidSumMCTruth += 1.;
717 if(gCharge < 0) negPidSumMCTruth += 1.;
e405901b 718
719 }//if(fUsePid) ----
720
bfa00cce 721
e405901b 722 }//MC-Truth Track loop--
723
bfa00cce 724 Double_t fContainerChTruth[3] = { fCentrality, positiveSumMCTruth, negativeSumMCTruth };
e405901b 725 fTHnCentNplusNminusChTruth->Fill(fContainerChTruth);//MC -Truth ch. particles
726
727 if( fUsePid ){
728 gPid = (Int_t)fParticleSpecies;
bfa00cce 729 Double_t fContainerPidTruth[3] = { fCentrality, posPidSumMCTruth, negPidSumMCTruth};
e405901b 730 fTHnCentNplusNminusPidTruth[gPid]->Fill(fContainerPidTruth);//MC-Truth pid
731 }
732
733 fEventCounter->Fill(7);
bfa00cce 734 PostData(1, fListOfHistos);
e405901b 735 return;
8768ec6a 736
737}
738
e405901b 739//---------------------------------------------------------------------------------------
740Bool_t AliEbyEHigherMomentsTask::ProperVertex(){
741
742 Bool_t IsVtx = kFALSE;
743
744 const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
745
746 if(vertex) {
747 Double32_t fCov[6];
748 vertex->GetCovarianceMatrix(fCov);
749 if(vertex->GetNContributors() > 0) {
750 if(fCov[5] != 0) {
751
752 Double_t lvx = vertex->GetX();
753 Double_t lvy = vertex->GetY();
754 Double_t lvz = vertex->GetZ();
755
756 fEventCounter->Fill(5);
757
758 fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz);
759
760 if(TMath::Abs(lvx) < fVxMax) {
761 if(TMath::Abs(lvy) < fVyMax) {
762 if(TMath::Abs(lvz) < fVzMax) {
763
764 fEventCounter->Fill(6);
765 fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz);
766
767 IsVtx = kTRUE;
768
769 }//Z-Vertex cut---
770 }//Y-vertex cut--
771 }//X-vertex cut---
772 }//Covariance------
773 }//Contributors check---
774 }//If vertex-----------
775
776 return IsVtx;
777}
778//---------------------------------------------------------------------------------------
8768ec6a 779
780//------------------------------------------------------------------------
781void AliEbyEHigherMomentsTask::Terminate( Option_t * ){
782
783 Info("AliEbyEHigherMomentTask"," Task Successfully finished");
784
785}
786
787//------------------------------------------------------------------------