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