Fix in the last caall to CleanOwnPrimaryVertex
[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 // User Exec
125   AliVEvent*   event = InputEvent();
126   AliESDEvent* esd   = dynamic_cast<AliESDEvent*>(event);
127   if (!esd) {
128       AliError("No ESD Event");
129       return;
130   }
131   
132   Float_t beamEnergy = esd->GetBeamEnergy();
133   Int_t   nTracks    = event->GetNumberOfTracks();     
134   Int_t   nPmdTracks = esd->GetNumberOfPmdTracks();     
135     
136   // ***** V0 info
137   AliESDVZERO* esdV0 = esd->GetVZEROData();
138   Double_t multV0A = esdV0->GetMTotV0A();
139   Double_t multV0C = esdV0->GetMTotV0C();
140     
141   
142   // ***** vertex info
143   const AliESDVertex *vertex = esd->GetPrimaryVertexSPD();
144   Double_t xVertex = vertex->GetX();
145   Double_t yVertex = vertex->GetY();
146   Double_t zVertex = vertex->GetZ();
147   Bool_t vertexer3d;
148   
149   if(vertex->IsFromVertexer3D()) vertexer3d = kTRUE;
150   else vertexer3d = kFALSE;
151   Double_t vertex3[3];
152   vertex->GetXYZ(vertex3);
153   
154   // ***** CB info (tracklets, clusters, chips)
155   const AliMultiplicity *mult = esd->GetMultiplicity();
156   Int_t nTracklets = mult->GetNumberOfTracklets();
157   Int_t nSingleClusters;
158   Int_t nClusters[6];
159   
160   for(Int_t ilay = 0; ilay < 6; ilay++){
161     nClusters[ilay] = mult->GetNumberOfITSClusters(ilay);
162   }
163   nSingleClusters = mult->GetNumberOfSingleClusters();
164
165   Int_t nChips[2];
166   for(Int_t ilay = 0; ilay < 2; ilay++){
167     nChips[ilay] = mult->GetNumberOfFiredChips(ilay);
168   }
169   
170   // ***** FMD info
171   AliESDFMD *fmd = esd->GetFMDData();
172   Float_t totalMultA = 0;
173   Float_t totalMultC = 0;
174   const Float_t fmdLowCut = 0.4;
175   
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) {
188                       nParticles = 1.;
189                   }
190           
191                   if (det<3) totalMultA = totalMultA + nParticles;
192                   else totalMultC = totalMultC + nParticles;
193                   
194               }
195           }
196       }
197   }
198   Float_t multFMDA = totalMultA;
199   Float_t multFMDC = totalMultC;
200   
201   // ***** ZDC info
202   AliESDZDC *esdZDC = esd->GetESDZDC();
203   UInt_t esdFlag =  esdZDC->GetESDQuality();   
204   
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));
211   
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();
218   
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();
223   //
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
230
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]);  
236   }
237   
238   Double_t xyZNC[2] = {-99.,-99.};
239   Double_t xyZNA[2] = {-99.,-99.};
240
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];
245   }
246
247   // ***** MC info
248   Double_t bMC          = 0.;
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;
255   Int_t nnColl          = 0;
256   Int_t nnwColl         = 0;
257   Int_t nwNColl         = 0;
258   Int_t nwNwColl        = 0;
259
260   if(fIsMCInput){
261     
262     AliMCEvent* mcEvent = MCEvent();
263     if (!mcEvent) {
264       printf("   Could not retrieve MC event!!!\n");
265       return;
266     }
267     
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);
274       // Primaries only
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;
282       
283       nMyTracks_gen ++;
284     } 
285     
286     AliGenEventHeader* genHeader = mcEvent->GenEventHeader();
287     if(!genHeader){
288       printf("  Event generator header not available!!!\n");
289       return;
290     }
291         
292
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();
305     }  
306   }
307   
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         (
322                                            nClusters[0],
323                                            nClusters[1],
324                                            nClusters[2],
325                                            nClusters[3],
326                                            nClusters[4],
327                                            nClusters[5]);
328   fAODCentrality->SetNChips            (
329                                            nChips[0],
330                                            nChips[1]);
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          (
345                                            znCtower[0],
346                                            znCtower[1],
347                                            znCtower[2],
348                                            znCtower[3],
349                                            znCtower[4]);
350   fAODCentrality->SetZPCtower          (
351                                          zpCtower[0],
352                                          zpCtower[1],
353                                          zpCtower[2],
354                                          zpCtower[3],
355                                          zpCtower[4]);
356  
357   fAODCentrality-> SetZNAtower          (
358                              znAtower[0],  
359                              znAtower[1], 
360                              znAtower[2],  
361                              znAtower[3],  
362                              znAtower[4]); 
363   fAODCentrality-> SetZPAtower          (
364                                             zpAtower[0], 
365                                             zpAtower[1], 
366                                             zpAtower[2], 
367                                             zpAtower[3],
368                                             zpAtower[4]);
369   fAODCentrality-> SetCentrZNC          (
370                                             centrZNC[0],
371                                             centrZNC[1]);
372   fAODCentrality-> SetCentrZNA          (
373                                             centrZNA[0],
374                                             centrZNA[1]);
375   fAODCentrality-> SetNTracks           (nTracks);
376   fAODCentrality-> SetNPmdTracks        (nPmdTracks);
377   fAODCentrality-> SetMultV0A           (multV0A);
378   fAODCentrality-> SetMultV0C           (multV0C);
379   fAODCentrality-> SetMultFMDA          (multFMDA);
380   fAODCentrality-> SetMultFMDC          (multFMDC);
381
382 //
383 // Header Replication
384 //
385   AliAODHeader* hdr = AODEvent()->GetHeader();
386   *fAODHeader =  *hdr;
387 //
388   return;
389 }
390
391 //________________________________________________________________________
392 void AliAnalysisTaskAODCentralityMaker::Terminate(Option_t */*option*/)
393 {
394   // Terminate analysis
395   //
396   if(fDebug > 1) printf("AnalysisTaskAODCentralityMaker: Terminate() \n");
397 }