added error messages when task does not post data in UserCreateOutputObjects for...
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskAODCentralityMaker.cxx
CommitLineData
c94a2509 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// AliAnalysisTaskSE to make AOD centrality
19// Author: Alberica Toia, CERN, Alberica.Toia@cern.ch
20//
21/////////////////////////////////////////////////////////////
22
23#include <TROOT.h>
24#include <TSystem.h>
25
26#include "AliVEvent.h"
27#include "AliESDEvent.h"
28#include "AliAODEvent.h"
29#include "AliAODCentrality.h"
c94a2509 30#include "AliAnalysisManager.h"
c94a2509 31#include "AliESDInputHandler.h"
32#include "AliESDZDC.h"
33#include "AliESDFMD.h"
34#include "AliESDVZERO.h"
35#include "AliMultiplicity.h"
36#include "AliAODHandler.h"
42f0d071 37#include "AliAODHeader.h"
c94a2509 38#include "AliAODEvent.h"
39#include "AliAODVertex.h"
c94a2509 40#include "AliMCEvent.h"
41#include "AliMCEventHandler.h"
42#include "AliMCParticle.h"
43#include "AliStack.h"
44#include "AliHeader.h"
c94a2509 45#include "AliGenEventHeader.h"
46#include "AliGenHijingEventHeader.h"
c94a2509 47#include "AliAnalysisTaskAODCentralityMaker.h"
9daedfd1 48#include "AliLog.h"
c94a2509 49
50ClassImp(AliAnalysisTaskAODCentralityMaker)
51
52
53//________________________________________________________________________
54AliAnalysisTaskAODCentralityMaker::AliAnalysisTaskAODCentralityMaker():
55AliAnalysisTaskSE(),
56fAODCentrality(0),
57fDeltaAODFileName("AliAOD.Centrality.root"),
42f0d071 58fAODHeader (0),
59fIsMCInput (0)
c94a2509 60{
61 // Default constructor
c94a2509 62}
63
64//________________________________________________________________________
65AliAnalysisTaskAODCentralityMaker::AliAnalysisTaskAODCentralityMaker(const char *name):
66AliAnalysisTaskSE(name),
67fAODCentrality(0),
68fDeltaAODFileName("AliAOD.Centrality.root"),
42f0d071 69fAODHeader (0),
70fIsMCInput (0)
c94a2509 71{
72 // Standard constructor
c94a2509 73}
74
75
76//________________________________________________________________________
77AliAnalysisTaskAODCentralityMaker::~AliAnalysisTaskAODCentralityMaker()
78{
79 // Destructor
80}
81
82//________________________________________________________________________
83void AliAnalysisTaskAODCentralityMaker::Init()
84{
85 // Initialization
86 if(fDebug > 1) printf("AnalysisTaskAODCentralityMaker::Init() \n");
87 AliAnalysisManager::GetAnalysisManager()->RegisterExtraFile(fDeltaAODFileName.Data());
88
89 return;
90}
91
92//________________________________________________________________________
93void AliAnalysisTaskAODCentralityMaker::UserCreateOutputObjects()
94{
95
96 // Create the output container
97 //
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");
103 return;
104 }
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 = "";
109
110 fAODCentrality = new AliAODCentrality();
111 fAODCentrality->SetName("AODCentrality");
112 AddAODBranch("AliAODCentrality", &fAODCentrality, filename);
113
42f0d071 114
115 fAODHeader = new AliAODHeader();
116 AddAODBranch("AliAODHeader", &fAODHeader, filename);
c94a2509 117 return;
118}
119
42f0d071 120
c94a2509 121//________________________________________________________________________
122void AliAnalysisTaskAODCentralityMaker::UserExec(Option_t */*option*/)
123{
42f0d071 124 AliVEvent* event = InputEvent();
125 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
9daedfd1 126 if (!esd) {
127 AliError("No ESD Event");
128 return;
129 }
130
42f0d071 131 Float_t beamEnergy = esd->GetBeamEnergy();
132 Int_t nTracks = event->GetNumberOfTracks();
133 Int_t nPmdTracks = esd->GetNumberOfPmdTracks();
c94a2509 134
135 // ***** V0 info
136 AliESDVZERO* esdV0 = esd->GetVZEROData();
42f0d071 137 Double_t multV0A = esdV0->GetMTotV0A();
138 Double_t multV0C = esdV0->GetMTotV0C();
c94a2509 139
c94a2509 140
141 // ***** vertex info
142 const AliESDVertex *vertex = esd->GetPrimaryVertexSPD();
42f0d071 143 Double_t xVertex = vertex->GetX();
144 Double_t yVertex = vertex->GetY();
145 Double_t zVertex = vertex->GetZ();
146 Bool_t vertexer3d;
147
148 if(vertex->IsFromVertexer3D()) vertexer3d = kTRUE;
149 else vertexer3d = kFALSE;
c94a2509 150 Double_t vertex3[3];
151 vertex->GetXYZ(vertex3);
152
153 // ***** CB info (tracklets, clusters, chips)
154 const AliMultiplicity *mult = esd->GetMultiplicity();
42f0d071 155 Int_t nTracklets = mult->GetNumberOfTracklets();
156 Int_t nSingleClusters;
157 Int_t nClusters[6];
c94a2509 158
42f0d071 159 for(Int_t ilay = 0; ilay < 6; ilay++){
160 nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
c94a2509 161 }
42f0d071 162 nSingleClusters = mult->GetNumberOfSingleClusters();
163
164 Int_t nChips[2];
165 for(Int_t ilay = 0; ilay < 2; ilay++){
166 nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
c94a2509 167 }
168
169 // ***** FMD info
170 AliESDFMD *fmd = esd->GetFMDData();
171 Float_t totalMultA = 0;
172 Float_t totalMultC = 0;
42f0d071 173 const Float_t fmdLowCut = 0.4;
c94a2509 174
42f0d071 175 for(UShort_t det = 1;det <= 3; det++) {
176 Int_t nRings = (det==1 ? 1 : 2);
177 for (UShort_t ir = 0; ir < nRings; ir++) {
178 Char_t ring = (ir == 0 ? 'I' : 'O');
179 UShort_t nsec = (ir == 0 ? 20 : 40);
180 UShort_t nstr = (ir == 0 ? 512 : 256);
181 for(UShort_t sec =0; sec < nsec; sec++) {
182 for(UShort_t strip = 0; strip < nstr; strip++) {
183 Float_t fmdMult = fmd->Multiplicity(det,ring,sec,strip);
184 if(fmdMult == 0 || fmdMult == AliESDFMD::kInvalidMult) continue;
185 Float_t nParticles=0;
186 if(fmdMult > fmdLowCut) {
187 nParticles = 1.;
188 }
c94a2509 189
42f0d071 190 if (det<3) totalMultA = totalMultA + nParticles;
191 else totalMultC = totalMultC + nParticles;
192
193 }
c94a2509 194 }
c94a2509 195 }
c94a2509 196 }
42f0d071 197 Float_t multFMDA = totalMultA;
198 Float_t multFMDC = totalMultC;
c94a2509 199
200 // ***** ZDC info
201 AliESDZDC *esdZDC = esd->GetESDZDC();
42f0d071 202 UInt_t esdFlag = esdZDC->GetESDQuality();
c94a2509 203
42f0d071 204 Float_t znCEnergy = (Float_t) (esdZDC->GetZDCN1Energy());
205 Float_t zpCEnergy = (Float_t) (esdZDC->GetZDCP1Energy());
206 Float_t znAEnergy = (Float_t) (esdZDC->GetZDCN2Energy());
207 Float_t zpAEnergy = (Float_t) (esdZDC->GetZDCP2Energy());
208 Float_t zem1Energy = (Float_t) (esdZDC->GetZDCEMEnergy(0));
209 Float_t zem2Energy = (Float_t) (esdZDC->GetZDCEMEnergy(1));
c94a2509 210
42f0d071 211 Double_t bZDC = esdZDC->GetImpactParameter();
212 Int_t nPartZDC = esdZDC->GetZDCParticipants();
213 Double_t bZDCA = esdZDC->GetImpactParamSideA();
214 Int_t nPartZDCA = esdZDC->GetZDCPartSideA();
215 Double_t bZDCC = esdZDC->GetImpactParamSideC();
216 Int_t nPartZDCC = esdZDC->GetZDCPartSideC();
c94a2509 217
218 const Double_t * towZNC = esdZDC->GetZN1TowerEnergy();
219 const Double_t * towZPC = esdZDC->GetZP1TowerEnergy();
220 const Double_t * towZNA = esdZDC->GetZN2TowerEnergy();
221 const Double_t * towZPA = esdZDC->GetZP2TowerEnergy();
222 //
42f0d071 223 Float_t znCtower[5]; // ZNC 5 tower signals
224 Float_t zpCtower[5]; // ZPC 5 tower signals
225 Float_t znAtower[5]; // ZNA 5 tower signals
226 Float_t zpAtower[5]; // ZPA 5 tower signals
227 Float_t centrZNC[2]; // centroid over ZNC
228 Float_t centrZNA[2]; // centroid over ZNA
229
230 for(Int_t it = 0; it < 5; it++){
231 znCtower[it] = (Float_t) (towZNC[it]);
232 zpCtower[it] = (Float_t) (towZPC[it]);
233 znAtower[it] = (Float_t) (towZNA[it]);
234 zpAtower[it] = (Float_t) (towZPA[it]);
c94a2509 235 }
236
42f0d071 237 Double_t xyZNC[2] = {-99.,-99.};
238 Double_t xyZNA[2] = {-99.,-99.};
239
240 esdZDC->GetZNCentroidInPbPb(beamEnergy, xyZNC, xyZNA);
241 for(Int_t it = 0; it < 2; it++){
242 centrZNC[it] = xyZNC[it];
243 centrZNA[it] = xyZNA[it];
c94a2509 244 }
245
246 // ***** MC info
42f0d071 247 Double_t bMC = 0.;
248 Int_t specNeutronProj = 0;
249 Int_t specProtonProj = 0;
250 Int_t specNeutronTarg = 0;
251 Int_t specProtonTarg = 0;
252 Int_t nPartTargMC = 0;
253 Int_t nPartProjMC = 0;
254 Int_t nnColl = 0;
255 Int_t nnwColl = 0;
256 Int_t nwNColl = 0;
257 Int_t nwNwColl = 0;
258
c94a2509 259 if(fIsMCInput){
260
261 AliMCEvent* mcEvent = MCEvent();
262 if (!mcEvent) {
263 printf(" Could not retrieve MC event!!!\n");
264 return;
265 }
266
42f0d071 267 Int_t nMyTracks_gen = 0;
c94a2509 268 AliStack *stack = 0x0; // needed for MC studies
269 stack = MCEvent()->Stack();
270 for (Int_t iTrack = 0; iTrack < MCEvent()->GetNumberOfTracks(); iTrack++) {
271 //get properties of mc particle
272 AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(iTrack);
273 // Primaries only
274 if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;
275 //charged tracks only
276 if (mcP->Particle()->GetPDG()->Charge() == 0) continue;
277 //same cuts as on ESDtracks
278 // if(TMath::Abs(mcP->Eta())>0.9)continue;
279 // if(mcP->Pt()<0.2)continue;
280 // if(mcP->Pt()>200)continue;
281
42f0d071 282 nMyTracks_gen ++;
c94a2509 283 }
284
285 AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
286 if(!genHeader){
287 printf(" Event generator header not available!!!\n");
288 return;
289 }
290
42f0d071 291
c94a2509 292 if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
42f0d071 293 bMC = ((AliGenHijingEventHeader*) genHeader)->ImpactParameter();
294 specNeutronProj = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsn();
295 specProtonProj = ((AliGenHijingEventHeader*) genHeader)->ProjSpectatorsp();
296 specNeutronTarg = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsn();
297 specProtonTarg = ((AliGenHijingEventHeader*) genHeader)->TargSpectatorsp();
298 nPartTargMC = Int_t (208.-(specNeutronTarg+specProtonTarg));
299 nPartProjMC = Int_t (208.-(specNeutronProj+specProtonProj));
300 nnColl = ((AliGenHijingEventHeader*) genHeader)->NN();
301 nnwColl = ((AliGenHijingEventHeader*) genHeader)->NNw();
302 nwNColl = ((AliGenHijingEventHeader*) genHeader)->NwN();
303 nwNwColl = ((AliGenHijingEventHeader*) genHeader)->NwNw();
c94a2509 304 }
c94a2509 305 }
306
42f0d071 307 fAODCentrality->SetxVertex (xVertex );
308 fAODCentrality->SetyVertex (yVertex );
309 fAODCentrality->SetzVertex (zVertex );
310 fAODCentrality->SetVertexer3d (vertexer3d );
311 fAODCentrality->SetbMC (bMC);
312 fAODCentrality->SetNpartTargMC (nPartTargMC);
313 fAODCentrality->SetNpartProjMC (nPartProjMC);
314 fAODCentrality->SetNNColl (nnColl);
315 fAODCentrality->SetNNwColl (nnwColl);
316 fAODCentrality->SetNwNColl (nwNColl);
317 fAODCentrality->SetNwNwColl (nwNwColl);
318 fAODCentrality->SetNTracklets (nTracklets);
319 fAODCentrality->SetNSingleClusters (nSingleClusters);
c94a2509 320 fAODCentrality->SetNClusters (
42f0d071 321 nClusters[0],
322 nClusters[1],
323 nClusters[2],
324 nClusters[3],
325 nClusters[4],
326 nClusters[5]);
c94a2509 327 fAODCentrality->SetNChips (
42f0d071 328 nChips[0],
329 nChips[1]);
330 fAODCentrality->SetbZDC (bZDC);
331 fAODCentrality->SetNpartZDC (nPartZDC);
332 fAODCentrality->SetbZDCA (bZDCA);
333 fAODCentrality->SetNpartZDCA (nPartZDCA);
334 fAODCentrality->SetbZDCC (bZDCC);
335 fAODCentrality->SetNpartZDCC (nPartZDCC);
336 fAODCentrality->SetESDFlag (esdFlag);
337 fAODCentrality->SetZNCEnergy (znCEnergy);
338 fAODCentrality->SetZPCEnergy (zpCEnergy);
339 fAODCentrality->SetZNAEnergy (znAEnergy);
340 fAODCentrality->SetZPAEnergy (zpAEnergy);
341 fAODCentrality->SetZEM1Energy (zem1Energy);
342 fAODCentrality->SetZEM2Energy (zem2Energy);
c94a2509 343 fAODCentrality->SetZNCtower (
42f0d071 344 znCtower[0],
345 znCtower[1],
346 znCtower[2],
347 znCtower[3],
348 znCtower[4]);
c94a2509 349 fAODCentrality->SetZPCtower (
42f0d071 350 zpCtower[0],
351 zpCtower[1],
352 zpCtower[2],
353 zpCtower[3],
354 zpCtower[4]);
c94a2509 355
356 fAODCentrality-> SetZNAtower (
42f0d071 357 znAtower[0],
358 znAtower[1],
359 znAtower[2],
360 znAtower[3],
361 znAtower[4]);
c94a2509 362 fAODCentrality-> SetZPAtower (
42f0d071 363 zpAtower[0],
364 zpAtower[1],
365 zpAtower[2],
366 zpAtower[3],
367 zpAtower[4]);
c94a2509 368 fAODCentrality-> SetCentrZNC (
42f0d071 369 centrZNC[0],
370 centrZNC[1]);
c94a2509 371 fAODCentrality-> SetCentrZNA (
42f0d071 372 centrZNA[0],
373 centrZNA[1]);
374 fAODCentrality-> SetNTracks (nTracks);
375 fAODCentrality-> SetNPmdTracks (nPmdTracks);
376 fAODCentrality-> SetMultV0A (multV0A);
377 fAODCentrality-> SetMultV0C (multV0C);
378 fAODCentrality-> SetMultFMDA (multFMDA);
379 fAODCentrality-> SetMultFMDC (multFMDC);
c94a2509 380
42f0d071 381//
382// Header Replication
383//
384 AliAODHeader* hdr = AODEvent()->GetHeader();
385 *fAODHeader = *hdr;
386//
c94a2509 387 return;
388}
389
390//________________________________________________________________________
391void AliAnalysisTaskAODCentralityMaker::Terminate(Option_t */*option*/)
392{
393 // Terminate analysis
394 //
395 if(fDebug > 1) printf("AnalysisTaskAODCentralityMaker: Terminate() \n");
396}