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