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 //
22 // Last Modified: 30/01/2014: only for net-charge part //
23 //=========================================================================//
31 #include "THnSparse.h"
34 #include "AliAnalysisTaskSE.h"
35 #include "AliAnalysisManager.h"
36 #include "AliAODEvent.h"
37 #include "AliVTrack.h"
38 #include "AliAODTrack.h"
39 #include "AliAODInputHandler.h"
40 #include "AliAODEvent.h"
42 #include "AliAODMCHeader.h"
43 #include "AliAODMCParticle.h"
44 #include "AliMCEventHandler.h"
45 #include "AliMCEvent.h"
47 #include "AliGenHijingEventHeader.h"
48 #include "AliGenEventHeader.h"
50 #include "AliAODPid.h"
51 #include "AliPIDResponse.h"
52 #include "AliAODpidUtil.h"
53 #include "AliPIDCombined.h"
54 #include "AliHelperPID.h"
56 #include "AliEbyEHigherMomentsTask.h"
61 ClassImp(AliEbyEHigherMomentsTask)
64 //-----------------------------------------------------------------------
65 AliEbyEHigherMomentsTask::AliEbyEHigherMomentsTask( const char *name )
66 : AliAnalysisTaskSE( name ),
70 fCentralityEstimator("V0M"),
84 fTHnCentNplusNminusCh(0),
85 fTHnCentNplusNminusChTruth(0),
86 fPtBinNplusNminusCh(0),
87 fPtBinNplusNminusChTruth(0)
90 for ( Int_t i = 0; i < 4; i++) {
94 DefineOutput(1, TList::Class()); // Outpput....
95 //DefineOutput(2, TList::Class());
98 AliEbyEHigherMomentsTask::~AliEbyEHigherMomentsTask() {
100 if(fListOfHistos) delete fListOfHistos;
103 //---------------------------------------------------------------------------------
104 void AliEbyEHigherMomentsTask::UserCreateOutputObjects() {
106 fListOfHistos = new TList();
107 fListOfHistos->SetOwner(kTRUE);
109 fEventCounter = new TH1D("fEventCounter","EventCounter", 10, 0.5,10.5);
110 fEventCounter->GetXaxis()->SetBinLabel(1,"Event Accesed");
111 fEventCounter->GetXaxis()->SetBinLabel(2,"Within 0-80% centrality");
112 fEventCounter->GetXaxis()->SetBinLabel(5,"Have a vertex");
113 fEventCounter->GetXaxis()->SetBinLabel(6,"After vertex Cut");
114 fEventCounter->GetXaxis()->SetBinLabel(7,"Event Analyzed");
115 fListOfHistos->Add(fEventCounter);
118 fHistQA[0] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
119 fHistQA[1] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
120 fHistQA[2] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
121 fHistQA[3] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
123 for(Int_t i = 0; i < 4; i++)
125 fListOfHistos->Add(fHistQA[i]);
128 fHistVxVy = new TH2D("fHistVxVy","Vertex-x Vs Vertex-y", 200, -1., 1., 200, -1., 1.);
129 fListOfHistos->Add(fHistVxVy);
132 const Int_t nDim = 3;
134 Int_t fBinsCh[nDim] = {100, 1500, 1500};
135 Double_t fMinCh[nDim] = { -0.5, -0.5, -0.5 };
136 Double_t fMaxCh[nDim] = { 99.5, 1499.5, 1499.5};
138 const Int_t dim = fNptBins*2 + 1;
139 //const Int_t dim = ;
143 bin[0] = 100; min[0] = -0.5; max[0] = 99.5;
145 for(Int_t i = 1; i < dim; i++){
153 fTHnCentNplusNminusCh = new THnSparseD("fTHnCentNplusNminusCh","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
154 fTHnCentNplusNminusCh->GetAxis(0)->SetTitle("Centrality");
155 fTHnCentNplusNminusCh->GetAxis(1)->SetTitle("Nplus");
156 fTHnCentNplusNminusCh->GetAxis(2)->SetTitle("Nminus");
157 fListOfHistos->Add(fTHnCentNplusNminusCh);
159 fPtBinNplusNminusCh = new THnSparseI("fPtBinNplusNminusCh","cent-nplus-nminus", dim, bin, min, max);
161 fListOfHistos->Add(fPtBinNplusNminusCh);
163 if( fAnalysisType == "MCAOD" ){
165 fTHnCentNplusNminusChTruth = new THnSparseD("fTHnCentNplusNminusChTruth","Cent-NplusChrg-NminusChrg", nDim, fBinsCh, fMinCh, fMaxCh);
166 fTHnCentNplusNminusChTruth->GetAxis(0)->SetTitle("Centrality");
167 fTHnCentNplusNminusChTruth->GetAxis(1)->SetTitle("Nplus");
168 fTHnCentNplusNminusChTruth->GetAxis(2)->SetTitle("Nminus");
169 fListOfHistos->Add(fTHnCentNplusNminusChTruth);
171 fPtBinNplusNminusChTruth = new THnSparseI("fPtBinNplusNminusChTruth","cent-nplus-nminus", dim, bin, min, max);
172 fListOfHistos->Add(fPtBinNplusNminusChTruth);
178 //PostData(1, fListOfHistosQA);
179 PostData(1, fListOfHistos);
185 //----------------------------------------------------------------------------------
186 void AliEbyEHigherMomentsTask::UserExec( Option_t * ){
188 fEventCounter->Fill(1);
190 if(fAnalysisType == "AOD") {
194 }//AOD--analysis-----
196 else if(fAnalysisType == "MCAOD") {
208 //--------------------------------------------------------------------------------------
209 void AliEbyEHigherMomentsTask::doAODEvent(){
211 Double_t positiveSum = 0.;
212 Double_t negativeSum = 0.;
213 const Int_t dim = fNptBins*2;
216 for(Int_t idx = 0; idx < dim; idx++){
221 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
223 cout<<"ERROR: Analysis manager not found."<<endl;
226 //coneect to the inputHandler------------
227 AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
229 cout<<"ERROR: Input handler not found."<<endl;
233 AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
237 cout<< "ERROR: AOD not found " <<endl;
241 if(!ProperVertex(fAOD)) return;
243 AliAODHeader *aodHeader = fAOD->GetHeader();
245 fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
247 if(fCentrality < 0 || fCentrality >= 81) return;
249 fEventCounter->Fill(2);
251 Int_t nTracks = fAOD->GetNumberOfTracks();
253 for(Int_t i = 0; i < nTracks; i++) {
255 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
258 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
262 if(!AcceptTrack(aodTrack) ) continue;
265 fBin = GetPtBin(aodTrack->Pt());
267 Short_t gCharge = aodTrack->Charge();
274 PtCh[fNptBins+fBin] += 1;
278 }//--------- Track Loop to select with filterbit
280 Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), positiveSum, negativeSum};
282 fTHnCentNplusNminusCh->Fill(fContainerCh);
284 //cout << fCentrality << " " << positiveSum <<" " << negativeSum << endl;
286 Double_t ptContainer[dim+1];
288 ptContainer[0] = fCentrality;
290 for(Int_t i = 1; i <= dim; i++){
291 ptContainer[i] = PtCh[i-1];
292 //cout << PtCh[i-1] <<" ,";
296 fPtBinNplusNminusCh->Fill(ptContainer);
299 fEventCounter->Fill(7);
300 PostData(1, fListOfHistos);
304 //--------------------------------------------------------------------------------------
305 //--------------------------------------------------------------------------------------
306 void AliEbyEHigherMomentsTask::doMCAODEvent(){
310 Double_t positiveSumMCRec = 0.;
311 Double_t negativeSumMCRec = 0.;
313 Double_t positiveSumMCTruth = 0.;
314 Double_t negativeSumMCTruth = 0.;
316 const Int_t dim = fNptBins*2;
318 Double_t ptChMC[dim];
319 for(Int_t idx = 0; idx < dim; idx++){
324 //Connect to Anlaysis manager------
325 AliAnalysisManager* manager = AliAnalysisManager::GetAnalysisManager();
327 cout<<"ERROR: Analysis manager not found."<<endl;
331 AliAODInputHandler* inputHandler = dynamic_cast<AliAODInputHandler*> (manager->GetInputEventHandler());
333 cout<<"ERROR: Input handler not found."<<endl;
338 AliAODEvent *fAOD = dynamic_cast<AliAODEvent*>(InputEvent());
342 cout<< "ERROR: AOD not found " <<endl;
347 // ------------------------------------------------------------------
349 fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
351 AliFatal("Error: MC particles branch not found!\n");
355 AliAODMCHeader *mcHdr=NULL;
356 mcHdr=(AliAODMCHeader*)fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
358 printf("MC header branch not found!\n");
362 if(!ProperVertex(fAOD)) return;
364 AliAODHeader *aodHeader = fAOD->GetHeader();
366 fCentrality = (Int_t)aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
368 if( fCentrality < 0 || fCentrality >= 81) return;
370 fEventCounter->Fill(2);
372 Int_t nTracks = fAOD->GetNumberOfTracks();
374 for(Int_t i = 0; i < nTracks; i++) {
376 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(fAOD->GetTrack(i));
379 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
383 if(!AcceptTrack(aodTrack) ) continue;
386 fBin = GetPtBin(aodTrack->Pt());
388 //cout << "Pt Bin " << fBin << endl;
390 Short_t gCharge = aodTrack->Charge();
392 positiveSumMCRec += 1.;
396 negativeSumMCRec += 1.;
397 PtCh[fNptBins+fBin] += 1;
400 }//--------- Track Loop to select with filterbit
402 //===============================================================================
403 //---------------------MC Truth--------------------------------------------------
404 //===============================================================================
406 Int_t nMCTrack = fArrayMC->GetEntriesFast();
408 for (Int_t iMC = 0; iMC < nMCTrack; iMC++) {
410 AliAODMCParticle *partMC = (AliAODMCParticle*) fArrayMC->At(iMC);
413 AliError(Form("ERROR: Could not retrieve AODMCtrack %d",iMC));
417 if(partMC->Charge() == 0) continue;
418 if(!partMC->IsPhysicalPrimary()) continue;
420 if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
421 if (partMC->Pt() < fPtLowerLimit || partMC->Pt() > fPtHigherLimit) continue;
423 Short_t gCharge = partMC->Charge();
424 Int_t mcbin = GetPtBin(partMC->Pt());
427 positiveSumMCTruth += 1.;
431 negativeSumMCTruth += 1.;
432 ptChMC[fNptBins+mcbin] += 1.;
436 }//MC-Truth Track loop--
438 Double_t fContainerCh[3] = { static_cast<Double_t>(fCentrality), positiveSumMCRec, negativeSumMCRec};//Reco. values ch. hadrons
439 Double_t fContainerChTruth[3] = { static_cast<Double_t>(fCentrality), positiveSumMCTruth, negativeSumMCTruth};
441 fTHnCentNplusNminusCh->Fill(fContainerCh);//Fill the rec. ch. particles---
442 fTHnCentNplusNminusChTruth->Fill(fContainerChTruth);//MC -Truth ch. particles
444 //cout << fCentrality << " " << positiveSumMCRec << " " << negativeSumMCRec <<endl;
445 //cout <<fCentrality<<" " << positiveSumMCTruth << " " << negativeSumMCTruth << endl;
446 //cout <<" " << posPidSumMCRec << " " << negPidSumMCRec << endl;
447 //cout <<" " << posPidSumMCTruth << " " << negPidSumMCTruth << endl;
448 //cout <<"---------------------------------" << endl;
450 Double_t ptContainer[dim+1];
451 Double_t ptContainerMC[dim+1];
452 ptContainer[0] = fCentrality;
453 ptContainerMC[0] = fCentrality;
455 for(Int_t i = 1; i <= dim; i++){
456 ptContainer[i] = PtCh[i-1];
457 ptContainerMC[i] = ptChMC[i-1];
458 //cout <<" "<< PtCh[i-1]<<endl;
459 // cout<< " MC=" << ptChMC[i-1];
464 fPtBinNplusNminusCh->Fill(ptContainer);
465 fPtBinNplusNminusChTruth->Fill(ptContainerMC);
467 fEventCounter->Fill(7);
468 PostData(1, fListOfHistos);
473 //---------------------------------------------------------------------------------------
474 Bool_t AliEbyEHigherMomentsTask::ProperVertex(AliAODEvent *fAOD) const{
476 Bool_t IsVtx = kFALSE;
478 const AliAODVertex *vertex = fAOD->GetPrimaryVertex();
482 vertex->GetCovarianceMatrix(fCov);
483 if(vertex->GetNContributors() > 0) {
486 Double_t lvx = vertex->GetX();
487 Double_t lvy = vertex->GetY();
488 Double_t lvz = vertex->GetZ();
490 fEventCounter->Fill(5);
492 if(TMath::Abs(lvx) < fVxMax) {
493 if(TMath::Abs(lvy) < fVyMax) {
494 if(TMath::Abs(lvz) < fVzMax) {
496 fEventCounter->Fill(6);
497 fHistQA[0]->Fill(lvz);
498 fHistVxVy->Fill(lvx,lvy);
505 }//Contributors check---
506 }//If vertex-----------
508 AliCentrality *centrality = fAOD->GetCentrality();
509 if (centrality->GetQuality() != 0) IsVtx = kFALSE;
513 //---------------------------------------------------------------------------------------
515 //---------------------------------------------------------------------------------------
516 Bool_t AliEbyEHigherMomentsTask::AcceptTrack(AliAODTrack* track) const{
518 if(!track) return kFALSE;
519 if( track->Charge() == 0 ) return kFALSE;
521 Double_t pt = track->Pt();
522 Double_t eta = track->Eta();
523 if(!track->TestFilterBit(fAODtrackCutBit)) return kFALSE;
524 if( pt < fPtLowerLimit || pt > fPtHigherLimit ) return kFALSE;
525 if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) return kFALSE;
527 fHistQA[1]->Fill(pt);
528 fHistQA[2]->Fill(eta);
529 fHistQA[3]->Fill(track->Phi());
533 //------------------------------------------------------------------------
535 //------------------------------------------------------------------------
536 Int_t AliEbyEHigherMomentsTask::GetPtBin(Double_t pt){
540 Double_t BinSize = (fPtHigherLimit - fPtLowerLimit)/fNptBins;
542 for(Int_t iBin = 0; iBin < fNptBins; iBin++){
544 Double_t xLow = fPtLowerLimit + iBin*BinSize;
545 Double_t xHigh = fPtLowerLimit + (iBin+1)*BinSize;
547 if( pt >= xLow && pt < xHigh){
549 //cout << "Pt Bin #" << bin <<"pt=" << pt << endl;
558 //------------------------------------------------------------------------
559 //------------------------------------------------------------------------
560 void AliEbyEHigherMomentsTask::Terminate( Option_t * ){
562 Info("AliEbyEHigherMomentTask"," Task Successfully finished");
566 //------------------------------------------------------------------------