Add forward declaration for ROOT
[u/mrichter/AliRoot.git] / ANALYSIS / AliAnalysisTaskMCParticleFilter.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, 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 //  Analysis task for Kinematic filtering
18 //  Fill AODtrackMC tracks from Kinematic stack
19 //
20  
21 #include <TChain.h>
22 #include <TFile.h>
23 #include "TParticle.h"
24 #include "TParameter.h"
25 #include "TString.h"
26 #include "TList.h"
27 #include "TProfile.h"
28 #include "TH1F.h"
29
30
31 #include "AliAnalysisTaskMCParticleFilter.h"
32 #include "AliAnalysisManager.h"
33 #include "AliAnalysisFilter.h"
34 #include "AliHeader.h"
35 #include "AliStack.h"
36 #include "AliMCEvent.h"
37 #include "AliMCEventHandler.h"
38 #include "AliAODEvent.h"
39 #include "AliAODHeader.h"
40 #include "AliAODMCHeader.h"
41 #include "AliAODHandler.h"
42 #include "AliAODVertex.h"
43 #include "AliAODMCParticle.h"
44 #include "AliCollisionGeometry.h"
45 #include "AliGenDPMjetEventHeader.h"
46 #include "AliGenHijingEventHeader.h"
47 #include "AliGenPythiaEventHeader.h"
48 #include "AliGenCocktailEventHeader.h"
49
50 #include "AliLog.h"
51
52 ClassImp(AliAnalysisTaskMCParticleFilter)
53
54 ////////////////////////////////////////////////////////////////////////
55
56 //____________________________________________________________________
57 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter():
58 AliAnalysisTaskSE(),
59   fTrackFilterMother(0x0),
60   fHistList(0x0)
61 {
62   // Default constructor
63 }
64
65 Bool_t AliAnalysisTaskMCParticleFilter::Notify()
66 {
67   //
68   // Implemented Notify() to read the cross sections
69   // from pyxsec.root
70   // 
71   TTree *tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
72   Double_t xsection = 0;
73   UInt_t   ntrials  = 0;
74   if(tree){
75     TFile *curfile = tree->GetCurrentFile();
76     if (!curfile) {
77       Error("Notify","No current file");
78       return kFALSE;
79     }
80
81     TString fileName(curfile->GetName());
82     if(fileName.Contains("AliESDs.root")){
83         fileName.ReplaceAll("AliESDs.root", "pyxsec.root");
84     }
85     else if(fileName.Contains("AliAOD.root")){
86         fileName.ReplaceAll("AliAOD.root", "pyxsec.root");
87     }
88     else if(fileName.Contains("galice.root")){
89         // for running with galice and kinematics alone...                      
90         fileName.ReplaceAll("galice.root", "pyxsec.root");
91     }
92
93
94     TFile *fxsec = TFile::Open(fileName.Data());
95     if(!fxsec){
96       AliInfo(Form("%s:%d %s not found in the Input",(char*)__FILE__,__LINE__,fileName.Data()));
97       // not a severe condition we just do not have the information...
98       return kTRUE;
99     }
100     TTree *xtree = (TTree*)fxsec->Get("Xsection");
101     if(!xtree){
102       AliWarning(Form("%s:%d tree not found in the pyxsec.root",(char*)__FILE__,__LINE__));
103       return kTRUE;
104     }
105     xtree->SetBranchAddress("xsection",&xsection);
106     xtree->SetBranchAddress("ntrials",&ntrials);
107     xtree->GetEntry(0);
108     ((TProfile*)(fHistList->FindObject("h1Xsec")))->Fill("<#sigma>",xsection);
109   }
110   return kTRUE;
111 }
112
113
114
115 //____________________________________________________________________
116 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const char* name):
117     AliAnalysisTaskSE(name),
118     fTrackFilterMother(0x0),
119     fHistList(0x0)
120 {
121   // Default constructor
122   DefineOutput(1, TList::Class());  
123 }
124
125 /*
126 //____________________________________________________________________
127 AliAnalysisTaskMCParticleFilter::AliAnalysisTaskMCParticleFilter(const AliAnalysisTaskMCParticleFilter& obj):
128     AliAnalysisTaskSE(obj),
129     fTrackFilterMother(obj.fTrackFilterMother)
130 {
131   // Copy constructor
132 }
133 */
134 //____________________________________________________________________
135 AliAnalysisTaskMCParticleFilter::~AliAnalysisTaskMCParticleFilter()
136 {
137   //  if( fTrackFilterMother ) delete fTrackFilterMother;
138 }
139
140 /*
141 //____________________________________________________________________
142 AliAnalysisTaskMCParticleFilter& AliAnalysisTaskMCParticleFilter::operator=(const AliAnalysisTaskMCParticleFilter& other)
143 {
144 // Assignment
145   if(this!=&other) {
146     AliAnalysisTaskSE::operator=(other);
147     fTrackFilterMother = other.fTrackFilterMother;
148   }
149   return *this;
150 }
151 */
152 //____________________________________________________________________
153 void AliAnalysisTaskMCParticleFilter::UserCreateOutputObjects()
154 {
155   // Create the output container
156     if (OutputTree()&&fTrackFilterMother) 
157         OutputTree()->GetUserInfo()->Add(fTrackFilterMother);
158
159     // this part is mainly needed to set the MCEventHandler
160     // to the AODHandler, this will not be needed when
161     // AODHandler goes to ANALYSISalice
162     // setting in the steering macro will not work on proof :(
163     // so we have to do it in a task
164
165     // the branch booking can also go into the AODHandler then
166
167
168     // mcparticles
169     TClonesArray *tca = new TClonesArray("AliAODMCParticle", 0);
170     tca->SetName(AliAODMCParticle::StdBranchName());
171     AddAODBranch("TClonesArray",&tca);
172
173     // MC header...
174     AliAODMCHeader *mcHeader = new AliAODMCHeader();
175     mcHeader->SetName(AliAODMCHeader::StdBranchName());
176     AddAODBranch("AliAODMCHeader",&mcHeader);    
177
178     
179
180     AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
181     AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
182     if(!aodH){
183       AliWarning("Could not get AODHandler");
184       return;
185     }
186     aodH->SetMCEventHandler(mcH);
187
188
189     // these are histograms, for averaging and summing
190     // do it via histograms to be PROOF-proof
191     // which merges the results from different workers
192     // histos are not saved reside only in memory
193     
194     
195
196     fHistList = new TList();
197     TProfile *h1Xsec = new TProfile("h1Xsec","xsec from pyxsec.root",1,0,1);
198     h1Xsec->GetXaxis()->SetBinLabel(1,"<#sigma>");
199     fHistList->Add(h1Xsec);
200     TH1F* h1Trials = new TH1F("h1Trials","trials from MC header",1,0,1);
201     h1Trials->GetXaxis()->SetBinLabel(1,"#sum{ntrials}");
202     fHistList->Add(h1Trials);
203
204     
205 }
206
207 //____________________________________________________________________
208 void AliAnalysisTaskMCParticleFilter::UserExec(Option_t */*option*/)
209 {
210 // Execute analysis for current event
211 //
212
213 // Fill AOD tracks from Kinematic stack
214     
215   // get AliAOD Event 
216   AliAODEvent* aod = AODEvent();
217   if (!aod) {
218       AliWarning("No Output Handler connected, doing nothing !") ;
219       return;
220   }
221
222   AliMCEventHandler *mcH = (AliMCEventHandler*) ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler()); 
223   if(!mcH){ 
224     AliWarning("No MC handler Found");
225     return;
226   }
227   
228   AliAODHandler *aodH = dynamic_cast<AliAODHandler*> ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
229   if(!aodH){
230     AliWarning("Could not get AODHandler");
231     return;
232   }
233
234
235
236   // fetch the output 
237   // Fill control histos
238
239   // tmp array for holding the mctracks
240   // need to set mother and duagther __before__
241   // filling in the tree
242
243   AliMCEvent *mcE = MCEvent();
244
245   if(!mcE){
246     AliWarning("No MC event Found");
247     return;
248   }
249
250   Int_t np    = mcE->GetNumberOfTracks();
251   Int_t nprim = mcE->GetNumberOfPrimaries();
252   // TODO ADD MC VERTEX
253
254   AliAODMCHeader *aodMCHo = (AliAODMCHeader *) aod->FindListObject("mcHeader");
255   // Get the proper MC Collision Geometry
256   AliGenEventHeader* mcEH = mcE->GenEventHeader();
257
258   AliGenPythiaEventHeader *pyH  = dynamic_cast<AliGenPythiaEventHeader*>(mcEH);
259   AliGenHijingEventHeader *hiH  = 0;
260   AliCollisionGeometry    *colG = 0;
261   AliGenDPMjetEventHeader *dpmH = 0;
262   // it can be only one save some casts
263   // assuming PYTHIA and HIJING are the most likely ones...
264   if(!pyH){
265       hiH = dynamic_cast<AliGenHijingEventHeader*>(mcEH);
266       if(!hiH){
267           dpmH = dynamic_cast<AliGenDPMjetEventHeader*>(mcEH);
268       }
269   }
270   
271   if (hiH || dpmH) colG = dynamic_cast<AliCollisionGeometry*>(mcEH);
272
273   // fetch the trials on a event by event basis, not from pyxsec.root otherwise 
274   // we will get a problem when running on proof since Notify may be called 
275   // more than once per file
276   // consider storing this information in the AOD output via AliAODHandler
277   Float_t ntrials = 0;
278   if (!colG) {
279     AliGenCocktailEventHeader *ccEH = dynamic_cast<AliGenCocktailEventHeader *>(mcEH);
280     if (ccEH) {
281       TList *genHeaders = ccEH->GetHeaders();
282       for (int imch=0; imch<genHeaders->GetEntries(); imch++) {
283         if(!pyH)dynamic_cast<AliGenPythiaEventHeader*>(genHeaders->At(imch));
284         if(!hiH)dynamic_cast<AliGenHijingEventHeader*>(genHeaders->At(imch));
285         if(!colG)colG = dynamic_cast<AliCollisionGeometry *>(genHeaders->At(imch));
286         if(!dpmH)dynamic_cast<AliGenDPMjetEventHeader*>(genHeaders->At(imch));
287       }
288     }
289   }
290
291   // take the trials from the p+p event
292   if(hiH)ntrials = hiH->Trials();
293   if(dpmH)ntrials = dpmH->Trials();
294   if(pyH)ntrials = pyH->Trials();
295   if(ntrials)((TH1F*)(fHistList->FindObject("h1Trials")))->Fill("#sum{ntrials}",ntrials); 
296   
297
298
299
300   if (colG) {
301     aodMCHo->SetReactionPlaneAngle(colG->ReactionPlaneAngle());
302     AliInfo(Form("Found Collision Geometry. Got Reaction Plane %lf\n", colG->ReactionPlaneAngle()));
303   }
304
305
306
307   // check varibales for charm need all daughters
308   static int  iTaken = 0;
309   static int  iAll = 0;
310   static int  iCharm = 0;
311
312
313   Int_t j=0;
314   for (Int_t ip = 0; ip < np; ip++){
315     AliMCParticle* mcpart = (AliMCParticle*) mcE->GetTrack(ip);
316     TParticle* part = mcpart->Particle();
317     Float_t xv = part->Vx();
318     Float_t yv = part->Vy();
319     Float_t zv = part->Vz();
320     Float_t rv = TMath::Sqrt(xv * xv + yv * yv);
321       
322     Bool_t write = kFALSE;
323     
324     if (ip < nprim) {
325       // Select the primary event
326       write = kTRUE;
327     } else if (part->GetUniqueID() == 4) {
328       // Particles from decay
329       // Check that the decay chain ends at a primary particle
330       AliMCParticle* mother = mcpart;
331       Int_t imo = mcpart->GetMother();
332       while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
333         mother =  (AliMCParticle*) mcE->GetTrack(imo);
334         imo =  mother->GetMother();
335       }
336       // Select according to pseudorapidity and production point of primary ancestor
337       if (imo < nprim && Select(((AliMCParticle*) mcE->GetTrack(imo))->Particle(), rv, zv))write = kTRUE;         
338     } else if (part->GetUniqueID() == 5) {
339       // Now look for pair production
340       Int_t imo = mcpart->GetMother();
341       if (imo < nprim) {
342         // Select, if the gamma is a primary
343         write = kTRUE;
344       } else {
345         // Check if the gamma comes from the decay chain of a primary particle
346         AliMCParticle* mother =  (AliMCParticle*) mcE->GetTrack(imo);
347         imo = mother->GetMother();
348         while((imo >= nprim) && (mother->GetUniqueID() == 4)) {
349           mother =   (AliMCParticle*) mcE->GetTrack(imo);
350           imo =  mother->GetMother();
351         }
352         // Select according to pseudorapidity and production point 
353         if (imo < nprim && Select(mother->Particle(), rv, zv)) 
354           write = kTRUE;
355       }
356     }
357     /*
358     else if (part->GetUniqueID() == 13){
359       // Evaporation
360       // Check that we end at a primary particle
361       TParticle* mother = part;
362       Int_t imo = part->GetFirstMother();
363       while((imo >= nprim) && ((mother->GetUniqueID() == 4 ) || ( mother->GetUniqueID() == 13))) {
364         mother =  mcE->Stack()->Particle(imo);
365         imo =  mother->GetFirstMother();
366       }
367       // Select according to pseudorapidity and production point 
368         if (imo < nprim && Select(mother, rv, zv)) 
369           write = kTRUE;
370     }
371     */    
372     if (write) {
373       if(mcH)mcH->SelectParticle(ip);
374       j++;
375       
376       // debug info to check fro charm daugthers
377       if((TMath::Abs(part->GetPdgCode()))==411){
378         iCharm++;
379         Int_t d0 =  mcpart->GetFirstDaughter();
380         Int_t d1 =  mcpart->GetLastDaughter();
381         if(d0>0&&d1>0){
382           for(int id = d0;id <= d1;id++){
383             iAll++;
384             if(mcH->IsParticleSelected(id))iTaken++;
385           }
386         }
387       }// if charm
388     }
389   }
390
391   AliInfo(Form("Taken daughters %d/%d of %d charm",iTaken,iAll,iCharm));
392
393   aodH->StoreMCParticles();
394   PostData(1,fHistList);
395
396   return;
397 }
398
399 void AliAnalysisTaskMCParticleFilter::Terminate(Option_t */*option*/){
400   //
401   // Terminate the execution save the PYTHIA x_section to the UserInfo()
402   //
403
404
405   fHistList = (TList*)GetOutputData(1);
406   if(!fHistList){
407     Printf("%s:%d Output list not found",(char*)__FILE__,__LINE__);
408     return;
409   }
410
411   TProfile *hXsec = ((TProfile*)(fHistList->FindObject("h1Xsec")));
412   TH1F* hTrials  = ((TH1F*)(fHistList->FindObject("h1Trials")));
413   if(!hXsec)return;
414   if(!hTrials)return;
415
416   Float_t xsec = hXsec->GetBinContent(1);
417   Float_t trials = hTrials->Integral();
418   AliInfo(Form("Average -section %.4E and total number of trials %E",xsec,trials));
419
420 }
421
422 Bool_t AliAnalysisTaskMCParticleFilter::Select(TParticle* part, Float_t rv, Float_t zv)
423 {
424   // Selection accoring to eta of the mother and production point
425   // This has to be refined !!!!!!
426
427   // Esp if we don't have collisison in the central barrel but e.g. beam gas
428   //  return kTRUE;
429
430   Float_t eta = part->Eta();
431
432   // central barrel consider everything in the ITS...
433   // large eta window for smeared vertex and SPD acceptance (2 at z = 0)
434   // larger for V0s in the TPC
435   //  if (TMath::Abs(eta) < 2.5 && rv < 250. && TMath::Abs(zv)<255)return kTRUE;
436
437   if (TMath::Abs(eta) < 2.5 && rv < 170)return kTRUE;   
438
439   // Andreas' Cuts
440   //  if (TMath::Abs(eta) < 1. && rv < 170)return kTRUE;   
441
442
443
444   // Muon arm
445   if(eta > -4.2 && eta < -2.3 && zv > -500.)return kTRUE; // Muon arms
446
447   // PMD acceptance 2.3 <= eta < = 3.5
448   //  if(eta>2.0&&eta<3.7)return kTRUE; 
449
450   return kFALSE;
451  
452 }
453