added error messages when task does not post data in UserCreateOutputObjects for...
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskAODCentralityMaker.cxx
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"
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"
43 #include "AliStack.h"
44 #include "AliHeader.h"
45 #include "AliGenEventHeader.h"
46 #include "AliGenHijingEventHeader.h"
47 #include "AliAnalysisTaskAODCentralityMaker.h"
48 #include "AliLog.h"
49
50 ClassImp(AliAnalysisTaskAODCentralityMaker)
51
52
53 //________________________________________________________________________
54 AliAnalysisTaskAODCentralityMaker::AliAnalysisTaskAODCentralityMaker():
55 AliAnalysisTaskSE(),
56 fAODCentrality(0),
57 fDeltaAODFileName("AliAOD.Centrality.root"),
58 fAODHeader        (0),
59 fIsMCInput        (0)
60 {
61   // Default constructor
62 }
63
64 //________________________________________________________________________
65 AliAnalysisTaskAODCentralityMaker::AliAnalysisTaskAODCentralityMaker(const char *name):
66 AliAnalysisTaskSE(name),
67 fAODCentrality(0),
68 fDeltaAODFileName("AliAOD.Centrality.root"),
69 fAODHeader        (0),
70 fIsMCInput        (0)
71 {
72   // Standard constructor
73 }
74
75
76 //________________________________________________________________________
77 AliAnalysisTaskAODCentralityMaker::~AliAnalysisTaskAODCentralityMaker()
78 {
79   // Destructor
80 }  
81
82 //________________________________________________________________________
83 void 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 //________________________________________________________________________
93 void 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   
114
115   fAODHeader = new AliAODHeader();
116   AddAODBranch("AliAODHeader", &fAODHeader, filename);
117   return;
118 }
119
120
121 //________________________________________________________________________
122 void AliAnalysisTaskAODCentralityMaker::UserExec(Option_t */*option*/)
123 {
124   AliVEvent*   event = InputEvent();
125   AliESDEvent* esd   = dynamic_cast<AliESDEvent*>(event);
126   if (!esd) {
127       AliError("No ESD Event");
128       return;
129   }
130   
131   Float_t beamEnergy = esd->GetBeamEnergy();
132   Int_t   nTracks    = event->GetNumberOfTracks();     
133   Int_t   nPmdTracks = esd->GetNumberOfPmdTracks();     
134     
135   // ***** V0 info
136   AliESDVZERO* esdV0 = esd->GetVZEROData();
137   Double_t multV0A = esdV0->GetMTotV0A();
138   Double_t multV0C = esdV0->GetMTotV0C();
139     
140   
141   // ***** vertex info
142   const AliESDVertex *vertex = esd->GetPrimaryVertexSPD();
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;
150   Double_t vertex3[3];
151   vertex->GetXYZ(vertex3);
152   
153   // ***** CB info (tracklets, clusters, chips)
154   const AliMultiplicity *mult = esd->GetMultiplicity();
155   Int_t nTracklets = mult->GetNumberOfTracklets();
156   Int_t nSingleClusters;
157   Int_t nClusters[6];
158   
159   for(Int_t ilay = 0; ilay < 6; ilay++){
160     nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
161   }
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);
167   }
168   
169   // ***** FMD info
170   AliESDFMD *fmd = esd->GetFMDData();
171   Float_t totalMultA = 0;
172   Float_t totalMultC = 0;
173   const Float_t fmdLowCut = 0.4;
174   
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                   }
189           
190                   if (det<3) totalMultA = totalMultA + nParticles;
191                   else totalMultC = totalMultC + nParticles;
192                   
193               }
194           }
195       }
196   }
197   Float_t multFMDA = totalMultA;
198   Float_t multFMDC = totalMultC;
199   
200   // ***** ZDC info
201   AliESDZDC *esdZDC = esd->GetESDZDC();
202   UInt_t esdFlag =  esdZDC->GetESDQuality();   
203   
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));
210   
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();
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   //
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]);  
235   }
236   
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];
244   }
245
246   // ***** MC info
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
259   if(fIsMCInput){
260     
261     AliMCEvent* mcEvent = MCEvent();
262     if (!mcEvent) {
263       printf("   Could not retrieve MC event!!!\n");
264       return;
265     }
266     
267     Int_t nMyTracks_gen = 0;
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       
282       nMyTracks_gen ++;
283     } 
284     
285     AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
286     if(!genHeader){
287       printf("  Event generator header not available!!!\n");
288       return;
289     }
290         
291
292     if(genHeader->InheritsFrom(AliGenHijingEventHeader::Class())){
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();
304     }  
305   }
306   
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);
320   fAODCentrality->SetNClusters         (
321                                            nClusters[0],
322                                            nClusters[1],
323                                            nClusters[2],
324                                            nClusters[3],
325                                            nClusters[4],
326                                            nClusters[5]);
327   fAODCentrality->SetNChips            (
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);
343   fAODCentrality->SetZNCtower          (
344                                            znCtower[0],
345                                            znCtower[1],
346                                            znCtower[2],
347                                            znCtower[3],
348                                            znCtower[4]);
349   fAODCentrality->SetZPCtower          (
350                                          zpCtower[0],
351                                          zpCtower[1],
352                                          zpCtower[2],
353                                          zpCtower[3],
354                                          zpCtower[4]);
355  
356   fAODCentrality-> SetZNAtower          (
357                              znAtower[0],  
358                              znAtower[1], 
359                              znAtower[2],  
360                              znAtower[3],  
361                              znAtower[4]); 
362   fAODCentrality-> SetZPAtower          (
363                                             zpAtower[0], 
364                                             zpAtower[1], 
365                                             zpAtower[2], 
366                                             zpAtower[3],
367                                             zpAtower[4]);
368   fAODCentrality-> SetCentrZNC          (
369                                             centrZNC[0],
370                                             centrZNC[1]);
371   fAODCentrality-> SetCentrZNA          (
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);
380
381 //
382 // Header Replication
383 //
384   AliAODHeader* hdr = AODEvent()->GetHeader();
385   *fAODHeader =  *hdr;
386 //
387   return;
388 }
389
390 //________________________________________________________________________
391 void AliAnalysisTaskAODCentralityMaker::Terminate(Option_t */*option*/)
392 {
393   // Terminate analysis
394   //
395   if(fDebug > 1) printf("AnalysisTaskAODCentralityMaker: Terminate() \n");
396 }