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 //=========================================================================//
33 #include "AliAnalysisTaskSE.h"
34 #include "AliAnalysisManager.h"
36 #include "AliESDVertex.h"
37 #include "AliESDEvent.h"
38 #include "AliESDInputHandler.h"
39 #include "AliAODEvent.h"
40 #include "AliAODTrack.h"
41 #include "AliAODInputHandler.h"
43 #include "AliESDEvent.h"
44 #include "AliAODEvent.h"
46 #include "AliESDtrackCuts.h"
47 #include "AliAODMCHeader.h"
48 #include "AliAODMCParticle.h"
49 #include "AliMCEventHandler.h"
50 #include "AliMCEvent.h"
52 #include "AliGenHijingEventHeader.h"
53 #include "AliGenEventHeader.h"
55 #include "AliAODPid.h"
56 #include "AliPIDResponse.h"
57 #include "AliAODpidUtil.h"
58 #include "AliPIDCombined.h"
60 #include "AliEbyEHigherMomentsTask.h"
65 ClassImp(AliEbyEHigherMomentsTask)
68 //-----------------------------------------------------------------------
69 AliEbyEHigherMomentsTask::AliEbyEHigherMomentsTask( const char *name )
70 : AliAnalysisTaskSE( name ),
72 fCentralityEstimator("V0M"),
87 for(Int_t i = 0; i < 91; i++)
89 fhNplusNminus[i] = NULL;
92 for ( Int_t i = 0; i < 11; i++) {
96 DefineOutput(1, TList::Class()); // Outpput....
99 AliEbyEHigherMomentsTask::~AliEbyEHigherMomentsTask() {
100 if(fListOfHistos) delete fListOfHistos;
103 //---------------------------------------------------------------------------------
104 void AliEbyEHigherMomentsTask::UserCreateOutputObjects() {
106 fListOfHistos = new TList();
107 fListOfHistos->SetOwner(kTRUE);
108 fEventCounter = new TH1D("fEventCounter","EventCounter", 150, 0.5,150.5);
109 fListOfHistos->Add(fEventCounter);
111 fHistQA[0] = new TH1D("fHistQAvx", "Histo Vx After Cut", 400, -4., 4.);
112 fHistQA[1] = new TH1D("fHistQAvy", "Histo Vy After Cut", 400, -4., 4.);
113 fHistQA[2] = new TH1D("fHistQAvz", "Histo Vz After Cut", 500, -25., 25.);
114 fHistQA[3] = new TH1D("fHistQAvxA", "Histo Vx All", 500, -5., 5.);
115 fHistQA[4] = new TH1D("fHistQAvyA", "Histo Vy All", 500, -5., 5.);
116 fHistQA[5] = new TH1D("fHistQAvzA", "Histo Vz All", 500, -25., 25.);
117 fHistQA[6] = new TH1D("fHistQADcaXyC", "Histo DCAxy after cut", 500, -5., 5.);
118 fHistQA[7] = new TH1D("fHistQADcaZC", "Histo DCAz after cut", 500, -5., 5.);
119 fHistQA[8] = new TH1D("fHistQAPt","p_{T} distribution",600,0,6);
120 fHistQA[9] = new TH1D("fHistQAEta","#eta distribution",240,-1.2,1.2);
121 fHistQA[10] = new TH1D("fHistQAPhi","#phi distribution",340,0,6.8);
123 for(Int_t i = 0; i < 11; i++)
125 fListOfHistos->Add(fHistQA[i]);
132 for(Int_t i = 0; i < 25; i++) {
133 hname = "hPlusMinusCentBin"; hname+=i;
134 htitle = "N+ and N- in Cent Bin "; htitle+=i;
135 fhNplusNminus[i] = new TH2F(hname.Data(),htitle.Data(),1400,0.5,1400.5,1400,0.5,1400.5);
136 fListOfHistos->Add(fhNplusNminus[i]);
139 for(Int_t i = 25; i < 50; i++) {
140 hname = "hPlusMinusCentBin"; hname+=i;
141 htitle = "N+ and N- in Cent Bin "; htitle+=i;
142 fhNplusNminus[i] = new TH2F(hname.Data(),htitle.Data(),800,0.5,800.5,800,0.5,800.5);
143 fListOfHistos->Add(fhNplusNminus[i]);
146 for(Int_t i = 50; i < 91; i++) {
147 hname = "hPlusMinusCentBin"; hname+=i;
148 htitle = "N+ and N- in Cent Bin "; htitle+=i;
149 fhNplusNminus[i] = new TH2F(hname.Data(),htitle.Data(),500,0.5,500.5,500,0.5,500.5);
150 fListOfHistos->Add(fhNplusNminus[i]);
153 PostData(1, fListOfHistos);
157 //----------------------------------------------------------------------------------
158 void AliEbyEHigherMomentsTask::UserExec( Option_t * ){
159 fEventCounter->Fill(1);
160 Double_t nPlusCharge = 0.;
161 Double_t nMinusCharge = 0.;
163 AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent());
165 cout<< "ERROR 01: AOD not found " <<endl;
170 if(fAnalysisType == "AOD") {
171 AliAODHeader *aodHeader = gAOD->GetHeader();
172 Double_t cent = aodHeader->GetCentralityP()->GetCentralityPercentile(fCentralityEstimator.Data());
173 Double_t count = (Int_t)cent + 20;
174 fEventCounter->Fill(count);
176 if(cent < 0 || cent >= 99.99) return;
178 fEventCounter->Fill(2);
180 Int_t nCentrality = (Int_t)cent;
182 const AliAODVertex *vertex = gAOD->GetPrimaryVertex();
186 vertex->GetCovarianceMatrix(fCov);
187 if(vertex->GetNContributors() > 0){
190 Double_t lvx = vertex->GetX();
191 Double_t lvy = vertex->GetY();
192 Double_t lvz = vertex->GetZ();
194 fEventCounter->Fill(5);
196 fHistQA[3]->Fill(lvx);fHistQA[4]->Fill(lvy);fHistQA[5]->Fill(lvz);
198 if(TMath::Abs(lvx) < fVxMax) {
199 if(TMath::Abs(lvy) < fVyMax) {
200 if(TMath::Abs(lvz) < fVzMax) {
202 fEventCounter->Fill(6);
203 fHistQA[0]->Fill(lvx);fHistQA[1]->Fill(lvy);fHistQA[2]->Fill(lvz);
205 Int_t ntracks = gAOD->GetNumberOfTracks();
207 for(Int_t i = 0; i < ntracks; i++) {
208 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(i));
210 AliError(Form("ERROR: Could not retrieve AODtrack %d",i));
214 if(!aodTrack->TestFilterBit(nAODtrackCutBit)) continue;
218 dxy = aodTrack->DCA();
219 dz = aodTrack->ZAtDCA();
221 Double_t pt = aodTrack->Pt();
222 Double_t eta = aodTrack->Eta();
223 //Double_t nclus = aodTrack->GetTPCNcls();
224 Double_t nclus = aodTrack->GetTPCClusterInfo(2,1);//important for 2011 data
226 if( dxy > fDCAxy ) continue;
227 if( dz > fDCAz ) continue;
228 if( pt < fPtLowerLimit || pt > fPtHigherLimit ) continue;
229 if( eta < fEtaLowerLimit || eta > fEtaHigherLimit ) continue;
230 if( nclus < fTPCNClus ) continue;
232 fHistQA[8]->Fill(pt);
233 fHistQA[9]->Fill(eta);
234 fHistQA[10]->Fill(aodTrack->Phi());
235 fHistQA[6]->Fill(dxy);
236 fHistQA[7]->Fill(dz);
238 Short_t gCharge = aodTrack->Charge();
241 if(gCharge > 0) nPlusCharge += 1.;
242 if(gCharge < 0) nMinusCharge += 1.;
244 } //--------- Track Loop
245 // cout << "nCentrality "<<nCentrality<<" nPlusCharge "<< nPlusCharge << " nMinusCharge " << nMinusCharge << endl;
246 fhNplusNminus[nCentrality]->Fill(nPlusCharge,nMinusCharge);
247 fEventCounter->Fill(7);
254 }//AOD--analysis-----
256 else if(fAnalysisType == "MCAOD") {
257 TClonesArray *arrayMC = (TClonesArray*) gAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
259 AliFatal("Error: MC particles branch not found!\n");
261 AliAODMCHeader *mcHdr=0;
262 mcHdr=(AliAODMCHeader*)gAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
264 printf("MC header branch not found!\n");
268 AliAODHeader *aodHeader = gAOD->GetHeader();
269 AliCentrality* fcentrality = aodHeader->GetCentralityP();
270 Double_t cent =fcentrality ->GetCentralityPercentile("V0M");
272 Int_t nCentrality = (Int_t)cent;
273 if(cent < 0 || cent >= 91) return;
276 const AliAODVertex *vertex = gAOD->GetPrimaryVertex();
279 vertex->GetCovarianceMatrix(fCov);
280 if(vertex->GetNContributors() > 0) {
283 if(TMath::Abs(vertex->GetX()) < fVxMax) {
284 if(TMath::Abs(vertex->GetY()) < fVyMax) {
285 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
295 Int_t nMC = arrayMC->GetEntries();
296 for (Int_t iMC = 0; iMC < nMC; iMC++) {
297 AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);
298 if(!partMC->IsPhysicalPrimary())continue;
299 if (partMC->Eta() < fEtaLowerLimit || partMC->Eta() > fEtaHigherLimit) continue;
300 if (partMC->Pt() < fPtLowerLimit || partMC->Pt() > fPtHigherLimit) continue;
302 if(partMC->Charge() > 0) nPlusCharge += 1.;
303 if(partMC->Charge() < 0) nMinusCharge += 1.;
305 }//MC Track loop-----
306 fhNplusNminus[nCentrality]->Fill(nPlusCharge,nMinusCharge);
307 //cout << "Cent=" << nCentrality << " MC-PlusChrg=" << nPlusCharge << " MC-MinusChrg=" << nMinusCharge << endl;
309 }//MCAOD--analysis----
315 fEventCounter->Fill(8);
316 PostData(1, fListOfHistos);
321 //------------------------------------------------------------------------
322 void AliEbyEHigherMomentsTask::Terminate( Option_t * ){
324 Info("AliEbyEHigherMomentTask"," Task Successfully finished");
328 //------------------------------------------------------------------------