]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/EVCHAR/AliAnalysisTaskCentralityTreeMaker.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
follow the old and nasty naming for consistency: FilterEvents_Trees.root
[u/mrichter/AliRoot.git] / PWGPP / EVCHAR / AliAnalysisTaskCentralityTreeMaker.cxx
CommitLineData
c5297977 1/**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
16/////////////////////////////////////////////////////////////
17// //
18// Class to analyze centrality measurements //
19// //
20/////////////////////////////////////////////////////////////
21
22#include <TTree.h>
23#include <TList.h>
24#include <TH2F.h>
25#include <TFile.h>
26#include <TString.h>
27#include <TCanvas.h>
28
29#include "AliAnalysisManager.h"
30#include "AliVEvent.h"
31#include "AliESD.h"
32#include "AliESDEvent.h"
33#include "AliESDHeader.h"
34#include "AliESDInputHandler.h"
35#include "AliESDZDC.h"
36#include "AliESDFMD.h"
37#include "AliESDVZERO.h"
38#include "AliMultiplicity.h"
39#include "AliAODHandler.h"
40#include "AliAODEvent.h"
41#include "AliAODVertex.h"
42#include "AliAODMCHeader.h"
43#include "AliMCEvent.h"
44#include "AliMCEventHandler.h"
45#include "AliMCParticle.h"
46#include "AliStack.h"
47#include "AliHeader.h"
48#include "AliAODMCParticle.h"
49#include "AliAnalysisTaskSE.h"
50#include "AliGenEventHeader.h"
51#include "AliGenHijingEventHeader.h"
52#include "AliPhysicsSelectionTask.h"
53#include "AliPhysicsSelection.h"
54#include "AliBackgroundSelection.h"
55#include "AliAnalysisTaskCentralityTreeMaker.h"
56
57ClassImp(AliAnalysisTaskCentralityTreeMaker)
58
59
60//________________________________________________________________________
61AliAnalysisTaskCentralityTreeMaker::AliAnalysisTaskCentralityTreeMaker():
62 AliAnalysisTaskSE(),
63 fDebug(0),
64 fAnalysisInput("ESD"),
65 fIsMCInput(kFALSE),
66 fOutput(0x0),
67 fEZDCvsEZEM(0x0),
68 fEZDCvsNtracklets(0x0),
69 fTreeFilling(kTRUE),
70 fCentralityTree(0x0),
71 fNev(0),
72 fBeamEnergy(0),
73 fNmyTracks_gen(0),
74 fxVertex(0),
75 fyVertex(0),
76 fzVertex(0),
77 fVertexer3d(kFALSE),
78 fbMC(0.),
79 fNpartTargMC(0),
80 fNpartProjMC(0),
81 fNNColl(0),
82 fNNwColl(0),
83 fNwNColl(0),
84 fNwNwColl(0),
85 fNTracklets(0),
86 fNSingleClusters(0),
87 fbZDC(0),
88 fNpartZDC(0),
89 fbZDCA(0),
90 fNpartZDCA(0),
91 fbZDCC(0),
92 fNpartZDCC(0),
93 fESDFlag(0),
94 fZNCEnergy(0),
95 fZPCEnergy(0),
96 fZNAEnergy(0),
97 fZPAEnergy(0),
98 fZEM1Energy(0),
99 fZEM2Energy(0),
100 fNTracks(0),
101 fNPmdTracks(0),
102 fMultV0A(0),
103 fMultV0C(0),
104 fMultFMDA(0),
105 fMultFMDC(0)
106{
107 // Default constructor
108}
109
110//________________________________________________________________________
111AliAnalysisTaskCentralityTreeMaker::AliAnalysisTaskCentralityTreeMaker(const char *name):
112 AliAnalysisTaskSE(name),
113 fDebug(0),
114 fAnalysisInput("ESD"),
115 fIsMCInput(kFALSE),
116 fOutput(0x0),
117 fEZDCvsEZEM(0x0),
118 fEZDCvsNtracklets(0x0),
119 fTreeFilling(kTRUE),
120 fCentralityTree(0x0),
121 fNev(0),
122 fBeamEnergy(0),
123 fNmyTracks_gen(0),
124 fxVertex(0),
125 fyVertex(0),
126 fzVertex(0),
127 fVertexer3d(kFALSE),
128 fbMC(0.),
129 fNpartTargMC(0),
130 fNpartProjMC(0),
131 fNNColl(0),
132 fNNwColl(0),
133 fNwNColl(0),
134 fNwNwColl(0),
135 fNTracklets(0),
136 fNSingleClusters(0),
137 fbZDC(0),
138 fNpartZDC(0),
139 fbZDCA(0),
140 fNpartZDCA(0),
141 fbZDCC(0),
142 fNpartZDCC(0),
143 fESDFlag(0),
144 fZNCEnergy(0),
145 fZPCEnergy(0),
146 fZNAEnergy(0),
147 fZPAEnergy(0),
148 fZEM1Energy(0),
149 fZEM2Energy(0),
150 fNTracks(0),
151 fNPmdTracks(0),
152 fMultV0A(0),
153 fMultV0C(0),
154 fMultFMDA(0),
155 fMultFMDC(0)
156{
157 // Default constructor
158
159 // Output slot #1 writes into a TList container
160 DefineOutput(1, TList::Class());
161 // Output slot #2 writes into a TTree container
162 if(fTreeFilling) DefineOutput(2, TTree::Class());
163
164}
165
166//________________________________________________________________________
167AliAnalysisTaskCentralityTreeMaker& AliAnalysisTaskCentralityTreeMaker::operator=(const AliAnalysisTaskCentralityTreeMaker& c)
168{
169 //
170 // Assignment operator
171 //
172 if (this!=&c) {
173 AliAnalysisTaskSE::operator=(c);
174 }
175 return *this;
176}
177
178//________________________________________________________________________
179AliAnalysisTaskCentralityTreeMaker::AliAnalysisTaskCentralityTreeMaker(const AliAnalysisTaskCentralityTreeMaker& ana):
180 AliAnalysisTaskSE(ana),
181 fDebug(ana.fDebug),
182 fAnalysisInput(ana.fDebug),
183 fIsMCInput(ana.fIsMCInput),
184 fOutput(ana.fOutput),
185 fEZDCvsEZEM(ana.fEZDCvsEZEM),
186 fEZDCvsNtracklets(ana.fEZDCvsNtracklets),
187 fTreeFilling(ana.fTreeFilling),
188 fCentralityTree(ana.fCentralityTree),
189 fNev(ana.fDebug),
190 fBeamEnergy(ana.fBeamEnergy),
191 fNmyTracks_gen(ana.fNmyTracks_gen),
192 fxVertex(ana.fxVertex),
193 fyVertex(ana.fyVertex),
194 fzVertex(ana.fzVertex),
195 fVertexer3d(ana.fVertexer3d),
196 fbMC(ana.fbMC),
197 fNpartTargMC(ana.fNpartTargMC),
198 fNpartProjMC(ana.fNpartProjMC),
199 fNNColl(ana.fNNColl),
200 fNNwColl(ana.fNNwColl),
201 fNwNColl(ana.fNwNColl),
202 fNwNwColl(ana.fNwNwColl),
203 fNTracklets(ana.fNTracklets),
204 fNSingleClusters(ana.fNSingleClusters),
205 fbZDC(ana.fbZDC),
206 fNpartZDC(ana.fNpartZDC),
207 fbZDCA(ana.fbZDCA),
208 fNpartZDCA(ana.fNpartZDCA),
209 fbZDCC(ana.fbZDCC),
210 fNpartZDCC(ana.fNpartZDCC),
211 fESDFlag(ana.fESDFlag),
212 fZNCEnergy(ana.fZNCEnergy),
213 fZPCEnergy(ana.fZPCEnergy),
214 fZNAEnergy(ana.fZNAEnergy),
215 fZPAEnergy(ana.fZPAEnergy),
216 fZEM1Energy(ana.fZEM1Energy),
217 fZEM2Energy(ana.fZEM2Energy),
218 fNTracks(ana.fNTracks),
219 fNPmdTracks(ana.fNPmdTracks),
220 fMultV0A(ana.fMultV0A),
221 fMultV0C(ana.fMultV0C),
222 fMultFMDA(ana.fMultFMDA),
223 fMultFMDC(ana.fMultFMDC)
224{
225 //
226 // Copy Constructor
227 //
228}
229
230//________________________________________________________________________
231 AliAnalysisTaskCentralityTreeMaker::~AliAnalysisTaskCentralityTreeMaker()
232 {
233 // Destructor
234 if(fOutput){
235 delete fOutput; fOutput=0;
236 }
237 if(fCentralityTree){
238 delete fCentralityTree; fCentralityTree=0;
239 }
240
241 }
242
243//________________________________________________________________________
244void AliAnalysisTaskCentralityTreeMaker::UserCreateOutputObjects()
245{
246
247 // Create the output containers
248 if(fDebug>1) printf("AnalysisTaskZDCpp::UserCreateOutputObjects() \n");
249
250 // Several histograms are more conveniently managed in a TList
251 fOutput = new TList();
252 fOutput->SetOwner();
253 fOutput->SetName("OutputHistos");
254
255 fEZDCvsEZEM = new TH2F("fEZDCvsEZEM", "E_{ZDC} vs. E_{ZEM}", 100,0.,600.,100,0.,5000.);
256 fEZDCvsEZEM->GetXaxis()->SetTitle("ZEM signal (GeV)");
257 fEZDCvsEZEM->GetYaxis()->SetTitle("ZDC signal (TeV)");
258 fOutput->Add(fEZDCvsEZEM);
259
260 fEZDCvsNtracklets = new TH2F("fEZDCvsNtracklets", "E_{ZDC} vs. N_{tracklets}", 100,0.,600.,100,0.,6000.);
261 fEZDCvsNtracklets->GetXaxis()->SetTitle("N_{tracklets}");
262 fEZDCvsNtracklets->GetYaxis()->SetTitle("ZDC signal (TeV)");
263 fOutput->Add(fEZDCvsNtracklets);
264
265 if(fTreeFilling){
266 OpenFile(2);
267 fCentralityTree = new TTree("fCentralityTree", "Centrality vs. multiplicity tree");
268 //
269 fCentralityTree->Branch("nev", &fNev,"nev/I");
270 fCentralityTree->Branch("beamEnergy", &fBeamEnergy,"beamEnergy/F");
271 fCentralityTree->Branch("nmyTracks_gen", &fNmyTracks_gen,"nmyTracks_gen/I");
272 fCentralityTree->Branch("trigClass",&fTrigClass,"trigClass/C");
273 fCentralityTree->Branch("xVertex", &fxVertex,"xVertex/D");
274 fCentralityTree->Branch("yVertex", &fyVertex,"yVertex/D");
275 fCentralityTree->Branch("zVertex", &fzVertex,"zVertex/D");
276 fCentralityTree->Branch("vertexer3d", &fVertexer3d,"vertexer3d/O");
277 fCentralityTree->Branch("bMC", &fbMC,"bMC/D");
278 fCentralityTree->Branch("npartTargMC", &fNpartTargMC,"npartTargMC/I");
279 fCentralityTree->Branch("npartProjMC", &fNpartProjMC,"npartProjMC/I");
280 fCentralityTree->Branch("NNColl", &fNNColl,"NNColl/I");
281 fCentralityTree->Branch("NwNColl", &fNwNColl,"NwNColl/I");
282 fCentralityTree->Branch("NNwColl", &fNNwColl,"NNwColl/I");
283 fCentralityTree->Branch("NwNwColl", &fNwNwColl,"NwNwColl/I");
284 fCentralityTree->Branch("nTracklets", &fNTracklets,"nTracklets/I");
285 fCentralityTree->Branch("nSingleClusters", &fNSingleClusters,"nSingleClusters/I");
286 fCentralityTree->Branch("nClusters", fNClusters,"nClusters[6]/I");
287 fCentralityTree->Branch("nChips", fNChips,"nChips[2]/I");
288 fCentralityTree->Branch("bZDC", &fbZDC,"bZDC/D");
289 fCentralityTree->Branch("npartZDC", &fNpartZDC,"npartZDC/I");
290 fCentralityTree->Branch("bZDCA", &fbZDCA,"bZDCA/D");
291 fCentralityTree->Branch("npartZDCA", &fNpartZDCA,"npartZDCA/I");
292 fCentralityTree->Branch("bZDCC", &fbZDCC,"bZDCC/D");
293 fCentralityTree->Branch("npartZDCC", &fNpartZDCC,"npartZDCC/I");
294 fCentralityTree->Branch("esdFlag", &fESDFlag,"esdFlag/i");
295 fCentralityTree->Branch("zncEnergy", &fZNCEnergy, "zncEnergy/F");
296 fCentralityTree->Branch("zpcEnergy", &fZPCEnergy, "zpcEnergy/F");
297 fCentralityTree->Branch("znaEnergy", &fZNAEnergy, "znaEnergy/F");
298 fCentralityTree->Branch("zpaEnergy", &fZPAEnergy, "zpaEnergy/F");
299 fCentralityTree->Branch("zem1Energy", &fZEM1Energy, "zem1Energy/F");
300 fCentralityTree->Branch("zem2Energy", &fZEM2Energy, "zem2Energy/F");
301 fCentralityTree->Branch("znctower", fZNCtower, "znctower[5]/F");
302 fCentralityTree->Branch("zpctower", fZPCtower, "zpctower[5]/F");
303 fCentralityTree->Branch("znatower", fZNAtower, "znatower[5]/F");
304 fCentralityTree->Branch("zpatower", fZPAtower, "zpatower[5]/F");
305 fCentralityTree->Branch("centrZNC", fCentrZNC, "centrZNC[2]/F");
306 fCentralityTree->Branch("centrZNA", fCentrZNA, "centrZNA[2]/F");
307 fCentralityTree->Branch("nTracks", &fNTracks,"nTracks/I");
308 fCentralityTree->Branch("nPmdTracks", &fNPmdTracks,"nPmdTracks/I");
309 fCentralityTree->Branch("multV0A", &fMultV0A,"multV0A/F");
310 fCentralityTree->Branch("multV0C", &fMultV0C,"multV0C/F");
311 fCentralityTree->Branch("multFMDA", &fMultFMDA,"multFMDA/F");
312 fCentralityTree->Branch("multFMDC", &fMultFMDC,"multFMDC/F");
313 }
314}
315
316//________________________________________________________________________
317void AliAnalysisTaskCentralityTreeMaker::UserExec(Option_t */*option*/)
318{
319 // Execute analysis for current event:
320 if(fDebug>1) printf(" **** AliAnalysisTaskCentralityTreeMaker::UserExec() \n");
321
322 if(fAnalysisInput.CompareTo("ESD")==0){
323
324 // AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
325 AliVEvent* event = InputEvent();
326 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
327
328 fNev++;
329
330 fNTracks = event->GetNumberOfTracks();
331 fNPmdTracks = esd->GetNumberOfPmdTracks();
332
333 AliESDVZERO* esdV0 = esd->GetVZEROData();
334 fMultV0A=esdV0->GetMTotV0A();
335 fMultV0C=esdV0->GetMTotV0C();
336
337 if(fIsMCInput){
338
339 AliMCEvent* mcEvent = MCEvent();
340 if (!mcEvent) {
341 printf(" Could not retrieve MC event!!!\n");
342 return;
343 }
344
345 fNmyTracks_gen = 0;
346 AliStack *stack = 0x0; // needed for MC studies
347 stack = MCEvent()->Stack();
348 for (Int_t iTrack = 0; iTrack < MCEvent()->GetNumberOfTracks(); iTrack++) {
349 //get properties of mc particle
350 AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(iTrack);
351 // Primaries only
352 if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;
353 //charged tracks only
354 if (mcP->Particle()->GetPDG()->Charge() == 0) continue;
355 //same cuts as on ESDtracks
356// if(TMath::Abs(mcP->Eta())>0.9)continue;
357// if(mcP->Pt()<0.2)continue;
358// if(mcP->Pt()>200)continue;
359
360 fNmyTracks_gen ++;
361 }
362
363 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
364 if(!genHeader){
365 printf(" Event generator header not available!!!\n");
366 return;
367 }
368
369 if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
370 fbMC = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
371 Int_t specNeutronProj = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsn();
372 Int_t specProtonProj = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsp();
373 Int_t specNeutronTarg = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsn();
374 Int_t specProtonTarg = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsp();
375 fNpartTargMC = 208.-(specNeutronTarg+specProtonTarg);
376 fNpartProjMC = 208.-(specNeutronProj+specProtonProj);
377 fNNColl = ((AliGenHijingEventHeader*) genHeader)->NN();
378 fNNwColl = ((AliGenHijingEventHeader*) genHeader)->NNw();
379 fNwNColl = ((AliGenHijingEventHeader*) genHeader)->NwN();
380 fNwNwColl = ((AliGenHijingEventHeader*) genHeader)->NwNw();
381 }
382
383 }
384
385 fBeamEnergy = esd->GetBeamEnergy();
386
387 // ***** Trigger selection
388 TString triggerClass = esd->GetFiredTriggerClasses();
389 sprintf(fTrigClass,"%s",triggerClass.Data());
390
391 const AliESDVertex *vertex = esd->GetPrimaryVertexSPD();
392 fxVertex = vertex->GetX();
393 fyVertex = vertex->GetY();
394 fzVertex = vertex->GetZ();
395 if(vertex->IsFromVertexer3D()) fVertexer3d = kTRUE;
396 else fVertexer3d = kFALSE;
397 Double_t vertex3[3];
398 vertex->GetXYZ(vertex3);
399
400 const AliMultiplicity *mult = esd->GetMultiplicity();
401 fNTracklets = mult->GetNumberOfTracklets();
402
403 for(Int_t ilay=0; ilay<6; ilay++){
404 fNClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
405 }
406 fNSingleClusters = mult->GetNumberOfSingleClusters();
407
408 for(Int_t ilay=0; ilay<2; ilay++){
409 fNChips[ilay] = mult->GetNumberOfFiredChips(ilay);
410 }
411
412
413 AliESDFMD *fmd = esd->GetFMDData();
414 //fmd->Print();
415 Float_t totalMultA = 0;
416 Float_t totalMultC = 0;
417 const Float_t fFMDLowCut = 0.4;
418
419 for(UShort_t det=1;det<=3;det++) {
420 Int_t nRings = (det==1 ? 1 : 2);
421 for (UShort_t ir = 0; ir < nRings; ir++) {
422 Char_t ring = (ir == 0 ? 'I' : 'O');
423 UShort_t nsec = (ir == 0 ? 20 : 40);
424 UShort_t nstr = (ir == 0 ? 512 : 256);
425 for(UShort_t sec =0; sec < nsec; sec++) {
426 for(UShort_t strip = 0; strip < nstr; strip++) {
427
428 Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
429 if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
430
431 Float_t nParticles=0;
432
433 if(FMDmult > fFMDLowCut) {
434 nParticles = 1.;
435 }
436
437 if (det<3) totalMultA = totalMultA + nParticles;
438 else totalMultC = totalMultC + nParticles;
439
440 }
441 }
442 }
443 }
444 fMultFMDA = totalMultA;
445 fMultFMDC = totalMultC;
446
447
448 AliESDZDC *esdZDC = esd->GetESDZDC();
449
450 fESDFlag = esdZDC->GetESDQuality();
451
452 fZNCEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
453 fZPCEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
454 fZNAEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
455 fZPAEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
456 fZEM1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
457 fZEM2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
458 //
459 //printf(" E_ZEM = %1.2f GeV, E_ZDC = %1.2f TeV\n",fZEM1Energy+fZEM2Energy, (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.);
460 fEZDCvsEZEM->Fill(fZEM1Energy+fZEM2Energy, (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.);
461 fEZDCvsNtracklets->Fill(fNTracklets, (fZNCEnergy+fZPCEnergy+fZNAEnergy+fZPAEnergy)/1000.);
462 PostData(1, fOutput);
463
464 fbZDC = esdZDC->GetImpactParameter();
465 fNpartZDC = esdZDC->GetZDCParticipants();
466 fbZDCA = esdZDC->GetImpactParamSideA();
467 fNpartZDCA = esdZDC->GetZDCPartSideA();
468 fbZDCC = esdZDC->GetImpactParamSideC();
469 fNpartZDCC = esdZDC->GetZDCPartSideC();
470
471 const Double_t * towZNC = esdZDC->GetZN1TowerEnergy();
472 const Double_t * towZPC = esdZDC->GetZP1TowerEnergy();
473 const Double_t * towZNA = esdZDC->GetZN2TowerEnergy();
474 const Double_t * towZPA = esdZDC->GetZP2TowerEnergy();
475 //
476 for(Int_t it=0; it<5; it++){
477 fZNCtower[it] = (Float_t) (towZNC[it]);
478 fZPCtower[it] = (Float_t) (towZPC[it]);
479 fZNAtower[it] = (Float_t) (towZNA[it]);
480 fZPAtower[it] = (Float_t) (towZPA[it]);
481 }
482
483 Double_t xyZNC[2]={-99.,-99.}, xyZNA[2]={-99.,-99.};
484 esdZDC->GetZNCentroidInPbPb(fBeamEnergy, xyZNC, xyZNA);
485 for(Int_t it=0; it<2; it++){
486 fCentrZNC[it] = xyZNC[it];
487 fCentrZNA[it] = xyZNA[it];
488 }
489
490 if(fTreeFilling){
491 fCentralityTree->Fill();
492 PostData(2, fCentralityTree);
493 }
494 }
495 else if(fAnalysisInput.CompareTo("AOD")==0){
496 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
497 // to be implemented
498 printf(" AOD analysis not yet implemented!!!\n\n");
499 return;
500 }
501
502}
503
504
505
506//________________________________________________________________________
507void AliAnalysisTaskCentralityTreeMaker::Terminate(Option_t */*option*/)
508{
509 // Terminate analysis
510 //
511/* if(fDebug > 1) printf(" **** AliAnalysisTaskCentralityTreeMaker::Terminate() \n");
512
513 fCentralityTree = dynamic_cast<TTree*> (GetOutputData(1));
514 if(!fCentralityTree) printf("ERROR: fCentralityTree not available\n");
515
516 //fOutput = dynamic_cast<TList*> (GetOutputData(1));
517 //if(!fOutput) printf("ERROR: fOutput not available\n");
518*/
519}
520
521