1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: Nirbhay K. Behera *
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 // Toy Model 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 "AliHigherMomentsToyModel.h"
71 ClassImp(AliHigherMomentsToyModel)
74 //-----------------------------------------------------------------------
75 AliHigherMomentsToyModel::AliHigherMomentsToyModel( const char *name )
76 : AliAnalysisTaskSE( name ),
82 fParticleSpecies(AliPID::kProton),
84 fCentralityEstimator("V0M"),
99 fAODtrackCutBit(1024),
105 fTHnCentNplusNminusCh(0),
106 fTHnCentNplusNminusChTruth(0),
107 fTHnCentNplusNminus(0)
110 for ( Int_t i = 0; i < 13; i++) {
114 for ( Int_t i = 0; i < 5; i++ ){
115 fTHnCentNplusNminusPid[i] = NULL;
116 fTHnCentNplusNminusPidTruth[i] = NULL;
119 DefineOutput(1, TList::Class()); // Outpput....
120 DefineOutput(2, TList::Class());
124 AliHigherMomentsToyModel::~AliHigherMomentsToyModel() {
126 if(fListOfHistosQA) delete fListOfHistosQA;
127 if(fListOfHistos) delete fListOfHistos;
131 //---------------------------------------------------------------------------------
132 void AliHigherMomentsToyModel::UserCreateOutputObjects() {
134 fListOfHistosQA = new TList();
135 fListOfHistosQA->SetOwner(kTRUE);
137 fListOfHistos = new TList();
138 fListOfHistos->SetOwner(kTRUE);
140 fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
141 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
142 fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-90% centrality");
143 fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
144 fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
145 fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
146 fEventCounter->GetXaxis()->SetBinLabel(8,"Event Analysis finished");
147 fListOfHistosQA->Add(fEventCounter);
149 fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.);
150 fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.);
151 fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
152 fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.);
153 fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.);
154 fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.);
155 fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.);
156 fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.);
157 fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
158 fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
159 fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
160 fHistQA[11] = new TH1D("fHistQANCls","Number of TPC cluster",200,0,200);
161 fHistQA[12] = new TH1D("fHistQAChi2","Chi2 per NDF",100,0,10);
163 for(Int_t i = 0; i < 13; i++)
165 fListOfHistosQA->Add(fHistQA[i]);
168 fHistDCA = new TH2D("fHistDCA","DCAxy Vs DCAz", 500, -5., 5., 500, -5., 5.);
169 fTPCSig = new TH2D("fTPCSig","TPC signal",200, 0.0, 10. ,1000,0.,1000);
170 fTPCSig->SetMarkerColor(kRed);
171 fTPCSigA = new TH2D("fTPCSigA","TPC signal all ",200, 0.0, 10. ,1000,0.,1000);
172 fListOfHistosQA->Add(fHistDCA);
173 fListOfHistosQA->Add(fTPCSig);
174 fListOfHistosQA->Add(fTPCSigA);
177 const Int_t nDim = 3;
178 const Int_t nPid = 5;
180 TString Species[nPid] = {"Electron","Muon","Pion","Kaon","Proton"};
182 Int_t fBins[nPid][nDim];
183 Double_t fMin[nPid][nDim];
184 Double_t fMax[nPid][nDim];
186 for( Int_t i = 0; i < nPid; i++ ){
187 for( Int_t j = 0; j < nDim; j++ ){
194 Int_t fBinsCh[nDim] = {100, 1500, 1500};
195 Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
196 Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
198 fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
199 fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
200 fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
201 fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
202 fListOfHistos->Add(fTHnCentNplusNminusCh);
204 if( fAnalysisType == "MCAOD" ){
206 fTHnCentNplusNminusChTruth = new THnSparseD("fTHnCentNplusNminusChTruth","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
207 fTHnCentNplusNminusChTruth->GetAxis(0)->SetTitle("Centrality");
208 fTHnCentNplusNminusChTruth->GetAxis(1)->SetTitle("Nplus");
209 fTHnCentNplusNminusChTruth->GetAxis(2)->SetTitle("Nminus");
210 fListOfHistos->Add(fTHnCentNplusNminusChTruth);
214 TString hname1, hname11;
215 TString htitle1, axisTitle1, axisTitle2;
221 Int_t gPid = (Int_t)fParticleSpecies;
223 if( gPid > 1 && gPid < 5 ){
225 fBins[2][0] = 100; fBins[2][1] = 1000; fBins[2][2] = 1000;
226 fMin[2][0] = -0.5; fMin[2][1] = -0.5; fMin[2][2] = -0.5;
227 fMax[2][0] = 99.5; fMax[2][1] = 999.5; fMax[2][2] = 999.5;
229 fBins[3][0] = 100; fBins[3][1] = 600; fBins[3][2] = 600;
230 fMin[3][0] = -0.5; fMin[3][1] = -0.5; fMin[3][2] = -0.5;
231 fMax[3][0] = 99.5; fMax[3][1] = 599.5; fMax[3][2] = 599.5;
233 fBins[4][0] = 100; fBins[4][1] = 400; fBins[4][2] = 400;
234 fMin[4][0] = -0.5; fMin[4][1] = -0.5; fMin[4][2] = -0.5;
235 fMax[4][0] = 99.5; fMax[4][1] = 399.5; fMax[4][2] = 399.5;
237 hname1 = "fCentNplusNminusPid"; hname1 +=gPid;
238 htitle1 = Species[gPid]; htitle1 +=" And Neg-"; htitle1 +=Species[gPid];
239 axisTitle1 = Species[gPid];
240 axisTitle2 ="Neg-"; axisTitle2 += Species[gPid];
242 fTHnCentNplusNminusPid[gPid] = new THnSparseD(hname1.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
244 fTHnCentNplusNminusPid[gPid]->GetAxis(0)->SetTitle("Centrality");
245 fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
246 fTHnCentNplusNminusPid[gPid]->GetAxis(1)->SetTitle(axisTitle2.Data());
248 fListOfHistos->Add(fTHnCentNplusNminusPid[gPid]);
250 if( fAnalysisType == "MCAOD" ){
252 hname11 = "fCentNplusNminusPidTruth"; hname11 +=gPid;
253 fTHnCentNplusNminusPidTruth[gPid] = new THnSparseD(hname11.Data(),htitle1.Data(), nDim, fBins[gPid], fMin[gPid],fMax[gPid]);
255 fTHnCentNplusNminusPidTruth[gPid]->GetAxis(0)->SetTitle("Centrality");
256 fTHnCentNplusNminusPidTruth[gPid]->GetAxis(1)->SetTitle(axisTitle1.Data());
257 fTHnCentNplusNminusPidTruth[gPid]->GetAxis(1)->SetTitle(axisTitle2.Data());
259 fListOfHistos->Add(fTHnCentNplusNminusPidTruth[gPid]);
262 }//Pion, Koan and Proton-------
266 Int_t fBinsX[nDim] = {100, 1500, 1500};
267 Double_t fMinX[nDim] = { -0.5, -0.5, -0.5 };
268 Double_t fMaxX[nDim] = { 99.5, 1499.5, 1499.5};
270 fTHnCentNplusNminus = new THnSparseD("fTHnCentNplusNminus","Cent-NplusChrg-NminusChrg", nDim, fBinsX, fMinX, fMaxX);
271 fTHnCentNplusNminus->GetAxis(0)->SetTitle("Centrality");
272 fTHnCentNplusNminus->GetAxis(1)->SetTitle("UnwantedPlus");
273 fTHnCentNplusNminus->GetAxis(2)->SetTitle("UnwantedMinus");
274 fListOfHistos->Add(fTHnCentNplusNminus);
276 }//Unwanted particle -------
282 PostData(1, fListOfHistosQA);
283 PostData(2, fListOfHistos);
288 //----------------------------------------------------------------------------------
289 void AliHigherMomentsToyModel::UserExec( Option_t * ){
291 fEventCounter->Fill(1);
294 if(fAnalysisType == "AOD") {
298 }//AOD--analysis-----
300 else if(fAnalysisType == "MCAOD") {
311 fEventCounter->Fill(8);
313 PostData(1, fListOfHistosQA);
314 PostData(2, fListOfHistos);
318 //--------------------------------------------------------------------------------------
319 void AliHigherMomentsToyModel::doAODEvent(){
321 //-------------------
322 Double_t nPlusCharge = 0.;
323 Double_t nMinusCharge = 0.;
324 Double_t nPartile = 0.;
325 Double_t nAntiParticle = 0.;
329 //connect to the analysis mannager-----
331 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
333 cout<<"ERROR: Analysis manager not found."<<endl;
336 //coneect to the inputHandler------------
337 AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
339 cout<<"ERROR: Input handler not found."<<endl;
343 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
346 cout<< "ERROR 01: AOD not found " <<endl;
350 fPIDResponse =inputHandler->GetPIDResponse();
354 AliFatal("This Task needs the PID response attached to the inputHandler");
359 AliAODHeader *aodHeader = fAOD->GetHeader();
361 fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
363 cent = aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data());
366 fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data());
367 else if (cent == 10 || cent == -1.)
369 else if (cent > 0 && cent < 9)
370 fCentrality = cent + 1;
372 if(fCentrality < 0 || fCentrality >= 91) return;
374 fEventCounter->Fill(2);
376 if(!ProperVertex()) return;
378 Int_t nTracks = fAOD->GetNumberOfTracks();
380 TExMap *trackMap = new TExMap();//Mapping matrix----
382 for(Int_t i = 0; i < nTracks; i++) {
384 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
387 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
391 Double_t tpcSignalAll = aodTrack->GetTPCsignal();
392 fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
394 Int_t gID = aodTrack->GetID();
396 if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
397 }//1st track loop----
399 AliAODTrack* newAodTrack;
401 for( Int_t j = 0; j < nTracks; j++ )
404 AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
407 AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
412 if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
414 Int_t gID = aodTrack1->GetID();
416 //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
417 newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
419 Float_t dxy = 0., dz = 0.;
421 dxy = aodTrack1->DCA();
422 dz = aodTrack1->ZAtDCA();
424 Double_t pt = aodTrack1->Pt();
425 Double_t eta = aodTrack1->Eta();
426 Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
427 Double_t chi2ndf = aodTrack1->Chi2perNDF();
430 if( fabs(dxy) > fDCAxy ) continue;
431 if( fabs(dz) > fDCAz ) continue;
432 //Extra cut on DCA---( Similar to BF Task.. )
433 if( fDCAxy !=0 && fDCAz !=0 ){
434 if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
437 if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
438 if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
439 if( nclus < fTPCNClus ) continue;
440 if( chi2ndf > fChi2perNDF ) continue;
443 fHistQA[6]->Fill(dxy);
444 fHistQA[7]->Fill(dz);
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 fHistDCA->Fill(dxy,dz);
452 Short_t gCharge = aodTrack1->Charge();
454 if(gCharge > 0) nPlusCharge += 1.;
455 if(gCharge < 0) nMinusCharge += 1.;
459 gPid = (Int_t)fParticleSpecies;
461 Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
462 Double_t tpcSignal = newAodTrack->GetTPCsignal();
463 //Double_t rap = aodTrack1->Y(AliPID::ParticleMass(fParticleSpecies));
464 //Double_t tpcSignal = aodTrack1->GetTPCsignal();
466 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
468 fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
470 Float_t nsigmaTPCPID = -999.;
471 Float_t nsigmaTOFPID = -999.;
472 //Float_t nsigmaTPCTOFPID = -999.;
474 nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
475 nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
477 if ( nsigmaTPCPID < fNSigmaCut ){
479 if (gCharge > 0) nPartile +=1.;
480 if( gCharge < 0 ) nAntiParticle +=1.;
485 }//--------- Track Loop to select with filterbit
489 Double_t fContainerCh[3] = { fCentrality, nPlusCharge, nMinusCharge};
490 Double_t fContainerPid[3] = { fCentrality, nPartile, nAntiParticle};
493 fTHnCentNplusNminusCh->Fill(fContainerCh);
496 gPid = (Int_t)fParticleSpecies;
497 fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
499 // cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl;
503 //cout << "nCentrality "<< fCentrality <<", nPlusCharge="<< nPlusCharge << ", nMinusCharge=" << nMinusCharge << endl;
505 fEventCounter->Fill(7);
510 //--------------------------------------------------------------------------------------
513 //--------------------------------------------------------------------------------------
514 void AliHigherMomentsToyModel::doMCAODEvent(){
518 Double_t nPlusCharge = 0.;
519 Double_t nMinusCharge = 0.;
521 Double_t nPlusChargeTruth = 0.;
522 Double_t nMinusChargeTruth = 0.;
524 Double_t nPartile = 0.;
525 Double_t nAntiParticle = 0.;
526 Double_t nPartileTruth = 0.;
527 Double_t nAntiParticleTruth = 0.;
530 Int_t gPdgCode = AliPID::ParticleCode(fParticleSpecies);
532 //Connect to Anlaysis manager------
533 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
535 cout<<"ERROR: Analysis manager not found."<<endl;
539 AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
541 cout<<"ERROR: Input handler not found."<<endl;
546 fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
549 cout<< "ERROR 01: AOD not found " <<endl;
553 fPIDResponse =inputHandler->GetPIDResponse();
557 AliFatal("This Task needs the PID response attached to the inputHandler");
561 // ------------------------------------------------------------------
563 fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
565 AliFatal("Error: MC particles branch not found!\n");
569 AliAODMCHeader *mcHdr=NULL;
570 mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
572 printf("MC header branch not found!\n");
576 AliAODHeader *aodHeader = fAOD->GetHeader();
578 fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
582 cent = aodHeader->GetCentralityP()->GetCentralityClass10(fCentralityEstimator.Data());
585 fCentrality = aodHeader->GetCentralityP()->GetCentralityClass5(fCentralityEstimator.Data());
586 else if (cent == 10 || cent == -1.)
588 else if (cent > 0 && cent < 9)
589 fCentrality = cent + 1;
591 if( fCentrality < 0 || fCentrality >= 91) return;
593 fEventCounter->Fill(2);
597 if(!ProperVertex()) return;
599 Int_t nTracks = fAOD->GetNumberOfTracks();
602 TExMap *trackMap = new TExMap();//Mapping matrix----
604 for(Int_t i = 0; i < nTracks; i++) {
606 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
609 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
613 Double_t tpcSignalAll = aodTrack->GetTPCsignal();
614 fTPCSigA->Fill(aodTrack->GetTPCmomentum(),tpcSignalAll);
616 Int_t gID = aodTrack->GetID();
618 if (aodTrack->TestFilterBit(1)) trackMap->Add(gID, i);//Global tracks-----
620 }//1st track loop----
622 AliAODTrack* newAodTrack;
624 for( Int_t j = 0; j < nTracks; j++ ){
626 AliAODTrack* aodTrack1 = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(j));
629 AliError(Form("ERROR: Could not retrieve AODtrack %d",j));
634 if(!aodTrack1->TestFilterBit(fAODtrackCutBit)) continue;
637 Int_t gID = aodTrack1->GetID();
639 //if( aodTrack->GetID() != (-aodTrack1->GetID() -1) ) continue;
641 newAodTrack = gID >= 0 ? aodTrack1 : fAOD->GetTrack(trackMap->GetValue(-1-gID)); //Take those global track who corresponds to TPC only track
644 //cout << dxy << endl;
645 Double_t pt = aodTrack1->Pt();
646 Double_t eta = aodTrack1->Eta();
647 Double_t nclus = aodTrack1->GetTPCClusterInfo(2,1);
648 Double_t chi2ndf = aodTrack1->Chi2perNDF();
650 Float_t dxy = 0., dz = 0.;
652 if( fAODtrackCutBit == 128 ){
654 dxy = aodTrack1->DCA();
655 dz = aodTrack1->ZAtDCA();
657 if( fabs(dxy) > fDCAxy ) continue;
658 if( fabs(dz) > fDCAz ) continue;
659 //Extra cut on DCA---( Similar to BF Task.. )
660 if( fDCAxy !=0 && fDCAz !=0 ){
661 if( TMath::Sqrt((dxy*dxy)/(fDCAxy*fDCAxy)+(dz*dz)/(fDCAz*fDCAz)) > 1. ) continue;
664 fHistQA[6]->Fill(dxy);
665 fHistQA[7]->Fill(dz);
666 fHistDCA->Fill(dxy,dz);
671 if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
672 if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
673 if( nclus < fTPCNClus ) continue;
674 if( chi2ndf > fChi2perNDF ) continue;
678 fHistQA[8]->Fill(pt);
679 fHistQA[9]->Fill(eta);
680 fHistQA[10]->Fill(aodTrack1->Phi());
681 fHistQA[11]->Fill(nclus);
682 fHistQA[12]->Fill(chi2ndf);
685 Short_t gCharge = aodTrack1->Charge();
687 if( gCharge == 0 ) continue;
690 if(gCharge > 0) nPlusCharge += 1.;
691 if(gCharge < 0) nMinusCharge += 1.;
695 gPid = (Int_t)fParticleSpecies;;
696 Double_t rap = newAodTrack->Y(AliPID::ParticleMass(fParticleSpecies));
697 Double_t tpcSignal = newAodTrack->GetTPCsignal();
699 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
701 fTPCSig->Fill(newAodTrack->GetTPCmomentum(),tpcSignal);
703 Float_t nsigmaTPCPID = -999.;
704 Float_t nsigmaTOFPID = -999.;
705 //Float_t nsigmaTPCTOFPID = -999.;
707 nsigmaTPCPID = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(newAodTrack,fParticleSpecies));
708 nsigmaTOFPID = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(newAodTrack,fParticleSpecies));
710 if ( nsigmaTPCPID < fNSigmaCut ){
713 if (gCharge > 0) nPartile +=1.;
714 if( gCharge < 0 ) nAntiParticle +=1.;
719 }//--------- Track Loop to select with filterbit
721 //cout << "Cent=" << fCentrality << " MC-PlusChrg=" << nPlusCharge << " MC-MinusChrg=" << nMinusCharge << endl;
723 //cout << "nCentrality "<< fCentrality <<", nParticle="<< nPartile << ", nMinusParticle=" << nAntiParticle << endl;
725 Double_t fContainerCh[3] = { fCentrality, nPlusCharge, nMinusCharge};
726 Double_t fContainerPid[3] = { fCentrality, nPartile, nAntiParticle};
728 fTHnCentNplusNminusCh->Fill(fContainerCh);
732 gPid = (Int_t)fParticleSpecies;
733 fTHnCentNplusNminusPid[gPid]->Fill(fContainerPid);
737 //===========================================
738 //--------MC Truth---------------------------
739 //===========================================
741 Int_t nMCTrack = fArrayMC->GetEntriesFast();
743 for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
745 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
748 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
752 if(partMC->Charge() == 0) continue;
753 if(!partMC->IsPhysicalPrimary()) continue;
755 if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
756 if (partMC->Pt() < fPtLowerLimit || partMC->Pt() > fPtHigherLimit) continue;
758 Short_t gCharge = partMC->Charge();
760 if(gCharge > 0) nPlusChargeTruth += 1.;
761 if(gCharge < 0) nMinusChargeTruth += 1.;
765 Double_t rap = partMC->Y();
767 if( fabs(rap ) > fRapidityCut ) continue;//Rapidity cut
769 if(TMath::Abs(partMC->GetPdgCode()) != gPdgCode) continue;
771 if(gCharge > 0) nPartileTruth += 1.;
772 if(gCharge < 0) nAntiParticleTruth += 1.;
776 }//MC-Truth Track loop--
778 Double_t fContainerChTruth[3] = { fCentrality, nPlusChargeTruth, nMinusChargeTruth };
779 Double_t fContainerPidTruth[3] = { fCentrality, nPartileTruth, nAntiParticleTruth };
781 //cout << "Cent=" << fCentrality << " MC-PlusChrgT=" << nPlusChargeTruth << " MC-MinusChrgT=" << nMinusChargeTruth << endl;
783 //cout << "nCentrality "<< fCentrality <<", nParticleT="<< nPartileTruth << ", nMinusParticleT=" << nAntiParticleTruth << endl;
785 //cout << "----------------------------" << endl;
786 fTHnCentNplusNminusChTruth->Fill(fContainerChTruth);
789 gPid = (Int_t)fParticleSpecies;
790 fTHnCentNplusNminusPidTruth[gPid]->Fill(fContainerPidTruth);
795 fEventCounter->Fill(7);
800 //---------------------------------------------------------------------------------------
802 //---------------------------------------------------------------------------------------
803 Bool_t AliHigherMomentsToyModel::ProperVertex(){
805 Bool_t IsVtx = kFALSE;
807 const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
811 vertex->GetCovarianceMatrix(fCov);
812 if(vertex->GetNContributors() > 0) {
815 Double_t lvx = vertex->GetX();
816 Double_t lvy = vertex->GetY();
817 Double_t lvz = vertex->GetZ();
819 fEventCounter->Fill(5);
821 fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz);
823 if(TMath::Abs(lvx) < fVxMax) {
824 if(TMath::Abs(lvy) < fVyMax) {
825 if(TMath::Abs(lvz) < fVzMax) {
827 fEventCounter->Fill(6);
828 fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz);
836 }//Contributors check---
837 }//If vertex-----------
841 //---------------------------------------------------------------------------------------
845 //------------------------------------------------------------------------
846 void AliHigherMomentsToyModel::Terminate( Option_t * ){
848 Info("AliHigherMomentsToyModel"," Task Successfully finished");
852 //------------------------------------------------------------------------