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 // AliAnalysisTaskSE to make AOD centrality
19 // Author: Alberica Toia, CERN, Alberica.Toia@cern.ch
21 /////////////////////////////////////////////////////////////
26 #include "AliVEvent.h"
27 #include "AliESDEvent.h"
28 #include "AliAODEvent.h"
29 #include "AliAODCentrality.h"
30 #include "AliAnalysisManager.h"
31 #include "AliESDInputHandler.h"
32 #include "AliESDZDC.h"
33 #include "AliESDFMD.h"
34 #include "AliESDVZERO.h"
35 #include "AliMultiplicity.h"
36 #include "AliAODHandler.h"
37 #include "AliAODHeader.h"
38 #include "AliAODEvent.h"
39 #include "AliAODVertex.h"
40 #include "AliMCEvent.h"
41 #include "AliMCEventHandler.h"
42 #include "AliMCParticle.h"
44 #include "AliHeader.h"
45 #include "AliGenEventHeader.h"
46 #include "AliGenHijingEventHeader.h"
47 #include "AliAnalysisTaskAODCentralityMaker.h"
50 ClassImp(AliAnalysisTaskAODCentralityMaker)
53 //________________________________________________________________________
54 AliAnalysisTaskAODCentralityMaker::AliAnalysisTaskAODCentralityMaker():
57 fDeltaAODFileName("AliAOD.Centrality.root"),
61 // Default constructor
64 //________________________________________________________________________
65 AliAnalysisTaskAODCentralityMaker::AliAnalysisTaskAODCentralityMaker(const char *name):
66 AliAnalysisTaskSE(name),
68 fDeltaAODFileName("AliAOD.Centrality.root"),
72 // Standard constructor
76 //________________________________________________________________________
77 AliAnalysisTaskAODCentralityMaker::~AliAnalysisTaskAODCentralityMaker()
82 //________________________________________________________________________
83 void AliAnalysisTaskAODCentralityMaker::Init()
86 if(fDebug > 1) printf("AnalysisTaskAODCentralityMaker::Init() \n");
87 AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFileName.Data());
92 //________________________________________________________________________
93 void AliAnalysisTaskAODCentralityMaker::UserCreateOutputObjects()
96 // Create the output container
98 if(fDebug > 1) printf("AnalysisTaskAODCentralityMaker::UserCreateOutPutData() \n");
99 // Support both the case when the AOD + deltaAOD are produced in an ESD
100 // analysis or if the deltaAOD is produced on an analysis on AOD's. (A.G. 27/04/09)
101 if(!AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()) {
102 Fatal("UserCreateOutputObjects", "This task needs an AOD handler");
105 TString filename = fDeltaAODFileName;
106 // When running on standard AOD to produce deltas, IsStandardAOD is never set,
107 // If AODEvent is NULL, new branches have to be added to the new file(s) (A.G. 15/01/10)
108 if(!IsStandardAOD() && AODEvent()) filename = "";
110 fAODCentrality = new AliAODCentrality();
111 fAODCentrality->SetName("AODCentrality");
112 AddAODBranch("AliAODCentrality", &fAODCentrality, filename);
115 fAODHeader = new AliAODHeader();
116 AddAODBranch("AliAODHeader", &fAODHeader, filename);
121 //________________________________________________________________________
122 void AliAnalysisTaskAODCentralityMaker::UserExec(Option_t */*option*/)
125 AliVEvent* event = InputEvent();
126 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
128 AliError("No ESD Event");
132 Float_t beamEnergy = esd->GetBeamEnergy();
133 Int_t nTracks = event->GetNumberOfTracks();
134 Int_t nPmdTracks = esd->GetNumberOfPmdTracks();
137 AliESDVZERO* esdV0 = esd->GetVZEROData();
138 Double_t multV0A = esdV0->GetMTotV0A();
139 Double_t multV0C = esdV0->GetMTotV0C();
143 const AliESDVertex *vertex = esd->GetPrimaryVertexSPD();
144 Double_t xVertex = vertex->GetX();
145 Double_t yVertex = vertex->GetY();
146 Double_t zVertex = vertex->GetZ();
149 if(vertex->IsFromVertexer3D()) vertexer3d = kTRUE;
150 else vertexer3d = kFALSE;
152 vertex->GetXYZ(vertex3);
154 // ***** CB info (tracklets, clusters, chips)
155 const AliMultiplicity *mult = esd->GetMultiplicity();
156 Int_t nTracklets = mult->GetNumberOfTracklets();
157 Int_t nSingleClusters;
160 for(Int_t ilay = 0; ilay < 6; ilay++){
161 nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
163 nSingleClusters = mult->GetNumberOfSingleClusters();
166 for(Int_t ilay = 0; ilay < 2; ilay++){
167 nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
171 AliESDFMD *fmd = esd->GetFMDData();
172 Float_t totalMultA = 0;
173 Float_t totalMultC = 0;
174 const Float_t fmdLowCut = 0.4;
176 for(UShort_t det = 1;det <= 3; det++) {
177 Int_t nRings = (det==1 ? 1 : 2);
178 for (UShort_t ir = 0; ir < nRings; ir++) {
179 Char_t ring = (ir == 0 ? 'I' : 'O');
180 UShort_t nsec = (ir == 0 ? 20 : 40);
181 UShort_t nstr = (ir == 0 ? 512 : 256);
182 for(UShort_t sec =0; sec < nsec; sec++) {
183 for(UShort_t strip = 0; strip < nstr; strip++) {
184 Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
185 if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
186 Float_t nParticles=0;
187 if(fmdMult > fmdLowCut) {
191 if (det<3) totalMultA = totalMultA + nParticles;
192 else totalMultC = totalMultC + nParticles;
198 Float_t multFMDA = totalMultA;
199 Float_t multFMDC = totalMultC;
202 AliESDZDC *esdZDC = esd->GetESDZDC();
203 UInt_t esdFlag = esdZDC->GetESDQuality();
205 Float_t znCEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
206 Float_t zpCEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
207 Float_t znAEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
208 Float_t zpAEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
209 Float_t zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
210 Float_t zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
212 Double_t bZDC = esdZDC->GetImpactParameter();
213 Int_t nPartZDC = esdZDC->GetZDCParticipants();
214 Double_t bZDCA = esdZDC->GetImpactParamSideA();
215 Int_t nPartZDCA = esdZDC->GetZDCPartSideA();
216 Double_t bZDCC = esdZDC->GetImpactParamSideC();
217 Int_t nPartZDCC = esdZDC->GetZDCPartSideC();
219 const Double_t * towZNC = esdZDC->GetZN1TowerEnergy();
220 const Double_t * towZPC = esdZDC->GetZP1TowerEnergy();
221 const Double_t * towZNA = esdZDC->GetZN2TowerEnergy();
222 const Double_t * towZPA = esdZDC->GetZP2TowerEnergy();
224 Float_t znCtower[5]; // ZNC 5 tower signals
225 Float_t zpCtower[5]; // ZPC 5 tower signals
226 Float_t znAtower[5]; // ZNA 5 tower signals
227 Float_t zpAtower[5]; // ZPA 5 tower signals
228 Float_t centrZNC[2]; // centroid over ZNC
229 Float_t centrZNA[2]; // centroid over ZNA
231 for(Int_t it = 0; it < 5; it++){
232 znCtower[it] = (Float_t) (towZNC[it]);
233 zpCtower[it] = (Float_t) (towZPC[it]);
234 znAtower[it] = (Float_t) (towZNA[it]);
235 zpAtower[it] = (Float_t) (towZPA[it]);
238 Double_t xyZNC[2] = {-99.,-99.};
239 Double_t xyZNA[2] = {-99.,-99.};
241 esdZDC->GetZNCentroidInPbPb(beamEnergy, xyZNC, xyZNA);
242 for(Int_t it = 0; it < 2; it++){
243 centrZNC[it] = xyZNC[it];
244 centrZNA[it] = xyZNA[it];
249 Int_t specNeutronProj = 0;
250 Int_t specProtonProj = 0;
251 Int_t specNeutronTarg = 0;
252 Int_t specProtonTarg = 0;
253 Int_t nPartTargMC = 0;
254 Int_t nPartProjMC = 0;
262 AliMCEvent* mcEvent = MCEvent();
264 printf(" Could not retrieve MC event!!!\n");
268 Int_t nMyTracks_gen = 0;
269 AliStack *stack = 0x0; // needed for MC studies
270 stack = MCEvent()->Stack();
271 for (Int_t iTrack = 0; iTrack < MCEvent()->GetNumberOfTracks(); iTrack++) {
272 //get properties of mc particle
273 AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(iTrack);
275 if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;
276 //charged tracks only
277 if (mcP->Particle()->GetPDG()->Charge() == 0) continue;
278 //same cuts as on ESDtracks
279 // if(TMath::Abs(mcP->Eta())>0.9)continue;
280 // if(mcP->Pt()<0.2)continue;
281 // if(mcP->Pt()>200)continue;
286 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
288 printf(" Event generator header not available!!!\n");
293 if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
294 bMC = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
295 specNeutronProj = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsn();
296 specProtonProj = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsp();
297 specNeutronTarg = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsn();
298 specProtonTarg = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsp();
299 nPartTargMC = Int_t (208.-(specNeutronTarg+specProtonTarg));
300 nPartProjMC = Int_t (208.-(specNeutronProj+specProtonProj));
301 nnColl = ((AliGenHijingEventHeader*) genHeader)->NN();
302 nnwColl = ((AliGenHijingEventHeader*) genHeader)->NNw();
303 nwNColl = ((AliGenHijingEventHeader*) genHeader)->NwN();
304 nwNwColl = ((AliGenHijingEventHeader*) genHeader)->NwNw();
308 fAODCentrality->SetxVertex (xVertex );
309 fAODCentrality->SetyVertex (yVertex );
310 fAODCentrality->SetzVertex (zVertex );
311 fAODCentrality->SetVertexer3d (vertexer3d );
312 fAODCentrality->SetbMC (bMC);
313 fAODCentrality->SetNpartTargMC (nPartTargMC);
314 fAODCentrality->SetNpartProjMC (nPartProjMC);
315 fAODCentrality->SetNNColl (nnColl);
316 fAODCentrality->SetNNwColl (nnwColl);
317 fAODCentrality->SetNwNColl (nwNColl);
318 fAODCentrality->SetNwNwColl (nwNwColl);
319 fAODCentrality->SetNTracklets (nTracklets);
320 fAODCentrality->SetNSingleClusters (nSingleClusters);
321 fAODCentrality->SetNClusters (
328 fAODCentrality->SetNChips (
331 fAODCentrality->SetbZDC (bZDC);
332 fAODCentrality->SetNpartZDC (nPartZDC);
333 fAODCentrality->SetbZDCA (bZDCA);
334 fAODCentrality->SetNpartZDCA (nPartZDCA);
335 fAODCentrality->SetbZDCC (bZDCC);
336 fAODCentrality->SetNpartZDCC (nPartZDCC);
337 fAODCentrality->SetESDFlag (esdFlag);
338 fAODCentrality->SetZNCEnergy (znCEnergy);
339 fAODCentrality->SetZPCEnergy (zpCEnergy);
340 fAODCentrality->SetZNAEnergy (znAEnergy);
341 fAODCentrality->SetZPAEnergy (zpAEnergy);
342 fAODCentrality->SetZEM1Energy (zem1Energy);
343 fAODCentrality->SetZEM2Energy (zem2Energy);
344 fAODCentrality->SetZNCtower (
350 fAODCentrality->SetZPCtower (
357 fAODCentrality-> SetZNAtower (
363 fAODCentrality-> SetZPAtower (
369 fAODCentrality-> SetCentrZNC (
372 fAODCentrality-> SetCentrZNA (
375 fAODCentrality-> SetNTracks (nTracks);
376 fAODCentrality-> SetNPmdTracks (nPmdTracks);
377 fAODCentrality-> SetMultV0A (multV0A);
378 fAODCentrality-> SetMultV0C (multV0C);
379 fAODCentrality-> SetMultFMDA (multFMDA);
380 fAODCentrality-> SetMultFMDC (multFMDC);
383 // Header Replication
385 AliAODHeader* hdr = AODEvent()->GetHeader();
391 //________________________________________________________________________
392 void AliAnalysisTaskAODCentralityMaker::Terminate(Option_t */*option*/)
394 // Terminate analysis
396 if(fDebug > 1) printf("AnalysisTaskAODCentralityMaker: Terminate() \n");