changes for Vertex and Tracks classes
[u/mrichter/AliRoot.git] / PWGCF / EBYE / Fluctuations / AliEbyEFluctuationAnalysisTask.cxx
CommitLineData
3eacba55 1#include "TChain.h"
2#include "TTree.h"
3#include "TH1F.h"
4#include "TH2F.h"
5#include "TCanvas.h"
d9649a65 6#include "TParticle.h"
8d9f1931 7#include "TObjArray.h"
8
d9649a65 9#include "AliAnalysisTask.h"
3eacba55 10#include "AliAnalysisManager.h"
d9649a65 11
3eacba55 12#include "AliStack.h"
d9649a65 13#include "AliMCEvent.h"
3eacba55 14#include "AliESDEvent.h"
15#include "AliESDInputHandler.h"
d9649a65 16#include "AliESDtrackCuts.h"
3eacba55 17#include "AliCentrality.h"
3eacba55 18
d9649a65 19#include "AliEbyEFluctuationAnalysisTask.h"
3eacba55 20
d9649a65 21// Event by event charge fluctuation analysis
22// Authors: Satyajit Jena and Panos Cristakoglou
23//
3eacba55 24
d9649a65 25ClassImp(AliEbyEFluctuationAnalysisTask)
3eacba55 26
27//________________________________________________________________________
28AliEbyEFluctuationAnalysisTask::AliEbyEFluctuationAnalysisTask(const char *name)
d9649a65 29 : AliAnalysisTaskSE(name), fESD(0), fOutputList(0),
30 fHistEventStats(0), fHistCentrality(0),
31 fHistNMultMC(0), fHistNPlusNMinusMC(0),
32 fESDtrackCuts(0),
33 fAnalysisType("ESD"), fAnalysisMode("TPC"),
34 fCentralityEstimator("V0M"), fCentralityBins20(kFALSE),
bfd18345 35 fVxMax(3.), fVyMax(3.), fVzMax(10.) {
d9649a65 36 // Constructor
37
38 for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
39 fHistNMult[iBin-1] = NULL;
40 fHistNPlusNMinus[iBin-1] = NULL;
3eacba55 41 }
d9649a65 42
43 // Define input and output slots here
44 // Input slot #0 works with a TChain
bfd18345 45 DefineInput(0, TChain::Class());
d9649a65 46 // Output slot #0 id reserved by the base class for AOD
47 // Output slot #1 writes into a TH1 container
bfd18345 48 DefineOutput(1, TList::Class());
3eacba55 49}
50
51//________________________________________________________________________
52void AliEbyEFluctuationAnalysisTask::UserCreateOutputObjects() {
d9649a65 53 // Create histograms
54 // Called once
55
bfd18345 56 fOutputList = new TList();
57 fOutputList->SetOwner();
d9649a65 58
59 //Event stats.
60 TString gCutName[4] = {"Total","Offline trigger",
61 "Vertex","Analyzed"};
62 fHistEventStats = new TH1F("fHistEventStats",
63 "Event statistics;;N_{events}",
64 4,0.5,4.5);
3eacba55 65 for(Int_t i = 1; i <= 4; i++)
66 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
d9649a65 67 fOutputList->Add(fHistEventStats);
68
69 //ESD analysis
70 if(fAnalysisType == "ESD") {
71 fHistCentrality = new TH1F("fHistCentrality",";Centrality bin;Events",
72 nCentralityBins,0.5,nCentralityBins+0.5);
73 fOutputList->Add(fHistCentrality);
74
75 TString histName;
76 for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
77 histName = "fHistNMult"; histName += "Centrality"; histName += iBin;
78 fHistNMult[iBin-1] = new TH1F(histName.Data(),
79 ";N_{mult.}",
80 500,0,3000);
81 fOutputList->Add(fHistNMult[iBin-1]);
82 }
83 for(Int_t iBin = 1; iBin <= nCentralityBins; iBin++) {
84 histName = "fHistNPlusNMinus"; histName += "Centrality"; histName += iBin;
85 fHistNPlusNMinus[iBin-1] = new TH2F(histName.Data(),
86 ";N_{+};N_{-}",
87 2000,0.5,2000.5,2000,0.5,2000.5);
88 fOutputList->Add(fHistNPlusNMinus[iBin-1]);
89 }
90 }//ESD analysis
91 else if(fAnalysisType == "MC") {
92 TString histName = "fHistNMultMC";
93 fHistNMultMC = new TH1F(histName.Data(),
94 ";N_{mult.}",
95 600,0,6000);
96 fOutputList->Add(fHistNMultMC);
97
98 histName = "fHistNPlusNMinusMC";
99 fHistNPlusNMinusMC = new TH2F(histName.Data(),
3eacba55 100 ";N_{+};N_{-}",
d9649a65 101 3000,0.5,3000.5,3000,0.5,3000.5);
102 fOutputList->Add(fHistNPlusNMinusMC);
103 }//MC analysis
3eacba55 104}
105
106//________________________________________________________________________
107void AliEbyEFluctuationAnalysisTask::UserExec(Option_t *) {
d9649a65 108 // Main loop
109 // Called for each event
3eacba55 110
d9649a65 111 Int_t nPlus = 0, nMinus = 0;
112
113 // Post output data.
3eacba55 114 //ESD analysis
d9649a65 115 if(fAnalysisType == "ESD") {
116 fESD = dynamic_cast<AliESDEvent*>(InputEvent());
117 if (!fESD) {
118 printf("ERROR: fESD not available\n");
119 return;
120 }
121
122 fHistEventStats->Fill(1); //all events
123
124 //Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
125 //if(isSelected) {
126 fHistEventStats->Fill(2); //triggered + centrality
127 Printf("Event accepted");
128
129 //Centrality stuff
130 AliCentrality *centrality = fESD->GetCentrality();
131 Int_t nCentrality = 0;
132 if(fCentralityBins20)
133 nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());
134 else
135 nCentrality = centrality->GetCentralityClass10(fCentralityEstimator.Data());
136
137 if((nCentrality < 0)||(nCentrality > 19)) return;
138
699d900c 139 if(fAnalysisMode == "TPC") {
d9649a65 140 const AliESDVertex *vertex = fESD->GetPrimaryVertexTPC();
141 if(vertex) {
142 if(vertex->GetNContributors() > 0) {
143 if(vertex->GetZRes() != 0) {
144 fHistEventStats->Fill(3); //events with a proper vertex
e690d4d0 145 if(TMath::Abs(vertex->GetX()) < fVxMax) {
146 if(TMath::Abs(vertex->GetY()) < fVyMax) {
147 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
d9649a65 148 fHistEventStats->Fill(4); //analyzed events
149
150 //Printf("Centrality percentile: %lf - Centrality: %d - Total tracks: %d",
151 //centrality->GetCentralityPercentile(fCentralityEstimator.Data()),
152 //nCentrality,fESD->GetNumberOfTracks());
153
8d9f1931 154 Int_t nAcceptedTracks = 0;
155 TObjArray *gTrackArray = 0;
156 if(fESDtrackCuts)
157 gTrackArray = fESDtrackCuts->GetAcceptedTracks(fESD,kTRUE);
158 if(gTrackArray) {
159 nAcceptedTracks = gTrackArray->GetEntries();
160
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));
165 if (!track) {
166 printf("ERROR: Could not receive track %d\n", iTracks);
167 continue;
168 }
d9649a65 169
8d9f1931 170 Short_t gCharge = track->Charge();
d9649a65 171 if(gCharge > 0) nPlus += 1;
172 if(gCharge < 0) nMinus += 1;
8d9f1931 173 }//track loop
174 }//TObjArray valid object
d9649a65 175 //if((nCentrality >= 1)&&(nCentrality <= 20)) {
8d9f1931 176
d9649a65 177 fHistCentrality->Fill(nCentrality);
178 fHistNPlusNMinus[nCentrality-1]->Fill(nPlus,nMinus);
179 fHistNMult[nCentrality-1]->Fill(nPlus+nMinus);
180 //}
181 }//Vz cut
182 }//Vy cut
183 }//Vx cut
184 }//Vz resolution
185 }//number of contributors
186 }//valid vertex
187 }//TPC analysis mode
188 //}//physics selection
189 }//ESD analysis level
8d9f1931 190
d9649a65 191 //MC analysis
192 if(fAnalysisType == "MC") {
193 AliMCEvent* mcEvent = MCEvent();
194 if (!mcEvent) {
195 Printf("ERROR: Could not retrieve MC event");
196 return;
197 }
198 AliStack* stack = mcEvent->Stack();
199 if (!stack) {
200 Printf("ERROR: Could not retrieve MC stack");
3eacba55 201 return;
202 }
3eacba55 203
d9649a65 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);
3eacba55 210 if (!track) {
d9649a65 211 Printf("ERROR: Could not receive track %d (mc loop)", iTracks);
3eacba55 212 continue;
213 }
214
d9649a65 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;
218
219 Short_t gCharge = track->Charge();
220 if(gCharge > 0) nPlus += 1;
221 if(gCharge < 0) nMinus += 1;
222 }//particle loop
223 fHistNPlusNMinusMC->Fill(nPlus,nMinus);
224 fHistNMultMC->Fill(nPlus+nMinus);
225
226 }//MC analysis level
3eacba55 227
3eacba55 228
bfd18345 229 PostData(1, fOutputList);
d9649a65 230}
3eacba55 231
232//________________________________________________________________________
233void AliEbyEFluctuationAnalysisTask::Terminate(Option_t *) {
d9649a65 234 // Draw result to the screen
235 // Called once at the end of the query
3eacba55 236
d9649a65 237 fOutputList = dynamic_cast<TList*> (GetOutputData(1));
238 if (!fOutputList) {
239 printf("ERROR: Output list not available\n");
240 return;
3eacba55 241 }
3eacba55 242}