1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: Satyajit Jena. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //=========================================================================//
19 // Analysis Task for Net-Charge Higher Moment Analysis //
20 // Author: Satyajit Jena || Nirbhay K. Behera //
21 // sjena@cern.ch || nbehera@cern.ch //
23 //=========================================================================//
31 #include "THnSparse.h"
33 #include "TParticle.h"
34 #include <TDatabasePDG.h>
35 #include <TParticlePDG.h>
37 #include "AliAnalysisTaskSE.h"
38 #include "AliAnalysisManager.h"
40 #include "AliESDVertex.h"
41 #include "AliESDEvent.h"
42 #include "AliESDInputHandler.h"
43 #include "AliAODEvent.h"
44 #include "AliAODTrack.h"
45 #include "AliAODInputHandler.h"
47 #include "AliESDEvent.h"
48 #include "AliAODEvent.h"
50 #include "AliESDtrackCuts.h"
51 #include "AliAODMCHeader.h"
52 #include "AliAODMCParticle.h"
53 #include "AliMCEventHandler.h"
54 #include "AliMCEvent.h"
55 #include "AliMCParticle.h"
57 #include "AliGenHijingEventHeader.h"
58 #include "AliGenEventHeader.h"
60 #include "AliAODPid.h"
61 #include "AliPIDResponse.h"
62 #include "AliAODpidUtil.h"
63 #include "AliPIDCombined.h"
65 #include "AliEbyEHigherMomentsEffContTask.h"
71 ClassImp(AliEbyEHigherMomentsEffContTask)
74 //-----------------------------------------------------------------------
75 AliEbyEHigherMomentsEffContTask::AliEbyEHigherMomentsEffContTask( const char *name )
76 : AliAnalysisTaskSE( name ),
82 fParticleSpecies(AliPID::kProton),
84 fCentralityEstimator("V0M"),
99 fAODtrackCutBit(1024),
108 fTHnCentNplusNminusCh(0),
109 fTHnCentNplusNminus(0),
114 for ( Int_t i = 0; i < 13; i++) {
118 for ( Int_t i = 0; i < 5; i++ ){
119 fTHnCentNplusNminusPid[i] = NULL;
122 DefineOutput(1, TList::Class()); // Outpput....
123 DefineOutput(2, TList::Class());
127 AliEbyEHigherMomentsEffContTask::~AliEbyEHigherMomentsEffContTask() {
129 if(fListOfHistosQA) delete fListOfHistosQA;
130 if(fListOfHistos) delete fListOfHistos;
131 if(fLabel[0] ) delete [] (fLabel[0]);
132 if(fLabel[1] ) delete [] (fLabel[1]);
136 //---------------------------------------------------------------------------------
137 void AliEbyEHigherMomentsEffContTask::UserCreateOutputObjects() {
139 fListOfHistosQA = new TList();
140 fListOfHistosQA->SetOwner(kTRUE);
142 fListOfHistos = new TList();
143 fListOfHistos->SetOwner(kTRUE);
145 fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
146 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
147 fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-90% centrality");
148 fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
149 fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
150 fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
151 fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
152 fListOfHistosQA->Add(fEventCounter);
154 fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.);
155 fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.);
156 fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
157 fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.);
158 fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.);
159 fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.);
160 fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.);
161 fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.);
162 fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
163 fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
164 fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
165 fHistQA[11] = new TH1D("fHistQANCls","Number of TPC cluster",200,0,200);
166 fHistQA[12] = new TH1D("fHistQAChi2","Chi2 per NDF",100,0,10);
168 for(Int_t i = 0; i < 13; i++)
170 fListOfHistosQA->Add(fHistQA[i]);
173 fHistDCA = new TH2D("fHistDCA","DCAxy Vs DCAz", 500, -5., 5., 500, -5., 5.);
174 fTPCSig = new TH2D("fTPCSig","TPC signal",200, 0.0, 10. ,1000,0.,1000);
175 fTPCSig->SetMarkerColor(kRed);
176 fTPCSigA = new TH2D("fTPCSigA","TPC signal all ",200, 0.0, 10. ,1000,0.,1000);
177 fListOfHistosQA->Add(fHistDCA);
178 fListOfHistosQA->Add(fTPCSig);
179 fListOfHistosQA->Add(fTPCSigA);
181 const Int_t nDim = 3;
182 const Int_t nPid = 5;
184 TString Species[nPid] = {"Electron","Muon","Pion","Kaon","Proton"};
186 Int_t fBins[nPid][nDim];
187 Double_t fMin[nPid][nDim];
188 Double_t fMax[nPid][nDim];
190 for( Int_t i = 0; i < nPid; i++ ){
191 for( Int_t j = 0; j < nDim; j++ ){
199 // if(fAnalysisType == "AOD" || fAnalysisType == "MCAOD")
203 TString htitle1, axisTitle1, axisTitle2;
209 Int_t gPid = (Int_t)fParticleSpecies;
211 if( gPid > 1 && gPid < 5 ){
213 fBins[2][0] = 100; fBins[2][1] = 1000; fBins[2][2] = 1000;
214 fMin[2][0] = -0.5; fMin[2][1] = -0.5; fMin[2][2] = -0.5;
215 fMax[2][0] = 99.5; fMax[2][1] = 999.5; fMax[2][2] = 999.5;
217 fBins[3][0] = 100; fBins[3][1] = 600; fBins[3][2] = 600;
218 fMin[3][0] = -0.5; fMin[3][1] = -0.5; fMin[3][2] = -0.5;
219 fMax[3][0] = 99.5; fMax[3][1] = 599.5; fMax[3][2] = 599.5;
221 fBins[4][0] = 100; fBins[4][1] = 400; fBins[4][2] = 400;
222 fMin[4][0] = -0.5; fMin[4][1] = -0.5; fMin[4][2] = -0.5;
223 fMax[4][0] = 99.5; fMax[4][1] = 399.5; fMax[4][2] = 399.5;
225 hname1 = "fCentNplusNminusPid"; hname1 +=gPid;
226 htitle1 = Species[gPid]; htitle1 +=" And Neg-"; htitle1 +=Species[gPid];
227 axisTitle1 = Species[gPid];
228 axisTitle2 ="Neg-"; axisTitle2 += Species[gPid];
230 fTHnCentNplusNminusPid[gPid] = new THnSparseD(hname1.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
232 fTHnCentNplusNminusPid[gPid]->GetAxis(0)->SetTitle("Centrality");
233 fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
234 fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle2.Data());
236 fListOfHistos->Add(fTHnCentNplusNminusPid[gPid]);
238 }//Pion, Koan and Proton-------
242 Int_t fBinsX[nDim] = {100, 1500, 1500};
243 Double_t fMinX[nDim] = { -0.5, -0.5, -0.5 };
244 Double_t fMaxX[nDim] = { 99.5, 1499.5, 1499.5};
246 fTHnCentNplusNminus = new THnSparseD("fTHnCentNplusNminus","Cent-NplusChrg-NminusChrg", nDim, fBinsX, fMinX, fMaxX);
247 fTHnCentNplusNminus->GetAxis(0)->SetTitle("Centrality");
248 fTHnCentNplusNminus->GetAxis(1)->SetTitle("UnwantedPlus");
249 fTHnCentNplusNminus->GetAxis(2)->SetTitle("UnwantedMinus");
250 fListOfHistos->Add(fTHnCentNplusNminus);
252 }//Unwanted particle -------
258 Int_t fBinsCh[nDim] = {100, 1500, 1500};
259 Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
260 Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
262 fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
263 fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
264 fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
265 fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
266 fListOfHistos->Add(fTHnCentNplusNminusCh);
267 }//All charge hadrons---------
271 if( fAnalysisType == "MCAOD" ){
273 Double_t dCent[2] = {-0.5, 9.5};
276 Double_t dEta[2] = {-0.9, 0.9};
277 Int_t iEta = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ;
279 Double_t dRap[2] = {-0.5, 0.5};
280 Int_t iRap = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ;
282 Double_t dPhi[2] = {0.0, TMath::TwoPi()};
285 Double_t dPt[2] = {0.1, 3.0};
286 Int_t iPt = Int_t((dPt[1]-dPt[0]) / 0.075);
288 Double_t dSign[2] = {-1.5, 1.5};
293 Int_t nBinEff[15] = { iCent, iEta, iRap, iPhi, iPt, iSign, 2, 2, 2, iEta, iPhi, iPt, iEta, 2*iPhi+1, 2*iPt+1};
294 Double_t nMinEff[15] = { dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0], -0.5, -0.5, -0.5, dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1] };
295 Double_t nMaxEff[15] = { dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1], 1.5, 1.5, 1.5, dEta[1], dPhi[1], dPt[1], dEta[1], dPhi[1], dPt[1] };
297 fTHnEff = new THnSparseF("fHnEff", "cent:etaMC:yMC:phiMC:ptMC:sign:findable:recStatus:pidStatus:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt", 15, nBinEff, nMinEff, nMaxEff);
301 fTHnEff->GetAxis(0)->SetTitle("centrality");
302 fTHnEff->GetAxis(1)->SetTitle("#eta_{MC}");
303 fTHnEff->GetAxis(2)->SetTitle("#it{y}_{MC}");
304 fTHnEff->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");
305 fTHnEff->GetAxis(4)->SetTitle("p_{T,MC} (GeV/#it{c})");
306 fTHnEff->GetAxis(5)->SetTitle("sign");
307 fTHnEff->GetAxis(6)->SetTitle("findable");
308 fTHnEff->GetAxis(7)->SetTitle("recStatus");
309 fTHnEff->GetAxis(8)->SetTitle("recPid");
310 fTHnEff->GetAxis(9)->SetTitle("#eta_{Rec}");
311 fTHnEff->GetAxis(10)->SetTitle("#varphi_{Rec} (rad)");
312 fTHnEff->GetAxis(11)->SetTitle("p_{T,Rec} (GeV/#it{c})");
313 fTHnEff->GetAxis(12)->SetTitle("#eta_{MC}-#eta_{Rec}");
314 fTHnEff->GetAxis(13)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");
315 fTHnEff->GetAxis(14)->SetTitle("p_{T,MC}-p_{T,Rec} (GeV/#it{c})");
317 fListOfHistos->Add(fTHnEff);
322 Int_t nBinCont[14] = { iCent, iEta, iRap, iPhi, iPt, iSign, 8, iSign, iEta, iPhi, iPt, iEta, 2*iPhi+1, 2*iPt+1 };
323 Double_t nMinCont[14] = {dCent[0], dEta[0], dRap[0], dPhi[0], dPt[0], dSign[0], 0.5, dSign[0], dEta[0], dPhi[0], dPt[0], dEta[0], -dPhi[1], -dPt[1]};
324 Double_t nMaxCont[14] = { dCent[1], dEta[1], dRap[1], dPhi[1], dPt[1], dSign[1], 8.5, dSign[1], dEta[1], dPhi[1], dPt[1], dEta[1], dPhi[1], dPt[1] };
326 fTHnCont = new THnSparseF("fHnCont", "cent:etaMC:yMC:phiMC:ptMC:sign:contPart:contSign:etaRec:phiRec:ptRec:deltaEta:deltaPhi:deltaPt",14, nBinCont, nMinCont, nMaxCont);
329 fTHnCont->GetAxis(0)->SetTitle("centrality");
330 fTHnCont->GetAxis(1)->SetTitle("#eta_{MC}");
331 fTHnCont->GetAxis(2)->SetTitle("#it{y}_{MC}");
332 fTHnCont->GetAxis(3)->SetTitle("#varphi_{MC} (rad)");
333 fTHnCont->GetAxis(4)->SetTitle("p_{T,MC} (GeV/#it{c})");
334 fTHnCont->GetAxis(5)->SetTitle("sign");
335 fTHnCont->GetAxis(6)->SetTitle("contPart");
336 fTHnCont->GetAxis(7)->SetTitle("contSign");
337 fTHnCont->GetAxis(8)->SetTitle("#eta_{Rec}");
338 fTHnCont->GetAxis(9)->SetTitle("#varphi_{Rec} (rad)");
339 fTHnCont->GetAxis(10)->SetTitle("p_{T,Rec} (GeV/#it{c})");
340 fTHnCont->GetAxis(11)->SetTitle("#eta_{MC}-#eta_{Rec}");
341 fTHnCont->GetAxis(12)->SetTitle("#varphi_{MC}-#varphi_{Rec} (rad)");
342 fTHnCont->GetAxis(13)->SetTitle("p_{T,MC}-p_{T,Rec} (GeV/#it{c})");
344 fListOfHistos->Add(fTHnCont);
349 PostData(1, fListOfHistosQA);
350 PostData(2, fListOfHistos);
355 //----------------------------------------------------------------------------------
356 void AliEbyEHigherMomentsEffContTask::UserExec( Option_t * ){
358 fEventCounter->Fill(1);
361 if(fAnalysisType == "AOD") {
365 }//AOD--analysis-----
367 else if(fAnalysisType == "MCAOD") {
378 fEventCounter->Fill(8);
380 PostData(1, fListOfHistosQA);
381 PostData(2, fListOfHistos);
385 //--------------------------------------------------------------------------------------
386 void AliEbyEHigherMomentsEffContTask::doAODEvent(){
388 //-------------------
389 Double_t nPlusCharge = 0.;
390 Double_t nMinusCharge = 0.;
391 Double_t nPartile = 0.;
392 Double_t nAntiParticle = 0.;
396 //connect to the analysis mannager-----
398 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
400 cout<<"ERROR: Analysis manager not found."<<endl;
403 //coneect to the inputHandler------------
404 AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
406 cout<<"ERROR: Input handler not found."<<endl;
410 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
413 cout<< "ERROR 01: AOD not found " <<endl;
417 fPIDResponse =inputHandler->GetPIDResponse();
421 AliFatal("This Task needs the PID response attached to the inputHandler");
426 AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
427 if(!aodHeader) AliFatal("Not a standard AOD");
430 cent = aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data());
433 fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data());
434 else if (cent == 10 || cent == -1.)
436 else if (cent > 0 && cent < 9)
437 fCentrality = cent + 1;
439 if(fCentrality < 0 || fCentrality >= 10) return;
441 fEventCounter->Fill(2);
443 if(!ProperVertex()) return;
445 Int_t nTracks = fAOD->GetNumberOfTracks();
447 TExMap *trackMap = new TExMap();//Mapping matrix----
449 for(Int_t i = 0; i < nTracks; i++) {
451 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
454 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
458 Double_t tpcSignalAll = aodTrack->GetTPCsignal();
459 fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
461 Int_t gID = aodTrack->GetID();
463 if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
464 }//1st track loop----
466 AliAODTrack* newAodTrack;
468 for( Int_t j = 0; j < nTracks; j++ )
471 AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
474 AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
479 if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
481 Int_t gID = aodTrack1->GetID();
483 //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
484 newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
486 Float_t dxy = 0., dz = 0.;
488 dxy = aodTrack1->DCA();
489 dz = aodTrack1->ZAtDCA();
491 Double_t pt = aodTrack1->Pt();
492 Double_t eta = aodTrack1->Eta();
493 Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
494 Double_t chi2ndf = aodTrack1->Chi2perNDF();
497 if( fabs(dxy) > fDCAxy ) continue;
498 if( fabs(dz) > fDCAz ) continue;
499 //Extra cut on DCA---( Similar to BF Task.. )
500 if( fDCAxy !=0 && fDCAz !=0 ){
501 if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
504 if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
505 if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
506 if( nclus < fTPCNClus ) continue;
507 if( chi2ndf > fChi2perNDF ) continue;
510 fHistQA[6]->Fill(dxy);
511 fHistQA[7]->Fill(dz);
512 fHistQA[8]->Fill(pt);
513 fHistQA[9]->Fill(eta);
514 fHistQA[10]->Fill(aodTrack1->Phi());
515 fHistQA[11]->Fill(nclus);
516 fHistQA[12]->Fill(chi2ndf);
517 fHistDCA->Fill(dxy,dz);
519 Short_t gCharge = aodTrack1->Charge();
521 if(gCharge > 0) nPlusCharge += 1.;
522 if(gCharge < 0) nMinusCharge += 1.;
526 gPid = (Int_t)fParticleSpecies;
528 Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
529 Double_t tpcSignal = newAodTrack->GetTPCsignal();
530 //Double_t rap = aodTrack1->Y(AliPID::ParticleMass(fParticleSpecies));
531 //Double_t tpcSignal = aodTrack1->GetTPCsignal();
533 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
535 fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
537 Float_t nsigmaTPCPID = -999.;
538 //Float_t nsigmaTOFPID = -999.;
539 //Float_t nsigmaTPCTOFPID = -999.;
541 nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
542 //nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
544 if ( nsigmaTPCPID < fNSigmaCut ){
546 if (gCharge > 0) nPartile +=1.;
547 if( gCharge < 0 ) nAntiParticle +=1.;
552 }//--------- Track Loop to select with filterbit
556 Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), nPlusCharge, nMinusCharge};
557 Double_t fContainerPid[3] = { static_cast<Double_t>(fCentrality), nPartile, nAntiParticle};
561 fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
563 // cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl;
567 fTHnCentNplusNminusCh->Fill(fContainerCh);
570 //cout << "nCentrality "<< fCentrality <<", nPlusCharge="<< nPlusCharge << ", nMinusCharge=" << nMinusCharge << endl;
572 fEventCounter->Fill(7);
577 //--------------------------------------------------------------------------------------
580 //--------------------------------------------------------------------------------------
581 void AliEbyEHigherMomentsEffContTask::doMCAODEvent(){
585 Double_t nPlusCharge = 0.;
586 Double_t nMinusCharge = 0.;
587 Double_t nPartile = 0.;
588 Double_t nAntiParticle = 0.;
591 //Connect to Anlaysis manager------
592 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
594 cout<<"ERROR: Analysis manager not found."<<endl;
598 AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
600 cout<<"ERROR: Input handler not found."<<endl;
605 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
608 cout<< "ERROR 01: AOD not found " <<endl;
612 fPIDResponse =inputHandler->GetPIDResponse();
616 AliFatal("This Task needs the PID response attached to the inputHandler");
620 // ------------------------------------------------------------------
622 fArrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
624 AliFatal("Error: MC particles branch not found!\n");
627 AliAODMCHeader *mcHdr=NULL;
628 mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
630 printf("MC header branch not found!\n");
634 AliAODHeader *aodHeader = dynamic_cast<AliAODHeader*>(fAOD->GetHeader());
635 if(!aodHeader) AliFatal("Not a standard AOD");
638 cent = aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data());
641 fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data());
642 else if (cent == 10 || cent == -1.)
644 else if (cent > 0 && cent < 9)
645 fCentrality = cent + 1;
647 if( fCentrality < 0 || fCentrality >= 10) return;
649 fEventCounter->Fill(2);
653 if(!ProperVertex()) return;
655 Int_t nTracks = fAOD->GetNumberOfTracks();
658 fLabel = new Int_t*[2];
659 fLabel[0] = new Int_t[nTracks]; //All charged hadrons----
660 fLabel[1] = new Int_t[nTracks]; //For Pid-----------
661 //Initialize labels----
664 AliError("Can't Get fLabel[0] ");
668 AliError("Can't Get fLabel[1] ");
672 for(Int_t i=0; i < 2; i++){
673 for(Int_t j=0; j < nTracks; j++){
679 TExMap *trackMap = new TExMap();//Mapping matrix----
681 for(Int_t i = 0; i < nTracks; i++) {
683 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
686 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
690 Double_t tpcSignalAll = aodTrack->GetTPCsignal();
691 fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
693 Int_t gID = aodTrack->GetID();
695 if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
697 }//1st track loop----
699 AliAODTrack* newAodTrack;
701 for( Int_t j = 0; j < nTracks; j++ ){
703 AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
706 AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
711 if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
714 Int_t gID = aodTrack1->GetID();
716 //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
718 newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
720 Float_t dxy = 0., dz = 0.;
722 dxy = aodTrack1->DCA();
723 dz = aodTrack1->ZAtDCA();
725 //cout << dxy << endl;
726 Double_t pt = aodTrack1->Pt();
727 Double_t eta = aodTrack1->Eta();
728 Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
729 Double_t chi2ndf = aodTrack1->Chi2perNDF();
732 if( fabs(dxy) > fDCAxy ) continue;
733 if( fabs(dz) > fDCAz ) continue;
734 //Extra cut on DCA---( Similar to BF Task.. )
735 if( fDCAxy !=0 && fDCAz !=0 ){
736 if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
741 if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
742 if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
743 // if( nclus < fTPCNClus ) continue;
744 //if( chi2ndf > fChi2perNDF ) continue;
747 fHistQA[6]->Fill(dxy);
748 fHistQA[7]->Fill(dz);
749 fHistQA[8]->Fill(pt);
750 fHistQA[9]->Fill(eta);
751 fHistQA[10]->Fill(aodTrack1->Phi());
752 fHistQA[11]->Fill(nclus);
753 fHistQA[12]->Fill(chi2ndf);
754 fHistDCA->Fill(dxy,dz);
756 Short_t gCharge = aodTrack1->Charge();
758 if( gCharge == 0 ) continue;
760 Int_t label = TMath::Abs(aodTrack1->GetLabel());
761 //fill the labels--------
762 fLabel[0][j] = label;
764 if(gCharge > 0) nPlusCharge += 1.;
765 if(gCharge < 0) nMinusCharge += 1.;
769 gPid = (Int_t)fParticleSpecies;;
770 Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
771 Double_t tpcSignal = newAodTrack->GetTPCsignal();
773 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
775 fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
777 Float_t nsigmaTPCPID = -999.;
778 //Float_t nsigmaTOFPID = -999.;
779 //Float_t nsigmaTPCTOFPID = -999.;
781 nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
782 //nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
784 if ( nsigmaTPCPID < fNSigmaCut ){
786 //Fill the labels----
787 fLabel[1][j] = label;
789 if (gCharge > 0) nPartile +=1.;
790 if( gCharge < 0 ) nAntiParticle +=1.;
794 //Check Contamination-------------------
797 CheckContTrackAOD(label, gCharge, j);
800 }//--------- Track Loop to select with filterbit
802 if(fEff) FillEffSparse();
804 //cout << "Cent=" << fCentrality << " MC-PlusChrg=" << nPlusCharge << " MC-MinusChrg=" << nMinusCharge << endl;
806 //cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl;
808 Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), nPlusCharge, nMinusCharge};
809 Double_t fContainerPid[3] = { static_cast<Double_t>(fCentrality), nPartile, nAntiParticle};
813 fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
818 fTHnCentNplusNminusCh->Fill(fContainerCh);
821 fEventCounter->Fill(7);
826 //---------------------------------------------------------------------------------------
828 //---------------------------------------------------------------------------------------
829 Bool_t AliEbyEHigherMomentsEffContTask::ProperVertex(){
831 Bool_t IsVtx = kFALSE;
833 const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
837 vertex->GetCovarianceMatrix(fCov);
838 if(vertex->GetNContributors() > 0) {
841 Double_t lvx = vertex->GetX();
842 Double_t lvy = vertex->GetY();
843 Double_t lvz = vertex->GetZ();
845 fEventCounter->Fill(5);
847 fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz);
849 if(TMath::Abs(lvx) < fVxMax) {
850 if(TMath::Abs(lvy) < fVyMax) {
851 if(TMath::Abs(lvz) < fVzMax) {
853 fEventCounter->Fill(6);
854 fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz);
862 }//Contributors check---
863 }//If vertex-----------
867 //---------------------------------------------------------------------------------------
869 //---------------------------------------------------------------------------------------
870 void AliEbyEHigherMomentsEffContTask::FillEffSparse(){
872 //For Efficiency-------------------------------------
874 Int_t nTracks = fAOD->GetNumberOfTracks();
875 Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
877 Int_t nMCTrack = fArrayMC->GetEntriesFast();
880 for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
882 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
885 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
889 /* TDatabasePDG *db = TDatabasePDG::Instance();
890 TParticlePDG *part = 0x0;
892 if (!(part = db->GetParticle(partMC->PdgCode()))) continue;
893 if ((part = db->GetParticle(partMC->PdgCode()))->Charge() == 0.0) continue;
895 if(!(partMC->PdgCode())) continue;
896 if(partMC->Charge() == 0) continue;
897 //-(1) Check for primary----
898 if(!partMC->IsPhysicalPrimary()) continue;
900 //-(2) Basic track cuts--------
902 if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
903 if (partMC->Pt() < fPtLowerLimit || partMC->Pt() > fPtHigherLimit) continue;
905 Double_t rap = partMC->Y();
907 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
909 if(TMath::Abs(partMC->GetPdgCode()) != gPdgCode) continue;
911 Short_t gCharge = partMC->Charge();
913 Float_t sign = (gCharge < 0 ) ? -1. : 1.;
915 Float_t findable = 1.;// Float_t(fHelper->IsParticleFindable(idxMC)); XXX
917 // (3) Initialize rec. variables-----
918 Float_t recStatus = 0.;
920 // (4) Initialize track parameters------
927 // (5) Loop over all Labels--------
930 for( Int_t iRec = 0; iRec < nTracks; iRec++ ){
932 if( iMC != fLabel[0][iRec] ) continue;
935 if( iMC == fLabel[1][iRec] )
938 AliAODTrack* aodTrack = NULL;
942 aodTrack = fAOD->GetTrack(iRec);
946 etaRec = aodTrack->Eta();
947 phiRec = aodTrack->Phi();
948 ptRec = aodTrack->Pt();
957 }//all charged track rec. check---
962 Double_t EffContainer[15] = {static_cast<Double_t>(fCentrality), partMC->Eta(), partMC->Y(), partMC->Phi(), partMC->Pt(), static_cast<Double_t>(sign),static_cast<Double_t>(findable), static_cast<Double_t>(recStatus), static_cast<Double_t>(recPid), static_cast<Double_t>(etaRec), static_cast<Double_t>(phiRec), static_cast<Double_t>(ptRec), partMC->Eta()-etaRec, partMC->Phi()-phiRec, partMC->Pt()-ptRec};
964 fTHnEff->Fill(EffContainer);
966 }//iMC Track loop-----
970 //---------------------------------------------------------------------------------------
972 //---------------------------------------------------------------------------------------
973 void AliEbyEHigherMomentsEffContTask::CheckContTrackAOD(Int_t label, Float_t sign, Int_t idxTrack){
975 AliAODMCParticle* particle = (AliAODMCParticle*)fArrayMC->At(label);
980 Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
981 Bool_t isSecondaryFromWeakDecay = kFALSE;
982 Bool_t isSecondaryFromMaterial = kFALSE;
984 // -- Check if correctly identified
985 if (particle->GetPdgCode() == (sign*gPdgCode)) {
987 // Check if is physical primary -> all ok
988 if(particle->IsPhysicalPrimary())
991 // -- Check if secondaries from material or weak decay
992 isSecondaryFromWeakDecay = particle->IsSecondaryFromWeakDecay();
993 isSecondaryFromMaterial = particle->IsSecondaryFromMaterial();
998 Float_t contSign = 0.;
1000 // TDatabasePDG *db = TDatabasePDG::Instance();
1001 //TParticlePDG *part = 0x0;
1003 if(( part = db->GetParticle(particle->PdgCode()))->Charge() == 0.) contSign = 0.;
1004 else if((part = db->GetParticle(particle->PdgCode()))->Charge() < 0.) contSign = -1.;
1005 else if((part = db->GetParticle(particle->PdgCode()))->Charge() > 0.) contSign = 1.;
1008 if(particle->Charge() == 0.) contSign = 0.;
1009 if(particle->Charge() < 0. ) contSign = -1;
1010 if(particle->Charge() > 0. ) contSign = 1.;
1012 // -- Get contaminating particle
1013 Float_t contPart = 0;
1014 if (isSecondaryFromWeakDecay) contPart = 7; // probeParticle from WeakDecay
1015 else if (isSecondaryFromMaterial) contPart = 8; // probeParticle from Material
1017 if (TMath::Abs(particle->GetPdgCode()) == 211) contPart = 1; // pion
1018 else if (TMath::Abs(particle->GetPdgCode()) == 321) contPart = 2; // kaon
1019 else if (TMath::Abs(particle->GetPdgCode()) == 2212) contPart = 3; // proton
1020 else if (TMath::Abs(particle->GetPdgCode()) == 11) contPart = 4; // electron
1021 else if (TMath::Abs(particle->GetPdgCode()) == 13) contPart = 5; // muon
1022 else contPart = 6; // other
1025 // -- Get Reconstructed values
1026 Float_t etaRec = 0.;
1027 Float_t phiRec = 0.;
1030 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(idxTrack));
1033 // if no track present (which should not happen)
1034 // -> pt = 0. , which is not in the looked at range
1036 // -- Get Reconstructed eta/phi/pt
1037 etaRec = aodTrack->Eta();
1038 phiRec = aodTrack->Phi();
1039 ptRec = aodTrack->Pt();
1042 // -- Fill THnSparse
1043 Double_t ContContainer[14] = { static_cast<Double_t>(fCentrality), particle->Eta(), particle->Y(), particle->Phi(), particle->Pt(), static_cast<Double_t>(sign), static_cast<Double_t>(contPart), static_cast<Double_t>(contSign), static_cast<Double_t>(etaRec), static_cast<Double_t>(phiRec), static_cast<Double_t>(ptRec), particle->Eta()-etaRec, particle->Phi()-phiRec, particle->Pt()-ptRec };
1045 fTHnCont->Fill(ContContainer);
1050 //------------------------------------------------------------------------
1051 void AliEbyEHigherMomentsEffContTask::Terminate( Option_t * ){
1053 Info("AliEbyEHigherMomentsEffContTask"," Task Successfully finished");
1057 //------------------------------------------------------------------------