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