4 #include "TLorentzVector.h"
\r
5 #include "TGraphErrors.h"
\r
9 #include "AliAnalysisTaskSE.h"
\r
10 #include "AliAnalysisManager.h"
\r
12 #include "AliESDVertex.h"
\r
13 #include "AliESDEvent.h"
\r
14 #include "AliESDInputHandler.h"
\r
15 #include "AliAODEvent.h"
\r
16 #include "AliAODTrack.h"
\r
17 #include "AliAODInputHandler.h"
\r
18 #include "AliMCEventHandler.h"
\r
19 #include "AliMCEvent.h"
\r
20 #include "AliStack.h"
\r
21 #include "AliESDtrackCuts.h"
\r
23 #include "AliAnalysisTaskBF.h"
\r
24 #include "AliBalance.h"
\r
27 // Analysis task for the BF code
\r
28 // Authors: Panos Cristakoglou@cern.ch
\r
30 ClassImp(AliAnalysisTaskBF)
\r
32 //________________________________________________________________________
\r
33 AliAnalysisTaskBF::AliAnalysisTaskBF(const char *name)
\r
34 : AliAnalysisTaskSE(name),
\r
52 fCentralityEstimator("VOM"),
\r
53 fCentralityPercentileMin(0.),
\r
54 fCentralityPercentileMax(5.),
\r
55 fUseOfflineTrigger(kFALSE),
\r
64 for(Int_t i = 0; i < NUMBER_OF_ANALYSES; i++){
\r
65 for (Int_t j = 0; j < 3; j++){
\r
66 fHistBF[i][j] = NULL;
\r
70 // Define input and output slots here
\r
71 // Input slot #0 works with a TChain
\r
72 DefineInput(0, TChain::Class());
\r
73 // Output slot #0 writes into a TH1 container
\r
74 DefineOutput(1, AliBalance::Class());
\r
75 DefineOutput(2, TList::Class());
\r
76 DefineOutput(3, TList::Class());
\r
80 //________________________________________________________________________
\r
81 void AliAnalysisTaskBF::UserCreateOutputObjects() {
\r
82 // Create histograms
\r
85 fBalance = new AliBalance();
\r
86 fBalance->SetAnalysisLevel("ESD");
\r
87 fBalance->SetNumberOfBins(-1,9);
\r
88 fBalance->SetInterval(-1,0.,0.9);
\r
92 fList = new TList();
\r
93 fList->SetName("listQA");
\r
96 //Balance Function list
\r
97 fListBF = new TList();
\r
98 fListBF->SetName("listBF");
\r
99 fListBF->SetOwner();
\r
102 TString gCutName[4] = {"Total","Offline trigger",
\r
103 "Vertex","Analyzed"};
\r
104 fHistEventStats = new TH1F("fHistEventStats",
\r
105 "Event statistics;;N_{events}",
\r
107 for(Int_t i = 1; i <= 4; i++)
\r
108 fHistEventStats->GetXaxis()->SetBinLabel(i,gCutName[i-1].Data());
\r
109 fList->Add(fHistEventStats);
\r
111 fHistTrackStats = new TH1F("fHistTrackStats","Event statistics;TriggerBit;N_{events}",130,0,130);
\r
112 fList->Add(fHistTrackStats);
\r
114 // Vertex distributions
\r
115 fHistVx = new TH1F("fHistVx","Primary vertex distribution - x coordinate;V_{x} (cm);Entries",100,-0.5,0.5);
\r
116 fList->Add(fHistVx);
\r
117 fHistVy = new TH1F("fHistVy","Primary vertex distribution - y coordinate;V_{y} (cm);Entries",100,-0.5,0.5);
\r
118 fList->Add(fHistVy);
\r
119 fHistVz = new TH1F("fHistVz","Primary vertex distribution - z coordinate;V_{z} (cm);Entries",100,-20.,20.);
\r
120 fList->Add(fHistVz);
\r
123 fHistClus = new TH2F("fHistClus","# Cluster (TPC vs. ITS)",10,0,10,200,0,200);
\r
124 fList->Add(fHistClus);
\r
125 fHistChi2 = new TH1F("fHistChi2","Chi2/NDF distribution",200,0,10);
\r
126 fList->Add(fHistChi2);
\r
127 fHistDCA = new TH2F("fHistDCA","DCA (xy vs. z)",400,-5,5,400,-5,5);
\r
128 fList->Add(fHistDCA);
\r
129 fHistPt = new TH1F("fHistPt","p_{T} distribution",200,0,10);
\r
130 fList->Add(fHistPt);
\r
131 fHistEta = new TH1F("fHistEta","#eta distribution",200,-2,2);
\r
132 fList->Add(fHistEta);
\r
133 fHistPhi = new TH1F("fHistPhi","#phi distribution",200,-20,380);
\r
134 fList->Add(fHistPhi);
\r
135 fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",500, 0, 20000, 500, 0, 20000);
\r
136 fList->Add(fHistV0M);
\r
139 if(fESDtrackCuts) fList->Add(fESDtrackCuts);
\r
142 //balance function histograms
\r
144 for(Int_t i = 0; i < NUMBER_OF_ANALYSES; i++){
\r
145 for (Int_t j = 0; j < 3; j++){
\r
149 fHistBF[i][j] = new TH1F(hname,hname,fBalance->GetNumberOfBins(i),fBalance->GetP2Start(i),fBalance->GetP2Stop(i));
\r
150 fListBF->Add(fHistBF[i][j]);
\r
153 fHistN = new TH1F("fN","fN",2,0,1);
\r
154 fListBF->Add(fHistN);
\r
156 // Post output data.
\r
157 //PostData(1, fBalance);
\r
158 PostData(2, fList);
\r
159 PostData(3, fListBF);
\r
163 //________________________________________________________________________
\r
164 void AliAnalysisTaskBF::UserExec(Option_t *) {
\r
166 // Called for each event
\r
167 TString gAnalysisLevel = fBalance->GetAnalysisLevel();
\r
169 TObjArray *array = new TObjArray();
\r
172 if(gAnalysisLevel == "ESD") {
\r
173 AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE
\r
175 Printf("ERROR: gESD not available");
\r
179 fHistEventStats->Fill(1); //all events
\r
180 Bool_t isSelected = kTRUE;
\r
181 if(fUseOfflineTrigger)
\r
182 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
\r
184 fHistEventStats->Fill(2); //triggered events
\r
187 AliCentrality *centrality = gESD->GetCentrality();
\r
188 //Int_t nCentrality = 0;
\r
189 //nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());
\r
190 //cout<<nCentrality<<" "<<centrality->IsEventInCentralityClass(fCentralityPercentileMin,fCentralityPercentileMax,fCentralityEstimator.Data())<<endl;
\r
192 // take only events inside centrality class
\r
193 if(centrality->IsEventInCentralityClass(fCentralityPercentileMin,
\r
194 fCentralityPercentileMax,
\r
195 fCentralityEstimator.Data())){
\r
197 // centrality QA (V0M)
\r
198 fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());
\r
200 const AliESDVertex *vertex = gESD->GetPrimaryVertex();
\r
202 if(vertex->GetNContributors() > 0) {
\r
203 if(vertex->GetZRes() != 0) {
\r
204 fHistEventStats->Fill(3); //events with a proper vertex
\r
205 if(TMath::Abs(vertex->GetXv()) < fVxMax) {
\r
206 if(TMath::Abs(vertex->GetYv()) < fVyMax) {
\r
207 if(TMath::Abs(vertex->GetZv()) < fVzMax) {
\r
208 fHistEventStats->Fill(4); //analayzed events
\r
209 fHistVx->Fill(vertex->GetXv());
\r
210 fHistVy->Fill(vertex->GetYv());
\r
211 fHistVz->Fill(vertex->GetZv());
\r
213 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());
\r
214 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {
\r
215 AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));
\r
217 Printf("ERROR: Could not receive track %d", iTracks);
\r
221 // take only TPC only tracks (HOW IS THIS DONE IN ESDs???)
\r
222 //if(!track->IsTPCOnly()) continue;
\r
226 if(!fESDtrackCuts->AcceptTrack(track)) continue;
\r
228 // fill QA histograms
\r
231 track->GetImpactParameters(b,bCov);
\r
232 if (bCov[0]<=0 || bCov[2]<=0) {
\r
233 AliDebug(1, "Estimated b resolution lower or equal zero!");
\r
234 bCov[0]=0; bCov[2]=0;
\r
237 fHistClus->Fill(track->GetITSclusters(0),track->GetTPCclusters(0));
\r
238 fHistDCA->Fill(b[1],b[0]);
\r
239 fHistPt->Fill(track->Pt());
\r
240 fHistEta->Fill(track->Eta());
\r
241 fHistPhi->Fill(track->Phi()*TMath::RadToDeg());
\r
250 }//proper vertex resolution
\r
251 }//proper number of contributors
\r
252 }//vertex object valid
\r
254 }//triggered event
\r
258 //AOD analysis (vertex and track cuts also here!!!!)
\r
259 else if(gAnalysisLevel == "AOD") {
\r
260 AliAODEvent* gAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // from TaskSE
\r
262 Printf("ERROR: gAOD not available");
\r
266 fHistEventStats->Fill(1); //all events
\r
267 Bool_t isSelected = kTRUE;
\r
268 if(fUseOfflineTrigger)
\r
269 isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
\r
271 fHistEventStats->Fill(2); //triggered events
\r
274 //Centrality stuff (centrality in AOD header)
\r
275 AliAODHeader *aodHeader = gAOD->GetHeader();
\r
276 Float_t fCentrality = aodHeader->GetCentralityP()->GetCentralityPercentile("V0M");
\r
278 // take only events inside centrality class
\r
279 if(fCentrality > fCentralityPercentileMin && fCentrality < fCentralityPercentileMax){
\r
281 // centrality QA (V0M)
\r
282 fHistV0M->Fill(gAOD->GetVZEROData()->GetMTotV0A(), gAOD->GetVZEROData()->GetMTotV0C());
\r
285 //cout<<gAOD->GetNumberOfVertices()<<endl;
\r
286 const AliAODVertex *vertex = gAOD->GetPrimaryVertex();
\r
287 //const AliAODVertex *vertex1 = gAOD->GetPrimaryVertexSPD();
\r
288 // vertex->Print();
\r
289 // vertex1->Print();
\r
292 Double32_t fCov[6];
\r
293 Double32_t fPosition[3];
\r
294 Double32_t fChi2 = vertex->GetChi2();
\r
295 Int_t nCont = vertex->GetNContributors();
\r
296 vertex->GetXYZ(fPosition);
\r
297 vertex->GetCovarianceMatrix(fCov);
\r
299 const AliESDVertex *esdVertex = new AliESDVertex(fPosition, fCov, fChi2, nCont);
\r
301 if(vertex->GetNContributors() > 0) {
\r
302 //if(vertex->GetZRes() != 0) {
\r
303 fHistEventStats->Fill(3); //events with a proper vertex
\r
304 if(TMath::Abs(vertex->GetX()) < fVxMax) {
\r
305 if(TMath::Abs(vertex->GetY()) < fVyMax) {
\r
306 if(TMath::Abs(vertex->GetZ()) < fVzMax) {
\r
307 fHistEventStats->Fill(4); //analayzed events
\r
308 fHistVx->Fill(vertex->GetX());
\r
309 fHistVy->Fill(vertex->GetY());
\r
310 fHistVz->Fill(vertex->GetZ());
\r
312 //Printf("There are %d tracks in this event", gAOD->GetNumberOfTracks());
\r
313 for (Int_t iTracks = 0; iTracks < gAOD->GetNumberOfTracks(); iTracks++) {
\r
314 AliAODTrack* aodTrack = dynamic_cast<AliAODTrack *>(gAOD->GetTrack(iTracks));
\r
316 Printf("ERROR: Could not receive track %d", iTracks);
\r
322 // For ESD Filter Information: ANALYSIS/macros/AddTaskESDfilter.C
\r
323 // take only TPC only tracks
\r
324 fHistTrackStats->Fill(aodTrack->GetFilterMap());
\r
325 if(!aodTrack->TestFilterBit(128)) continue;
\r
330 aodTrack->GetPosition(p);
\r
331 //AliExternalTrackParam exParam;
\r
332 AliESDtrack esdTrack( aodTrack );
\r
333 esdTrack.RelateToVertex(esdVertex,gAOD->GetMagneticField(),100);
\r
334 //esdTrack.RelateToVertexTPC(esdVertex,gAOD->GetMagneticField(),100,&exParam);
\r
335 //AliExternalTrackParam* outer = (AliExternalTrackParam*)aodTrack->GetOuterParam();
\r
336 //if(outer) outer->PropagateToDCA(vertex,gAOD->GetMagneticField(),1000,b,bCov);
\r
338 //cout<<" ___ "<<b[0]<<" "<<b[1]<<endl;
\r
341 //only for QA --> Remove?
\r
342 esdTrack.GetImpactParameters(b,bCov);
\r
343 if (bCov[0]<=0 || bCov[2]<=0) {
\r
344 AliDebug(1, "Estimated b resolution lower or equal zero!");
\r
345 bCov[0]=0; bCov[2]=0;
\r
348 Float_t pt = aodTrack->Pt();
\r
349 Float_t eta = aodTrack->Eta();
\r
351 // Kinematics cuts from ESD track cuts
\r
352 if( pt < fPtMin || pt > fPtMax) continue;
\r
353 if( eta < fEtaMin || eta > fEtaMax) continue;
\r
355 // fill QA histograms
\r
356 fHistClus->Fill(aodTrack->GetITSNcls(),aodTrack->GetTPCNcls());
\r
357 fHistDCA->Fill(b[1],b[0]);
\r
358 //fHistDCA->Fill(p[1],p[0]);
\r
359 fHistChi2->Fill(aodTrack->Chi2perNDF());
\r
361 fHistEta->Fill(eta);
\r
362 fHistPhi->Fill(aodTrack->Phi()*TMath::RadToDeg());
\r
365 array->Add(aodTrack);
\r
371 //}//proper vertex resolution
\r
372 }//proper number of contributors
\r
373 }//vertex object valid
\r
375 }//triggered event
\r
379 else if(gAnalysisLevel == "MC") {
\r
381 AliMCEvent* mcEvent = MCEvent();
\r
383 Printf("ERROR: mcEvent not available");
\r
387 Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());
\r
388 for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {
\r
389 AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));
\r
391 Printf("ERROR: Could not receive particle %d", iTracks);
\r
398 fBalance->CalculateBalance(array);
\r
404 //________________________________________________________________________
\r
405 void AliAnalysisTaskBF::FinishTaskOutput(){
\r
406 //Printf("END BF");
\r
408 //fBalance = dynamic_cast<AliBalance*> (GetOutputData(1));
\r
410 Printf("ERROR: fBalance not available");
\r
414 for(Int_t a = 0; a < NUMBER_OF_ANALYSES; a++){
\r
415 for(Int_t iBin = 1; iBin <= fBalance->GetNumberOfBins(a); iBin++){
\r
417 //Printf("%d %d -> %f",a,iBin,fBalance->GetNnn(a,iBin-1));
\r
419 fHistBF[a][0]->SetBinContent(iBin,fBalance->GetNnn(a,iBin-1));
\r
420 fHistBF[a][1]->SetBinContent(iBin,fBalance->GetNpn(a,iBin-1));
\r
421 fHistBF[a][2]->SetBinContent(iBin,fBalance->GetNpp(a,iBin-1));
\r
424 fHistN->SetBinContent(1,fBalance->GetNn());
\r
425 fHistN->SetBinContent(2,fBalance->GetNp());
\r
429 //________________________________________________________________________
\r
430 void AliAnalysisTaskBF::Terminate(Option_t *) {
\r
431 // Draw result to the screen
\r
432 // Called once at the end of the query
\r
434 // not implemented ...
\r