Corrected protection.
[u/mrichter/AliRoot.git] / STEER / AliAODHandler.cxx
1
2 /**************************************************************************
3  * Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
4  *                                                                        *
5  * Author: The ALICE Off-line Project.                                    *
6  * Contributors are mentioned in the code where appropriate.              *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is          *
14  * provided "as is" without express or implied warranty.                  *
15  **************************************************************************/
16
17 /* $Id$ */
18
19 //-------------------------------------------------------------------------
20 //     Implementation of the Virtual Event Handler Interface for AOD
21 //     Author: Andreas Morsch, CERN
22 //-------------------------------------------------------------------------
23
24
25 #include <TTree.h>
26 #include <TFile.h>
27 #include <TString.h>
28 #include <TList.h>
29 #include <TROOT.h>
30
31 #include "AliLog.h"
32 #include "AliAODHandler.h"
33 #include "AliAODEvent.h"
34 #include "AliAODTracklets.h"
35 #include "AliStack.h"
36 #include "AliAODMCParticle.h"
37 #include "AliAODMCHeader.h"
38 #include "AliMCEventHandler.h"
39 #include "AliMCEvent.h"
40 #include "AliGenEventHeader.h"
41 #include "AliGenHijingEventHeader.h"
42 #include "AliGenDPMjetEventHeader.h"
43 #include "AliGenPythiaEventHeader.h"
44 #include "AliGenCocktailEventHeader.h"
45
46
47
48 ClassImp(AliAODHandler)
49
50 //______________________________________________________________________________
51 AliAODHandler::AliAODHandler() :
52     AliVEventHandler(),
53     fIsStandard(kTRUE),
54     fFillAOD(kTRUE),
55     fFillAODRun(kTRUE),
56     fNeedsHeaderReplication(kFALSE),
57     fNeedsTracksBranchReplication(kFALSE),
58     fNeedsVerticesBranchReplication(kFALSE),
59     fNeedsV0sBranchReplication(kFALSE),
60     fNeedsCascadesBranchReplication(kFALSE),
61     fNeedsTrackletsBranchReplication(kFALSE),
62     fNeedsPMDClustersBranchReplication(kFALSE),
63     fNeedsJetsBranchReplication(kFALSE),
64     fNeedsFMDClustersBranchReplication(kFALSE),
65     fNeedsCaloClustersBranchReplication(kFALSE),
66     fNeedsMCParticlesBranchReplication(kFALSE),
67     fNeedsDimuonsBranchReplication(kFALSE),
68     fAODIsReplicated(kFALSE),
69     fAODEvent(NULL),
70     fMCEventH(NULL),
71     fTreeA(NULL),
72     fFileA(NULL),
73     fFileName(""),
74     fExtensions(NULL),
75     fFilters(NULL)
76 {
77   // default constructor
78 }
79
80 //______________________________________________________________________________
81 AliAODHandler::AliAODHandler(const char* name, const char* title):
82     AliVEventHandler(name, title),
83     fIsStandard(kTRUE),
84     fFillAOD(kTRUE),
85     fFillAODRun(kTRUE),
86     fNeedsHeaderReplication(kFALSE),
87     fNeedsTracksBranchReplication(kFALSE),
88     fNeedsVerticesBranchReplication(kFALSE),
89     fNeedsV0sBranchReplication(kFALSE),
90     fNeedsCascadesBranchReplication(kFALSE),
91     fNeedsTrackletsBranchReplication(kFALSE),
92     fNeedsPMDClustersBranchReplication(kFALSE),
93     fNeedsJetsBranchReplication(kFALSE),
94     fNeedsFMDClustersBranchReplication(kFALSE),
95     fNeedsCaloClustersBranchReplication(kFALSE),
96     fNeedsMCParticlesBranchReplication(kFALSE),
97     fNeedsDimuonsBranchReplication(kFALSE),
98     fAODIsReplicated(kFALSE),
99     fAODEvent(NULL),
100     fMCEventH(NULL),
101     fTreeA(NULL),
102     fFileA(NULL),
103     fFileName(""),
104     fExtensions(NULL),
105     fFilters(NULL)
106 {
107 // Normal constructor.
108 }
109
110 //______________________________________________________________________________
111 AliAODHandler::~AliAODHandler() 
112 {
113  // Destructor.
114   if (fAODEvent) delete fAODEvent;
115   if(fFileA){
116     // is already handled in TerminateIO
117     fFileA->Close();
118     delete fFileA;
119     fTreeA = 0;
120   }
121   if (fTreeA) delete fTreeA;
122   if (fExtensions) {fExtensions->Delete(); delete fExtensions;}
123   if (fFilters)    {fFilters->Delete();    delete fFilters;}
124 }
125
126 //______________________________________________________________________________
127 Bool_t AliAODHandler::Init(Option_t* opt)
128 {
129   // Initialize IO
130   //
131   // Create the AODevent object
132     
133   Bool_t createStdAOD = fIsStandard || fFillAOD;
134   if(!fAODEvent && createStdAOD){
135     fAODEvent = new AliAODEvent();
136     if (fIsStandard) 
137       fAODEvent->CreateStdContent();
138   }
139   //
140   // File opening according to execution mode
141   TString option(opt);
142   option.ToLower();
143   if (createStdAOD) {
144     TDirectory *owd = gDirectory;
145     if (option.Contains("proof")) {
146       // proof
147       // Merging via files. Need to access analysis manager via interpreter.
148       gROOT->ProcessLine(Form("AliAnalysisDataContainer *c_common_out = AliAnalysisManager::GetAnalysisManager()->GetCommonOutputContainer();"));
149       gROOT->ProcessLine(Form("AliAnalysisManager::GetAnalysisManager()->OpenProofFile(c_common_out, \"RECREATE\");"));
150       fFileA = gFile;
151     } else {
152       // local and grid
153       fFileA = new TFile(fFileName.Data(), "RECREATE");
154     }
155     CreateTree(1);
156     owd->cd();
157   }  
158   if (fExtensions) {
159      TIter next(fExtensions);
160      AliAODExtension *ext;
161      while ((ext=(AliAODExtension*)next())) ext->Init(option);
162   }   
163   if (fFilters) {
164      TIter nextf(fFilters);
165      AliAODExtension *filteredAOD;
166      while ((filteredAOD=(AliAODExtension*)nextf())) {
167         filteredAOD->SetEvent(fAODEvent);
168         filteredAOD->Init(option);
169      }   
170   }   
171   return kTRUE;
172 }
173
174 //______________________________________________________________________________
175 void AliAODHandler::StoreMCParticles(){
176
177   // 
178   // Remap the labels from ESD stack and store
179   // the AODMCParticles, makes only sense if we have
180   // the mcparticles branch
181   // has to be done here since we cannot know in advance 
182   // which particles are needed (e.g. by the tracks etc.)
183   //
184   // Particles have been selected by AliMCEventhanlder->SelectParticle()
185   // To use the MCEventhandler here we need to set it from the outside
186   // can vanish when Handler go to the ANALYSISalice library
187   //
188   // The Branch booking for mcParticles and mcHeader has to happen 
189   // in an external task for now since the AODHandler does not have access
190   // the AnalysisManager. For the same reason the pointer t o the MCEventH
191   // has to passed to the AOD Handler by this task 
192   // (doing this in the steering macro would not work on PROOF)
193
194   if (!fAODEvent) return;
195   TClonesArray *mcarray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()); 
196   if(!mcarray)return;
197
198   AliAODMCHeader *mcHeader = (AliAODMCHeader*)fAODEvent->FindListObject(AliAODMCHeader::StdBranchName()); 
199   if(!mcHeader)return;
200
201   // Get the MC Infos.. Handler needs to be set before 
202   // while adding the branch
203   // This needs to be done, not to depend on the AnalysisManager
204
205   if(!fMCEventH)return;
206   if(!fMCEventH->MCEvent())return;
207   AliStack *pStack = fMCEventH->MCEvent()->Stack();
208   if(!pStack)return;
209
210   fMCEventH->CreateLabelMap();
211
212   //
213   // Get the Event Header 
214   // 
215
216   AliHeader* header = fMCEventH->MCEvent()->Header();
217    // get the MC vertex
218   AliGenEventHeader* genHeader = 0;
219   if (header) genHeader = header->GenEventHeader();
220   if (genHeader) {
221       TArrayF vtxMC(3);
222       genHeader->PrimaryVertex(vtxMC);
223       mcHeader->SetVertex(vtxMC[0],vtxMC[1],vtxMC[2]);
224
225       // we search the MCEventHeaders first 
226       // Two cases, cocktail or not...
227       AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
228       if(genCocktailHeader){
229           // we have a coktail header
230           mcHeader->AddGeneratorName(genHeader->GetName());
231           // Loop from the back so that the first one sets the process type
232           TList* headerList = genCocktailHeader->GetHeaders();
233           for(int i = headerList->GetEntries()-1;i>=0;--i){
234               AliGenEventHeader *headerEntry = dynamic_cast<AliGenEventHeader*>(headerList->At(i));
235               SetMCHeaderInfo(mcHeader,headerEntry);
236           }
237       }
238       else{
239           // No Cocktail just take the first one
240           SetMCHeaderInfo(mcHeader,genHeader);
241       }
242   }
243   
244
245
246
247
248   // Store the AliAODParticlesMC
249   AliMCEvent* mcEvent = fMCEventH->MCEvent();
250   
251   Int_t np    = mcEvent->GetNumberOfTracks();
252   Int_t nprim = mcEvent->GetNumberOfPrimaries();
253
254
255   Int_t j = 0;
256   TClonesArray& l = *mcarray;
257
258   for(int i = 0; i < np; ++i){
259       if(fMCEventH->IsParticleSelected(i)){
260           Int_t flag = 0;
261           AliMCParticle* mcpart =  (AliMCParticle*) mcEvent->GetTrack(i);
262           if(i<nprim)flag |= AliAODMCParticle::kPrimary;
263           
264           if(mcEvent->IsPhysicalPrimary(i))flag |= AliAODMCParticle::kPhysicalPrim;
265           
266           if(fMCEventH->GetNewLabel(i)!=j){
267               AliError(Form("MISMATCH New label %d j: %d",fMCEventH->GetNewLabel(i),j));
268           }
269
270           AliAODMCParticle mcpart_tmp(mcpart,i,flag);
271           
272           mcpart_tmp.SetStatus(mcpart->Particle()->GetStatusCode());
273           // 
274           Int_t d0 =  mcpart_tmp.GetDaughter(0);
275           Int_t d1 =  mcpart_tmp.GetDaughter(1);
276           Int_t m =   mcpart_tmp.GetMother();
277           
278           // other than for the track labels, negative values mean
279           // no daughter/mother so preserve it
280           
281           if(d0<0 && d1<0){
282               // no first daughter -> no second daughter
283               // nothing to be done
284               // second condition not needed just for sanity check at the end
285               mcpart_tmp.SetDaughter(0,d0);
286               mcpart_tmp.SetDaughter(1,d1);
287           } else if(d1 < 0 && d0 >= 0) {
288               // Only one daughter
289               // second condition not needed just for sanity check at the end
290               if(fMCEventH->IsParticleSelected(d0)){
291                   mcpart_tmp.SetDaughter(0,fMCEventH->GetNewLabel(d0));
292               } else {
293                   mcpart_tmp.SetDaughter(0,-1);
294               }
295               mcpart_tmp.SetDaughter(1,d1);
296           }
297           else if (d0 > 0 && d1 > 0 ){
298               // we have two or more daughters loop on the stack to see if they are
299               // selected
300               Int_t d0_tmp = -1;
301               Int_t d1_tmp = -1;
302               for(int id = d0; id<=d1;++id){
303                   if(fMCEventH->IsParticleSelected(id)){
304                       if(d0_tmp==-1){
305                           // first time
306                           d0_tmp = fMCEventH->GetNewLabel(id);
307                           d1_tmp = d0_tmp; // this is to have the same schema as on the stack i.e. with one daugther d0 and d1 are the same 
308                       }
309                       else d1_tmp = fMCEventH->GetNewLabel(id);
310                   }
311               }
312               mcpart_tmp.SetDaughter(0,d0_tmp);
313               mcpart_tmp.SetDaughter(1,d1_tmp);
314           } else {
315               AliError(Form("Unxpected indices %d %d",d0,d1));
316           }
317           
318           if(m<0){
319               mcpart_tmp.SetMother(m);
320           } else {
321               if(fMCEventH->IsParticleSelected(m))mcpart_tmp.SetMother(fMCEventH->GetNewLabel(m));
322               else AliError(Form("PROBLEM Mother not selected %d \n", m));
323           }
324
325           new (l[j++]) AliAODMCParticle(mcpart_tmp);
326           
327       }
328   }
329   AliInfo(Form("AliAODHandler::StoreMCParticles: Selected %d (Primaries %d / total %d) after validation \n",
330                j,nprim,np));
331   
332   // Set the labels in the AOD output...
333   // Remapping
334
335   // AODTracks
336   TClonesArray* tracks = fAODEvent->GetTracks();
337   if(tracks){
338     for(int it = 0; it < fAODEvent->GetNTracks();++it){
339       AliAODTrack *track = fAODEvent->GetTrack(it);
340       
341       Int_t sign = 1;
342       Int_t label = track->GetLabel();
343       if(label<0){ // preserve the sign for later usage
344         label *= -1;
345         sign  = -1;
346       }
347
348       if (label >= AliMCEvent::BgLabelOffset()) label =  mcEvent->BgLabelToIndex(label);
349       if(label > np || track->GetLabel() == 0){
350         AliWarning(Form("Wrong ESD track label %5d (%5d)",track->GetLabel(), label));
351       }
352       if(fMCEventH->GetNewLabel(label) == 0){
353         AliWarning(Form("New label not found for %5d (%5d)",track->GetLabel(), label));
354       }
355       track->SetLabel(sign*fMCEventH->GetNewLabel(label));
356     }
357   }
358   
359   // AOD calo cluster
360   TClonesArray *clusters = fAODEvent->GetCaloClusters();
361   if(clusters){
362     for (Int_t iClust = 0;iClust < fAODEvent->GetNumberOfCaloClusters(); ++iClust) {
363       AliAODCaloCluster * cluster = fAODEvent->GetCaloCluster(iClust);
364       UInt_t nLabel    = cluster->GetNLabels();
365       // Ugly but do not want to fragment memory by creating 
366       // new Int_t (nLabel)
367       Int_t* labels    = const_cast<Int_t*>(cluster->GetLabels());
368       if (labels){
369         for(UInt_t i = 0;i < nLabel;++i){
370           labels[i] = fMCEventH->GetNewLabel(cluster->GetLabelAt(i));
371         }
372       }
373       //      cluster->SetLabels(labels,nLabel);
374     }// iClust
375   }// clusters
376
377   // AOD tracklets
378   AliAODTracklets *tracklets = fAODEvent->GetTracklets();
379   if(tracklets){
380     for(int it = 0;it < tracklets->GetNumberOfTracklets();++it){
381       int label0 = tracklets->GetLabel(it,0);
382       int label1 = tracklets->GetLabel(it,1);
383       if(label0>=0)label0 = fMCEventH->GetNewLabel(label0);      
384       if(label1>=0)label1 = fMCEventH->GetNewLabel(label1);
385       tracklets->SetLabel(it,0,label0);
386       tracklets->SetLabel(it,1,label1);
387     }
388   }
389
390 }
391
392 //______________________________________________________________________________
393 Bool_t AliAODHandler::FinishEvent()
394 {
395   // Fill data structures
396   if(fFillAOD && fFillAODRun && fAODEvent){
397       fAODEvent->MakeEntriesReferencable();
398       fTreeA->BranchRef();
399       FillTree();
400   }
401
402   if (fFillAOD && fFillAODRun) {      
403       if (fExtensions) {
404           TIter next(fExtensions);
405           AliAODExtension *ext;
406           while ((ext=(AliAODExtension*)next())) ext->FinishEvent();
407       }
408       if (fFilters) {   
409           TIter nextf(fFilters);
410           AliAODExtension *ext;
411           while ((ext=(AliAODExtension*)nextf())) {
412               ext->FinishEvent();
413           }  
414       }       
415   }  
416   if (fIsStandard) fAODEvent->ResetStd();
417   if (fAODEvent) {
418     TClonesArray *mcarray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName()); 
419     if(mcarray) mcarray->Delete();
420     
421     AliAODMCHeader *mcHeader = (AliAODMCHeader*)fAODEvent->FindListObject(AliAODMCHeader::StdBranchName()); 
422     if(mcHeader) mcHeader->Reset();
423   }
424   // Reset AOD replication flag
425   fAODIsReplicated = kFALSE;
426   return kTRUE;
427 }
428
429 //______________________________________________________________________________
430 Bool_t AliAODHandler::Terminate()
431 {
432   // Terminate 
433   AddAODtoTreeUserInfo();
434   if (fExtensions) {
435     TIter next(fExtensions);
436     AliAODExtension *ext;
437     while ((ext=(AliAODExtension*)next())) ext->GetTree()->GetUserInfo()->Add(ext->GetAOD());
438   }  
439   return kTRUE;
440 }
441
442 //______________________________________________________________________________
443 Bool_t AliAODHandler::TerminateIO()
444 {
445   // Terminate IO
446   if (fFileA) {
447     fFileA->Write();
448     fFileA->Close();
449     delete fFileA;
450     fFileA = 0;
451     // When closing the file, the tree is also deleted.
452     fTreeA = 0;
453   }
454   if (fExtensions) {
455     TIter next(fExtensions);
456     AliAODExtension *ext;
457     while ((ext=(AliAODExtension*)next())) ext->TerminateIO();
458   }  
459   if (fFilters) {
460     TIter nextf(fFilters);
461     AliAODExtension *ext;
462     while ((ext=(AliAODExtension*)nextf())) ext->TerminateIO();
463   }  
464   return kTRUE;
465 }
466
467 //______________________________________________________________________________
468 void AliAODHandler::CreateTree(Int_t flag)
469 {
470     // Creates the AOD Tree
471     fTreeA = new TTree("aodTree", "AliAOD tree");
472     fTreeA->Branch(fAODEvent->GetList());
473     if (flag == 0) fTreeA->SetDirectory(0);
474 }
475
476 //______________________________________________________________________________
477 void AliAODHandler::FillTree()
478 {
479  
480     // Fill the AOD Tree
481     fTreeA->Fill();
482 }
483
484 //______________________________________________________________________________
485 void AliAODHandler::AddAODtoTreeUserInfo()
486 {
487   // Add aod event to tree user info
488   if (fTreeA) fTreeA->GetUserInfo()->Add(fAODEvent);
489   // Now the tree owns our fAODEvent...
490   fAODEvent = 0;
491 }
492
493 //______________________________________________________________________________
494 void AliAODHandler::AddBranch(const char* cname, void* addobj, const char* filename)
495 {
496     // Add a new branch to the aod. Added optional filename parameter if the
497     // branch should be written to a separate file.
498     if (strlen(filename)) {
499        AliAODExtension *ext = AddExtension(filename);
500        ext->AddBranch(cname, addobj);
501        return;
502     } else {
503        // Add branch to all filters
504       if (fFilters) {
505          TIter next(fFilters);
506          AliAODExtension *ext;
507          while ((ext=(AliAODExtension*)next())) ext->AddBranch(cname, addobj);
508       }
509     }
510     TDirectory *owd = gDirectory;
511     if (fFileA) {
512       fFileA->cd();
513     }
514     char** apointer = (char**) addobj;
515     TObject* obj = (TObject*) *apointer;
516
517     fAODEvent->AddObject(obj);
518  
519     const Int_t kSplitlevel = 99; // default value in TTree::Branch()
520     const Int_t kBufsize = 32000; // default value in TTree::Branch()
521
522     if (!fTreeA->FindBranch(obj->GetName())) {
523       // Do the same as if we book via 
524       // TTree::Branch(TCollection*)
525       
526       fTreeA->Bronch(obj->GetName(), cname, fAODEvent->GetList()->GetObjectRef(obj),
527                      kBufsize, kSplitlevel - 1);
528     }
529     owd->cd();
530 }
531
532 //______________________________________________________________________________
533 AliAODExtension *AliAODHandler::AddExtension(const char *filename, const char *title)
534 {
535 // Add an AOD extension with some branches in a different file.
536    TString fname(filename);
537    if (!fname.EndsWith(".root")) fname += ".root";
538    if (!fExtensions) {
539       fExtensions = new TObjArray();
540       fExtensions->SetOwner();
541    }   
542    AliAODExtension *ext = (AliAODExtension*)fExtensions->FindObject(fname);
543    if (!ext) {
544       ext = new AliAODExtension(fname, title);
545       fExtensions->Add(ext);
546    }   
547    return ext;
548 }
549
550 //______________________________________________________________________________
551 AliAODExtension *AliAODHandler::GetExtension(const char *filename) const
552 {
553 // Getter for AOD extensions via file name.
554    if (!fExtensions) return NULL;
555    return (AliAODExtension*)fExtensions->FindObject(filename);
556 }   
557
558 //______________________________________________________________________________
559 AliAODExtension *AliAODHandler::AddFilteredAOD(const char *filename, const char *filtername)
560 {
561 // Add an AOD extension that can write only AOD events that pass a user filter.
562    if (!fFilters) {
563       fFilters = new TObjArray();
564       fFilters->SetOwner();
565    } 
566    AliAODExtension *filter = (AliAODExtension*)fFilters->FindObject(filename);
567    if (!filter) {
568       filter = new AliAODExtension(filename, filtername, kTRUE);
569       fFilters->Add(filter);
570    }
571    return filter;
572 }      
573
574 //______________________________________________________________________________
575 AliAODExtension *AliAODHandler::GetFilteredAOD(const char *filename) const
576 {
577 // Getter for AOD filters via file name.
578    if (!fFilters) return NULL;
579    return (AliAODExtension*)fFilters->FindObject(filename);
580 }   
581    
582 //______________________________________________________________________________
583 void AliAODHandler::SetOutputFileName(const char* fname)
584 {
585 // Set file name.
586    fFileName = fname;
587 }
588
589 //______________________________________________________________________________
590 const char *AliAODHandler::GetOutputFileName()
591 {
592 // Get file name.
593    return fFileName.Data();
594 }
595
596 //______________________________________________________________________________
597 const char *AliAODHandler::GetExtraOutputs() const
598 {
599 // Get extra outputs as a string separated by commas.
600    static TString eoutputs;
601    eoutputs = "";
602    TObject *obj;
603    if (fExtensions) {
604       TIter next1(fExtensions);
605       while ((obj=next1())) {
606          if (!eoutputs.IsNull()) eoutputs += ",";
607          eoutputs += obj->GetName();
608       }
609    }
610    if (fFilters) {
611       TIter next2(fFilters);
612       while ((obj=next2())) {
613          if (!eoutputs.IsNull()) eoutputs += ",";
614          eoutputs += obj->GetName();
615       }
616    }
617    return eoutputs.Data();
618 }
619
620 //______________________________________________________________________________
621 void  AliAODHandler::SetMCHeaderInfo(AliAODMCHeader *mcHeader,AliGenEventHeader *genHeader){
622
623
624   // Utility function to cover different cases for the AliGenEventHeader
625   // Needed since different ProcessType and ImpactParamter are not 
626   // in the base class...
627   // We don't encode process types for event cocktails yet
628   // could be done e.g. by adding offsets depnding on the generator
629
630   mcHeader->AddGeneratorName(genHeader->GetName());
631   if(!genHeader)return;
632   AliGenPythiaEventHeader *pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
633     if (pythiaGenHeader) {
634       mcHeader->SetEventType(pythiaGenHeader->ProcessType());
635       mcHeader->SetPtHard(pythiaGenHeader->GetPtHard());
636       return;
637     }
638     
639     AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
640
641   if (dpmJetGenHeader){
642     mcHeader->SetEventType(dpmJetGenHeader->ProcessType());
643     return;
644   } 
645
646   AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
647   if(hijingGenHeader){
648     mcHeader->SetImpactParameter(hijingGenHeader->ImpactParameter());
649     return;
650   }
651   
652   AliWarning(Form("MC Eventheader not known: %s",genHeader->GetName()));
653
654 }
655
656 ClassImp(AliAODExtension)
657
658 //-------------------------------------------------------------------------
659 //     Support class for AOD extensions. This is created by the user analysis
660 //     that requires a separate file for some AOD branches. The name of the 
661 //     AliAODExtension object is the file name where the AOD branches will be
662 //     stored.
663 //-------------------------------------------------------------------------
664
665 //______________________________________________________________________________
666 AliAODExtension::AliAODExtension(const char* name, const char* title, Bool_t isfilter)
667                 :TNamed(name,title), 
668                  fAODEvent(0), 
669                  fTreeE(0), 
670                  fFileE(0), 
671                  fNtotal(0), 
672                  fNpassed(0),
673                  fSelected(kFALSE)
674 {
675 // Constructor.
676   if (isfilter) {
677     TObject::SetBit(kFilteredAOD);
678     printf("####### Added AOD filter %s\n", name);
679   } else printf("####### Added AOD extension %s\n", name);
680 }   
681
682 //______________________________________________________________________________
683 AliAODExtension::~AliAODExtension()
684 {
685 // Destructor.
686   if(fFileE){
687     // is already handled in TerminateIO
688     fFileE->Close();
689     delete fFileE;
690     fTreeE = 0;
691     fAODEvent = 0;
692   }
693   if (fTreeE) delete fTreeE;
694 }
695
696 //______________________________________________________________________________
697 void AliAODExtension::AddBranch(const char* cname, void* addobj)
698 {
699     // Add a new branch to the aod 
700 //    if (IsFilteredAOD()) {
701 //       Error("AddBranch", "Not allowed to add branched to filtered AOD's.");
702 //       return;
703 //    }   
704     if (!fAODEvent) {
705        char type[20];
706        gROOT->ProcessLine(Form("TString s_tmp; AliAnalysisManager::GetAnalysisManager()->GetAnalysisTypeString(s_tmp); sprintf((char*)%p, \"%%s\", s_tmp.Data());", type));
707        Init(type);
708     }
709     TDirectory *owd = gDirectory;
710     if (fFileE) {
711       fFileE->cd();
712     }
713     char** apointer = (char**) addobj;
714     TObject* obj = (TObject*) *apointer;
715
716     fAODEvent->AddObject(obj);
717  
718     const Int_t kSplitlevel = 99; // default value in TTree::Branch()
719     const Int_t kBufsize = 32000; // default value in TTree::Branch()
720
721     if (!fTreeE->FindBranch(obj->GetName())) {
722       // Do the same as if we book via 
723       // TTree::Branch(TCollection*)
724       fTreeE->Bronch(obj->GetName(), cname, fAODEvent->GetList()->GetObjectRef(obj),
725                      kBufsize, kSplitlevel - 1);
726     }
727     owd->cd();
728 }
729
730 //______________________________________________________________________________
731 Bool_t AliAODExtension::FinishEvent()
732 {
733 // Fill current event.
734   fNtotal++;
735   if (!IsFilteredAOD()) {
736     fAODEvent->MakeEntriesReferencable();
737     fTreeE->Fill();
738     return kTRUE;
739   }  
740   // Filtered AOD. Fill only if event is selected.
741   if (!fSelected) return kTRUE;
742   fNpassed++;
743   fTreeE->Fill();
744   fSelected = kFALSE; // so that next event will not be selected unless demanded
745   return kTRUE;
746 }  
747
748 //______________________________________________________________________________
749 Bool_t AliAODExtension::Init(Option_t *option)
750 {
751 // Initialize IO.
752   if(!fAODEvent) fAODEvent = new AliAODEvent();
753   TDirectory *owd = gDirectory;
754   TString opt(option);
755   opt.ToLower();
756   if (opt.Contains("proof")) {
757     // proof
758     // Merging via files. Need to access analysis manager via interpreter.
759     gROOT->ProcessLine(Form("AliAnalysisManager::GetAnalysisManager()->OpenProofFile(\"%s\", \"RECREATE\");", fName.Data()));
760     fFileE = gFile;
761   } else {
762     fFileE = new TFile(GetName(), "RECREATE");
763   }  
764   fTreeE = new TTree("aodTree", "AliAOD tree");
765   fTreeE->Branch(fAODEvent->GetList());
766   fTreeE->BranchRef();
767   owd->cd();
768   return kTRUE;
769 }
770
771 //______________________________________________________________________________
772 void AliAODExtension::SetEvent(AliAODEvent *event)
773 {
774 // Connects to an external event
775    if (!IsFilteredAOD()) {
776       Error("SetEvent", "Not allowed to set external event for filtered AOD's");   
777       return;
778    }
779    // Use the copy constructor or assignment operator to synchronize with external event.
780 //   AliAODEvent &other = *event;
781 //   if (!fAODEvent)     fAODEvent = new AliAODEvent(other);
782 //   else if (fSelected) *fAODEvent = other;
783    fAODEvent = event;
784 }
785    
786 //______________________________________________________________________________
787 Bool_t AliAODExtension::TerminateIO()
788 {
789   // Terminate IO
790   if (TObject::TestBit(kFilteredAOD))
791     printf("AOD Filter %s: events processed: %d   passed: %d\n", GetName(), fNtotal, fNpassed);
792   else
793     printf("AOD extension %s: events processed: %d\n", GetName(), fNtotal);
794   if (fFileE) {
795     fFileE->Write();
796     fFileE->Close();
797     delete fFileE;
798     fFileE = 0;
799     fTreeE = 0;
800     fAODEvent = 0;
801   }
802   return kTRUE;
803 }