9 #include "AliAnalysisTask.h"
10 #include "AliAnalysisManager.h"
13 #include "AliMCEvent.h"
14 #include "AliESDEvent.h"
15 #include "AliESDInputHandler.h"
16 #include "AliESDtrackCuts.h"
17 #include "AliCentrality.h"
19 #include "AliEbyEFluctuationAnalysisTask.h"
21 // Event by event charge fluctuation analysis
22 // Authors: Satyajit Jena and Panos Cristakoglou
25 ClassImp(AliEbyEFluctuationAnalysisTask)
27 //________________________________________________________________________
28 AliEbyEFluctuationAnalysisTask::AliEbyEFluctuationAnalysisTask(const char *name)
29 : AliAnalysisTaskSE(name), fESD(0), fOutputList(0),
30 fHistEventStats(0), fHistCentrality(0),
31 fHistNMultMC(0), fHistNPlusNMinusMC(0),
33 fAnalysisType("ESD"), fAnalysisMode("TPC"),
34 fCentralityEstimator("V0M"), fCentralityBins20(kFALSE),
35 fVxMax(3.), fVyMax(3.), fVzMax(10.) {
38 for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
39 fHistNMult[iBin-1] = NULL;
40 fHistNPlusNMinus[iBin-1] = NULL;
43 // Define input and output slots here
44 // Input slot #0 works with a TChain
45 DefineInput(0, TChain::Class());
46 // Output slot #0 id reserved by the base class for AOD
47 // Output slot #1 writes into a TH1 container
48 DefineOutput(1, TList::Class());
51 //________________________________________________________________________
52 void AliEbyEFluctuationAnalysisTask::UserCreateOutputObjects() {
56 fOutputList = new TList();
57 fOutputList->SetOwner();
60 TString gCutName[4] = {"Total","Offline trigger",
62 fHistEventStats = new TH1F("fHistEventStats",
63 "Event statistics;;N_{events}",
65 for(Int_t i = 1; i <= 4; i++)
66 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
67 fOutputList->Add(fHistEventStats);
70 if(fAnalysisType == "ESD") {
71 fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
72 nCentralityBins,0.5,nCentralityBins+0.5);
73 fOutputList->Add(fHistCentrality);
76 for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
77 histName = "fHistNMult"; histName += "Centrality"; histName += iBin;
78 fHistNMult[iBin-1] = new TH1F(histName.Data(),
81 fOutputList->Add(fHistNMult[iBin-1]);
83 for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
84 histName = "fHistNPlusNMinus"; histName += "Centrality"; histName += iBin;
85 fHistNPlusNMinus[iBin-1] = new TH2F(histName.Data(),
87 2000,0.5,2000.5,2000,0.5,2000.5);
88 fOutputList->Add(fHistNPlusNMinus[iBin-1]);
91 else if(fAnalysisType == "MC") {
92 TString histName = "fHistNMultMC";
93 fHistNMultMC = new TH1F(histName.Data(),
96 fOutputList->Add(fHistNMultMC);
98 histName = "fHistNPlusNMinusMC";
99 fHistNPlusNMinusMC = new TH2F(histName.Data(),
101 3000,0.5,3000.5,3000,0.5,3000.5);
102 fOutputList->Add(fHistNPlusNMinusMC);
106 //________________________________________________________________________
107 void AliEbyEFluctuationAnalysisTask::UserExec(Option_t *) {
109 // Called for each event
111 Int_t nPlus = 0, nMinus = 0;
115 if(fAnalysisType == "ESD") {
116 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
118 printf("ERROR: fESD not available\n");
122 fHistEventStats->Fill(1); //all events
124 //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
126 fHistEventStats->Fill(2); //triggered + centrality
127 Printf("Event accepted");
130 AliCentrality *centrality = fESD->GetCentrality();
131 Int_t nCentrality = 0;
132 if(fCentralityBins20)
133 nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());
135 nCentrality = centrality->GetCentralityClass10(fCentralityEstimator.Data());
137 if((nCentrality < 0)||(nCentrality > 19)) return;
139 if(fAnalysisMode == "TPC") {
140 const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
142 if(vertex->GetNContributors() > 0) {
143 if(vertex->GetZRes() != 0) {
144 fHistEventStats->Fill(3); //events with a proper vertex
145 if(TMath::Abs(vertex->GetXv()) < fVxMax) {
146 if(TMath::Abs(vertex->GetYv()) < fVyMax) {
147 if(TMath::Abs(vertex->GetZv()) < fVzMax) {
148 fHistEventStats->Fill(4); //analyzed events
150 //Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",
151 //centrality->GetCentralityPercentile(fCentralityEstimator.Data()),
152 //nCentrality,fESD->GetNumberOfTracks());
154 Int_t nAcceptedTracks = 0;
155 TObjArray *gTrackArray = 0;
157 gTrackArray = fESDtrackCuts->GetAcceptedTracks(fESD,kTRUE);
159 nAcceptedTracks = gTrackArray->GetEntries();
161 AliESDtrack* track = 0;
162 // Track loop to fill a pT spectrum
163 for (Int_t iTracks = 0; iTracks < nAcceptedTracks; iTracks++) {
164 track = dynamic_cast<AliESDtrack *>(gTrackArray->At(iTracks));
166 printf("ERROR: Could not receive track %d\n", iTracks);
170 Short_t gCharge = track->Charge();
171 if(gCharge > 0) nPlus += 1;
172 if(gCharge < 0) nMinus += 1;
174 }//TObjArray valid object
175 //if((nCentrality >= 1)&&(nCentrality <= 20)) {
177 fHistCentrality->Fill(nCentrality);
178 fHistNPlusNMinus[nCentrality-1]->Fill(nPlus,nMinus);
179 fHistNMult[nCentrality-1]->Fill(nPlus+nMinus);
185 }//number of contributors
188 //}//physics selection
189 }//ESD analysis level
192 if(fAnalysisType == "MC") {
193 AliMCEvent* mcEvent = MCEvent();
195 Printf("ERROR: Could not retrieve MC event");
198 AliStack* stack = mcEvent->Stack();
200 Printf("ERROR: Could not retrieve MC stack");
204 fHistEventStats->Fill(1);
205 fHistEventStats->Fill(2);
206 fHistEventStats->Fill(3);
207 fHistEventStats->Fill(4); //analyzed events
208 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfTracks(); iTracks++) {
209 AliVParticle* track = mcEvent->GetTrack(iTracks);
211 Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
215 if(!stack->IsPhysicalPrimary(iTracks)) continue;
216 if((track->Pt() < 0.3) && (track->Pt() > 1.5)) continue;
217 if(TMath::Abs(track->Eta()) > 0.8) continue;
219 Short_t gCharge = track->Charge();
220 if(gCharge > 0) nPlus += 1;
221 if(gCharge < 0) nMinus += 1;
223 fHistNPlusNMinusMC->Fill(nPlus,nMinus);
224 fHistNMultMC->Fill(nPlus+nMinus);
229 PostData(1, fOutputList);
232 //________________________________________________________________________
233 void AliEbyEFluctuationAnalysisTask::Terminate(Option_t *) {
234 // Draw result to the screen
235 // Called once at the end of the query
237 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
239 printf("ERROR: Output list not available\n");