1 /**************************************************************************
2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
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 **************************************************************************/
16 /////////////////////////////////////////////////////////////
18 // Class to analyze centrality measurements //
20 /////////////////////////////////////////////////////////////
29 #include "AliAnalysisManager.h"
30 #include "AliVEvent.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"
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"
57 ClassImp(AliAnalysisTaskCentralityTreeMaker)
60 //________________________________________________________________________
61 AliAnalysisTaskCentralityTreeMaker::AliAnalysisTaskCentralityTreeMaker():
64 fAnalysisInput("ESD"),
68 fEZDCvsNtracklets(0x0),
107 // Default constructor
110 //________________________________________________________________________
111 AliAnalysisTaskCentralityTreeMaker::AliAnalysisTaskCentralityTreeMaker(const char *name):
112 AliAnalysisTaskSE(name),
114 fAnalysisInput("ESD"),
118 fEZDCvsNtracklets(0x0),
120 fCentralityTree(0x0),
157 // Default constructor
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());
166 //________________________________________________________________________
167 AliAnalysisTaskCentralityTreeMaker& AliAnalysisTaskCentralityTreeMaker::operator=(const AliAnalysisTaskCentralityTreeMaker& c)
170 // Assignment operator
173 AliAnalysisTaskSE::operator=(c);
178 //________________________________________________________________________
179 AliAnalysisTaskCentralityTreeMaker::AliAnalysisTaskCentralityTreeMaker(const AliAnalysisTaskCentralityTreeMaker& ana):
180 AliAnalysisTaskSE(ana),
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),
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),
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),
206 fNpartZDC(ana.fNpartZDC),
208 fNpartZDCA(ana.fNpartZDCA),
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)
230 //________________________________________________________________________
231 AliAnalysisTaskCentralityTreeMaker::~AliAnalysisTaskCentralityTreeMaker()
235 delete fOutput; fOutput=0;
238 delete fCentralityTree; fCentralityTree=0;
243 //________________________________________________________________________
244 void AliAnalysisTaskCentralityTreeMaker::UserCreateOutputObjects()
247 // Create the output containers
248 if(fDebug>1) printf("AnalysisTaskZDCpp::UserCreateOutputObjects() \n");
250 // Several histograms are more conveniently managed in a TList
251 fOutput = new TList();
253 fOutput->SetName("OutputHistos");
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);
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);
267 fCentralityTree = new TTree("fCentralityTree", "Centrality vs. multiplicity tree");
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");
316 //________________________________________________________________________
317 void AliAnalysisTaskCentralityTreeMaker::UserExec(Option_t */*option*/)
319 // Execute analysis for current event:
320 if(fDebug>1) printf(" **** AliAnalysisTaskCentralityTreeMaker::UserExec() \n");
322 if(fAnalysisInput.CompareTo("ESD")==0){
324 // AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
325 AliVEvent* event = InputEvent();
326 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
330 fNTracks = event->GetNumberOfTracks();
331 fNPmdTracks = esd->GetNumberOfPmdTracks();
333 AliESDVZERO* esdV0 = esd->GetVZEROData();
334 fMultV0A=esdV0->GetMTotV0A();
335 fMultV0C=esdV0->GetMTotV0C();
339 AliMCEvent* mcEvent = MCEvent();
341 printf(" Could not retrieve MC event!!!\n");
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);
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;
363 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
365 printf(" Event generator header not available!!!\n");
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();
385 fBeamEnergy = esd->GetBeamEnergy();
387 // ***** Trigger selection
388 TString triggerClass = esd->GetFiredTriggerClasses();
389 sprintf(fTrigClass,"%s",triggerClass.Data());
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;
398 vertex->GetXYZ(vertex3);
400 const AliMultiplicity *mult = esd->GetMultiplicity();
401 fNTracklets = mult->GetNumberOfTracklets();
403 for(Int_t ilay=0; ilay<6; ilay++){
404 fNClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
406 fNSingleClusters = mult->GetNumberOfSingleClusters();
408 for(Int_t ilay=0; ilay<2; ilay++){
409 fNChips[ilay] = mult->GetNumberOfFiredChips(ilay);
413 AliESDFMD *fmd = esd->GetFMDData();
415 Float_t totalMultA = 0;
416 Float_t totalMultC = 0;
417 const Float_t fFMDLowCut = 0.4;
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++) {
428 Float_t FMDmult = fmd->Multiplicity(det,ring,sec,strip);
429 if(FMDmult == 0 || FMDmult == AliESDFMD::kInvalidMult) continue;
431 Float_t nParticles=0;
433 if(FMDmult > fFMDLowCut) {
437 if (det<3) totalMultA = totalMultA + nParticles;
438 else totalMultC = totalMultC + nParticles;
444 fMultFMDA = totalMultA;
445 fMultFMDC = totalMultC;
448 AliESDZDC *esdZDC = esd->GetESDZDC();
450 fESDFlag = esdZDC->GetESDQuality();
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));
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);
464 fbZDC = esdZDC->GetImpactParameter();
465 fNpartZDC = esdZDC->GetZDCParticipants();
466 fbZDCA = esdZDC->GetImpactParamSideA();
467 fNpartZDCA = esdZDC->GetZDCPartSideA();
468 fbZDCC = esdZDC->GetImpactParamSideC();
469 fNpartZDCC = esdZDC->GetZDCPartSideC();
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();
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]);
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];
491 fCentralityTree->Fill();
492 PostData(2, fCentralityTree);
495 else if(fAnalysisInput.CompareTo("AOD")==0){
496 //AliAODEvent *aod = dynamic_cast<AliAODEvent*> (InputEvent());
498 printf(" AOD analysis not yet implemented!!!\n\n");
506 //________________________________________________________________________
507 void AliAnalysisTaskCentralityTreeMaker::Terminate(Option_t */*option*/)
509 // Terminate analysis
511 /* if(fDebug > 1) printf(" **** AliAnalysisTaskCentralityTreeMaker::Terminate() \n");
513 fCentralityTree = dynamic_cast<TTree*> (GetOutputData(1));
514 if(!fCentralityTree) printf("ERROR: fCentralityTree not available\n");
516 //fOutput = dynamic_cast<TList*> (GetOutputData(1));
517 //if(!fOutput) printf("ERROR: fOutput not available\n");