STEER/AliAODHandler.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, 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 /* $Id$ */
  17
  18 //-------------------------------------------------------------------------
  19 //     Implementation of the Virtual Event Handler Interface for AOD
  20 //     Author: Andreas Morsch, CERN
  21 //-------------------------------------------------------------------------
  22
  23
  24 #include <TTree.h>
  25 #include <TFile.h>
  26 #include <TString.h>
  27 #include <TList.h>
  28
  29 #include "AliLog.h"
  30 #include "AliAODHandler.h"
  31 #include "AliAODEvent.h"
  32 #include "AliAODTracklets.h"
  33 #include "AliStack.h"
  34 #include "AliAODMCParticle.h"
  35 #include "AliAODMCHeader.h"
  36 #include "AliMCEventHandler.h"
  37 #include "AliMCEvent.h"
  38 #include "AliGenEventHeader.h"
  39 #include "AliGenHijingEventHeader.h"
  40 #include "AliGenDPMjetEventHeader.h"
  41 #include "AliGenPythiaEventHeader.h"
  42 #include "AliGenCocktailEventHeader.h"
  43
  44
  45
  46 ClassImp(AliAODHandler)
  47
  48 //______________________________________________________________________________
  49 AliAODHandler::AliAODHandler() :
  50     AliVEventHandler(),
  51     fIsStandard(kTRUE),
  52     fFillAOD(kTRUE),
  53     fNeedsHeaderReplication(kFALSE),
  54     fNeedsTracksBranchReplication(kFALSE),
  55     fNeedsVerticesBranchReplication(kFALSE),
  56     fNeedsV0sBranchReplication(kFALSE),
  57     fNeedsTrackletsBranchReplication(kFALSE),
  58     fNeedsPMDClustersBranchReplication(kFALSE),
  59     fNeedsJetsBranchReplication(kFALSE),
  60     fNeedsFMDClustersBranchReplication(kFALSE),
  61     fNeedsCaloClustersBranchReplication(kFALSE),
  62     fAODIsReplicated(kFALSE),
  63     fAODEvent(NULL),
  64     fMCEventH(NULL),
  65     fTreeA(NULL),
  66     fFileA(NULL),
  67     fFileName("")
  68 {
  69   // default constructor
  70 }
  71
  72 //______________________________________________________________________________
  73 AliAODHandler::AliAODHandler(const char* name, const char* title):
  74     AliVEventHandler(name, title),
  75     fIsStandard(kTRUE),
  76     fFillAOD(kTRUE),
  77     fNeedsHeaderReplication(kFALSE),
  78     fNeedsTracksBranchReplication(kFALSE),
  79     fNeedsVerticesBranchReplication(kFALSE),
  80     fNeedsV0sBranchReplication(kFALSE),
  81     fNeedsTrackletsBranchReplication(kFALSE),
  82     fNeedsPMDClustersBranchReplication(kFALSE),
  83     fNeedsJetsBranchReplication(kFALSE),
  84     fNeedsFMDClustersBranchReplication(kFALSE),
  85     fNeedsCaloClustersBranchReplication(kFALSE),
  86     fAODIsReplicated(kFALSE),
  87     fAODEvent(NULL),
  88     fMCEventH(NULL),
  89     fTreeA(NULL),
  90     fFileA(NULL),
  91     fFileName("")
  92 {
  93 }
  94
  95 //______________________________________________________________________________
  96 AliAODHandler::~AliAODHandler()
  97 {
  98   delete fAODEvent;
  99   if(fFileA){
 100     // is already handled in TerminateIO
 101     fFileA->Close();
 102     delete fFileA;
 103   }
 104   delete fTreeA;
 105  // destructor
 106 }
 107
 108 //______________________________________________________________________________
 109 Bool_t AliAODHandler::Init(Option_t* opt)
 110 {
 111   // Initialize IO
 112   //
 113   // Create the AODevent object
 114   if(!fAODEvent){
 115     fAODEvent = new AliAODEvent();
 116     if (fIsStandard) fAODEvent->CreateStdContent();
 117   }
 118   //
 119   // File opening according to execution mode
 120   TString option(opt);
 121   option.ToLower();
 122   if (option.Contains("proof")) {
 123     // proof
 124     if (option.Contains("special")) {
 125        // File for tree already opened on slave -> merging via files
 126        fFileA = gFile;
 127        CreateTree(1);
 128     } else {
 129        // Merging in memory
 130        CreateTree(0);
 131     }
 132   } else {
 133     // local and grid
 134     TDirectory *owd = gDirectory;
 135     fFileA = new TFile(fFileName.Data(), "RECREATE");
 136     CreateTree(1);
 137     owd->cd();
 138   }
 139   return kTRUE;
 140 }
 141
 142
 143 void AliAODHandler::StoreMCParticles(){
 144
 145   //
 146   // Remap the labels from ESD stack and store
 147   // the AODMCParticles, makes only sense if we have
 148   // the mcparticles branch
 149   // has to be done here since we cannot know in advance
 150   // which particles are needed (e.g. by the tracks etc.)
 151   //
 152   // Particles have been selected by AliMCEventhanlder->SelectParticle()
 153   // To use the MCEventhandler here we need to set it from the outside
 154   // can vanish when Handler go to the ANALYSISalice library
 155   //
 156   // The Branch booking for mcParticles and mcHeader has to happen
 157   // in an external task for now since the AODHandler does not have access
 158   // the AnalysisManager. For the same reason the pointer t o the MCEventH
 159   // has to passed to the AOD Handler by this task
 160   // (doing this in the steering macro would not work on PROOF)
 161
 162   TClonesArray *mcarray = (TClonesArray*)fAODEvent->FindListObject(AliAODMCParticle::StdBranchName());
 163   if(!mcarray)return;
 164   mcarray->Delete();
 165
 166   AliAODMCHeader *mcHeader = (AliAODMCHeader*)fAODEvent->FindListObject(AliAODMCHeader::StdBranchName());
 167   if(!mcHeader)return;
 168
 169   mcHeader->Reset();
 170
 171   // Get the MC Infos.. Handler needs to be set before
 172   // while adding the branch
 173   // This needs to be done, not to depend on the AnalysisManager
 174
 175   if(!fMCEventH)return;
 176   if(!fMCEventH->MCEvent())return;
 177   AliStack *pStack = fMCEventH->MCEvent()->Stack();
 178   if(!pStack)return;
 179
 180   fMCEventH->CreateLabelMap();
 181
 182   //
 183   // Get the Event Header
 184   //
 185
 186   AliHeader* header = fMCEventH->MCEvent()->Header();
 187   if (!header)return;
 188
 189   // get the MC vertex
 190   AliGenEventHeader* genHeader = header->GenEventHeader();
 191   TArrayF vtxMC(3);
 192   genHeader->PrimaryVertex(vtxMC);
 193   mcHeader->SetVertex(vtxMC[0],vtxMC[1],vtxMC[2]);
 194
 195   // we search the MCEventHeaders first
 196   // Two cases, cocktail or not...
 197   AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(genHeader);
 198   if(genCocktailHeader){
 199     // we have a coktail header
 200     mcHeader->AddGeneratorName(genHeader->GetName());
 201     // Loop from the back so that the first one sets the process type
 202     TList* headerList = genCocktailHeader->GetHeaders();
 203     for(int i = headerList->GetEntries()-1;i>=0;--i){
 204       AliGenEventHeader *headerEntry = dynamic_cast<AliGenEventHeader*>(headerList->At(i));
 205       SetMCHeaderInfo(mcHeader,headerEntry);
 206     }
 207   }
 208   else{
 209     // No Cocktail just take the first one
 210     SetMCHeaderInfo(mcHeader,genHeader);
 211   }
 212
 213
 214
 215
 216   // Store the AliAODParticlesMC
 217
 218   Int_t np    = pStack->GetNtrack();
 219   Int_t nprim = pStack->GetNprimary();
 220
 221
 222   Int_t j = 0;
 223   TClonesArray& l = *mcarray;
 224
 225   for(int i = 0;i < np;++i){
 226     if(fMCEventH->IsParticleSelected(i)){
 227
 228       Int_t flag = 0;
 229       TParticle *part = pStack->Particle(i);
 230       if(i<nprim)flag |= AliAODMCParticle::kPrimary;
 231       if(pStack->IsPhysicalPrimary(i))flag |= AliAODMCParticle::kPhysicalPrim;
 232
 233       if(fMCEventH->GetNewLabel(i)!=j){
 234         AliError(Form("MISMATCH New label %d j: %d",fMCEventH->GetNewLabel(i),j));
 235       }
 236       AliAODMCParticle mcpart_tmp(part,i,flag);
 237
 238       //
 239       Int_t d0 =  mcpart_tmp.GetDaughter(0);
 240       Int_t d1 =  mcpart_tmp.GetDaughter(1);
 241       Int_t m =  mcpart_tmp.GetMother();
 242
 243       // other than for the track labels, negative values mean
 244       // no daughter/mother so preserve it
 245
 246       if(d0<0 && d1<0){
 247         // no first daughter -> no second daughter
 248         // nothing to be done
 249         // second condition not needed just for sanity check at the end
 250         mcpart_tmp.SetDaughter(0,d0);
 251         mcpart_tmp.SetDaughter(1,d1);
 252       }
 253       else if(d1 < 0 && d0 >= 0){
 254         // Only one daughter
 255         // second condition not needed just for sanity check at the end
 256         if(fMCEventH->IsParticleSelected(d0)){
 257           mcpart_tmp.SetDaughter(0,fMCEventH->GetNewLabel(d0));
 258         }
 259         else{
 260           mcpart_tmp.SetDaughter(0,-1);
 261         }
 262         mcpart_tmp.SetDaughter(1,d1);
 263       }
 264       else if (d0 > 0 && d1 > 0 ){
 265         // we have two or more daughters loop on the stack to see if they are
 266         // selected
 267         Int_t d0_tmp = -1;
 268         Int_t d1_tmp = -1;
 269         for(int id = d0; id<=d1;++id){
 270           if(fMCEventH->IsParticleSelected(id)){
 271             if(d0_tmp==-1){
 272               // first time
 273               d0_tmp = fMCEventH->GetNewLabel(id);
 274               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
 275             }
 276             else d1_tmp = fMCEventH->GetNewLabel(id);
 277           }
 278         }
 279         mcpart_tmp.SetDaughter(0,d0_tmp);
 280         mcpart_tmp.SetDaughter(1,d1_tmp);
 281       }
 282       else{
 283         AliError(Form("Unxpected indices %d %d",d0,d1));
 284       }
 285
 286       if(m<0){
 287         mcpart_tmp.SetMother(m);
 288       }
 289       else{
 290         if(fMCEventH->IsParticleSelected(m))mcpart_tmp.SetMother(fMCEventH->GetNewLabel(m));
 291         else AliError("PROBLEM Mother not selected");
 292       }
 293
 294       new (l[j++]) AliAODMCParticle(mcpart_tmp);
 295
 296     }
 297   }
 298   AliInfo(Form("AliAODHandler::StoreMCParticles: Selected %d (Primaries %d / total %d) after validation",
 299                j,nprim,np));
 300
 301   // Set the labels in the AOD output...
 302   // Remapping
 303
 304   // AODTracks
 305   TClonesArray* tracks = fAODEvent->GetTracks();
 306   if(tracks){
 307     for(int it = 0; it < fAODEvent->GetNTracks();++it){
 308       AliAODTrack *track = fAODEvent->GetTrack(it);
 309
 310       if(TMath::Abs(track->GetLabel())>np||track->GetLabel()==0){
 311         AliWarning(Form("Wrong ESD track label %d",track->GetLabel()));
 312       }
 313       if(fMCEventH->GetNewLabel(track->GetLabel())==0){
 314         AliWarning(Form("New label not found for %d",track->GetLabel()));
 315       }
 316       track->SetLabel(fMCEventH->GetNewLabel(track->GetLabel()));
 317     }
 318   }
 319
 320   // AOD calo cluster
 321   TClonesArray *clusters = fAODEvent->GetCaloClusters();
 322   if(clusters){
 323     for (Int_t iClust = 0;iClust < fAODEvent->GetNCaloClusters(); ++iClust) {
 324       AliAODCaloCluster * cluster = fAODEvent->GetCaloCluster(iClust);
 325       UInt_t nLabel    = cluster->GetNLabel();
 326       // Ugly but do not want to fragment memory by creating
 327       // new Int_t (nLabel)
 328       Int_t* labels    = const_cast<Int_t*>(cluster->GetLabels());
 329       if (labels){
 330         for(UInt_t i = 0;i < nLabel;++i){
 331           labels[i] = fMCEventH->GetNewLabel(cluster->GetLabel(i));
 332         }
 333       }
 334       //      cluster->SetLabels(labels,nLabel);
 335     }// iClust
 336   }// clusters
 337
 338   // AOD tracklets
 339   AliAODTracklets *tracklets = fAODEvent->GetTracklets();
 340   if(tracklets){
 341     for(int it = 0;it < tracklets->GetNumberOfTracklets();++it){
 342       int label0 = tracklets->GetLabel(it,0);
 343       int label1 = tracklets->GetLabel(it,1);
 344       if(label0>=0)label0 = fMCEventH->GetNewLabel(label0);
 345       if(label1>=0)label1 = fMCEventH->GetNewLabel(label1);
 346       tracklets->SetLabel(it,0,label0);
 347       tracklets->SetLabel(it,1,label1);
 348     }
 349   }
 350
 351 }
 352
 353 Bool_t AliAODHandler::FinishEvent()
 354 {
 355   // Fill data structures
 356   if(fFillAOD){
 357     fAODEvent->MakeEntriesReferencable();
 358     StoreMCParticles();
 359     FillTree();
 360   }
 361
 362   if (fIsStandard) fAODEvent->ResetStd();
 363   // Reset AOD replication flag
 364   fAODIsReplicated = kFALSE;
 365   return kTRUE;
 366 }
 367
 368 //______________________________________________________________________________
 369 Bool_t AliAODHandler::Terminate()
 370 {
 371     // Terminate
 372     AddAODtoTreeUserInfo();
 373     return kTRUE;
 374 }
 375
 376 //______________________________________________________________________________
 377 Bool_t AliAODHandler::TerminateIO()
 378 {
 379     // Terminate IO
 380     if (fFileA) {
 381         fFileA->Close();
 382         delete fFileA;
 383     }
 384     return kTRUE;
 385 }
 386
 387 //______________________________________________________________________________
 388 void AliAODHandler::CreateTree(Int_t flag)
 389 {
 390     // Creates the AOD Tree
 391     fTreeA = new TTree("aodTree", "AliAOD tree");
 392     fTreeA->Branch(fAODEvent->GetList());
 393     if (flag == 0) fTreeA->SetDirectory(0);
 394 }
 395
 396 //______________________________________________________________________________
 397 void AliAODHandler::FillTree()
 398 {
 399     // Fill the AOD Tree
 400   fTreeA->Fill();
 401 }
 402
 403 //______________________________________________________________________________
 404 void AliAODHandler::AddAODtoTreeUserInfo()
 405 {
 406   // Add aod event to tree user info
 407   fTreeA->GetUserInfo()->Add(fAODEvent);
 408 }
 409
 410 //______________________________________________________________________________
 411 void AliAODHandler::AddBranch(const char* cname, void* addobj)
 412 {
 413     // Add a new branch to the aod
 414     TDirectory *owd = gDirectory;
 415     if (fFileA) {
 416       fFileA->cd();
 417     }
 418     char** apointer = (char**) addobj;
 419     TObject* obj = (TObject*) *apointer;
 420
 421     fAODEvent->AddObject(obj);
 422
 423     const Int_t kSplitlevel = 99; // default value in TTree::Branch()
 424     const Int_t kBufsize = 32000; // default value in TTree::Branch()
 425
 426     if (!fTreeA->FindBranch(obj->GetName())) {
 427       // Do the same as if we book via
 428       // TTree::Branch(TCollection*)
 429
 430       fTreeA->Bronch(obj->GetName(), cname, fAODEvent->GetList()->GetObjectRef(obj),
 431                      kBufsize, kSplitlevel - 1);
 432       //    fTreeA->Branch(obj->GetName(), cname, addobj);
 433     }
 434     owd->cd();
 435 }
 436
 437 //______________________________________________________________________________
 438 void AliAODHandler::SetOutputFileName(const char* fname)
 439 {
 440 // Set file name.
 441    fFileName = fname;
 442 }
 443
 444 //______________________________________________________________________________
 445 const char *AliAODHandler::GetOutputFileName()
 446 {
 447 // Get file name.
 448    return fFileName.Data();
 449 }
 450
 451 void  AliAODHandler::SetMCHeaderInfo(AliAODMCHeader *mcHeader,AliGenEventHeader *genHeader){
 452
 453
 454   // Utility function to cover different cases for the AliGenEventHeader
 455   // Needed since different ProcessType and ImpactParamter are not
 456   // in the base class...
 457   // We don't encode process types for event cocktails yet
 458   // coul be done e.g. by adding offsets depnding on the generator
 459
 460
 461   mcHeader->AddGeneratorName(genHeader->GetName());
 462
 463
 464
 465   if(!genHeader)return;
 466   AliGenPythiaEventHeader *pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
 467     if (pythiaGenHeader) {
 468       mcHeader->SetEventType(pythiaGenHeader->ProcessType());
 469       return;
 470     }
 471
 472     AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(genHeader);
 473
 474   if (dpmJetGenHeader){
 475     mcHeader->SetEventType(dpmJetGenHeader->ProcessType());
 476     return;
 477   }
 478
 479   AliGenHijingEventHeader* hijingGenHeader = dynamic_cast<AliGenHijingEventHeader*>(genHeader);
 480   if(hijingGenHeader){
 481     mcHeader->SetImpactParameter(hijingGenHeader->ImpactParameter());
 482     return;
 483   }
 484
 485   AliWarning(Form("MC Eventheader not known: %s",genHeader->GetName()));
 486
 487 }