STEER/AliMCEventHandler.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 /* $Id$ */
  17 //---------------------------------------------------------------------------------
  18 //                          Class AliMCEventHandler
  19 // This class gives access to MC truth during the analysis.
  20 // Monte Carlo truth is containe in the kinematics tree (produced particles) and
  21 // the tree of reference hits.
  22 //
  23 // Origin: Andreas Morsch, CERN, andreas.morsch@cern.ch
  24 //---------------------------------------------------------------------------------
  25
  26
  27
  28 #include "AliMCEventHandler.h"
  29 #include "AliTrackReference.h"
  30 #include "AliHeader.h"
  31 #include "AliStack.h"
  32 #include "AliAnalysisManager.h"
  33
  34 #include <TTree.h>
  35 #include <TFile.h>
  36 #include <TParticle.h>
  37 #include <TClonesArray.h>
  38 #include <TDirectoryFile.h>
  39 #include <TArrow.h>
  40 #include <TMarker.h>
  41 #include <TH2F.h>
  42
  43
  44 ClassImp(AliMCEventHandler)
  45
  46 AliMCEventHandler::AliMCEventHandler() :
  47     AliVirtualEventHandler(),
  48     fFileE(0),
  49     fFileK(0),
  50     fFileTR(0),
  51     fTmpFileTR(0),
  52     fTreeE(0),
  53     fTreeK(0),
  54     fTreeTR(0),
  55     fTmpTreeTR(0),
  56     fStack(0),
  57     fHeader(0),
  58     fTrackReferences(0),
  59     fNEvent(-1),
  60     fEvent(-1),
  61     fNprimaries(-1),
  62     fNparticles(-1),
  63     fPathName("./")
  64 {
  65     // Default constructor
  66 }
  67
  68 AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
  69     AliVirtualEventHandler(name, title),
  70     fFileE(0),
  71     fFileK(0),
  72     fFileTR(0),
  73     fTmpFileTR(0),
  74     fTreeE(0),
  75     fTreeK(0),
  76     fTreeTR(0),
  77     fTmpTreeTR(0),
  78     fStack(0),
  79     fHeader(new AliHeader()),
  80     fTrackReferences(new TClonesArray("AliTrackReference", 200)),
  81     fNEvent(-1),
  82     fEvent(-1),
  83     fNprimaries(-1),
  84     fNparticles(-1),
  85     fPathName("./")
  86 {
  87     // Constructor
  88 }
  89 AliMCEventHandler::~AliMCEventHandler()
  90 {
  91     // Destructor
  92     delete fFileE;
  93     delete fFileK;
  94     delete fFileTR;
  95 }
  96
  97 Bool_t AliMCEventHandler::InitIO(Option_t* /*opt*/)
  98 {
  99     // Initialize input
 100     //
 101     fFileE = new TFile(Form("%sgalice.root", fPathName));
 102     if (!fFileE) printf("galice.root not found in directory %s ! \n", fPathName);
 103
 104     fFileE->GetObject("TE", fTreeE);
 105     fTreeE->SetBranchAddress("Header", &fHeader);
 106     fNEvent = fTreeE->GetEntries();
 107     //
 108     // Tree K
 109     fFileK = new TFile(Form("%sKinematics.root", fPathName));
 110     if (!fFileK) printf("Kinematics.root not found in directory %s ! \n", fPathName);
 111     //
 112     // Tree TR
 113     fFileTR = new TFile(Form("%sTrackRefs.root", fPathName));
 114     if (!fFileTR) printf("TrackRefs.root not found in directory %s ! \n", fPathName);
 115     //
 116     // Reset the event number
 117     fEvent = -1;
 118     printf("Number of events in this directory %5d \n", fNEvent);
 119     return kTRUE;
 120
 121 }
 122
 123 Bool_t AliMCEventHandler::BeginEvent()
 124 {
 125     // Read the next event
 126     printf("AliMCEventHandler::BeginEvent called \n");
 127
 128     fEvent++;
 129     char folder[20];
 130     sprintf(folder, "Event%d", fEvent);
 131     // TreeE
 132     fTreeE->GetEntry(fEvent);
 133     fStack = fHeader->Stack();
 134     // Tree K
 135     TDirectoryFile* dirK  = 0;
 136     fFileK->GetObject(folder, dirK);
 137     dirK->GetObject("TreeK", fTreeK);
 138     fStack->ConnectTree(fTreeK);
 139     fStack->GetEvent();
 140     //Tree TR
 141     TDirectoryFile* dirTR = 0;
 142     fFileTR->GetObject(folder, dirTR);
 143     dirTR->GetObject("TreeTR", fTreeTR);
 144     if (fTreeTR->GetBranch("AliRun")) {
 145         // This is an old format with one branch per detector not in synch with TreeK
 146         ReorderAndExpandTreeTR();
 147     } else {
 148         // New format
 149         fTreeTR->SetBranchAddress("TrackReferences", &fTrackReferences);
 150     }
 151
 152     //
 153     fNparticles = fStack->GetNtrack();
 154     fNprimaries = fStack->GetNprimary();
 155     printf("Event#: %5d Number of particles %5d \n", fEvent, fNparticles);
 156
 157     return kTRUE;
 158 }
 159
 160 Int_t AliMCEventHandler::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
 161 {
 162     // Retrieve entry i
 163     if (i > -1 && i < fNparticles) {
 164         fTreeTR->GetEntry(fStack->TreeKEntry(i));
 165     } else {
 166         printf("AliMCEventHandler::GetEntry: Index out of range");
 167     }
 168     particle = fStack->Particle(i);
 169     trefs    = fTrackReferences;
 170     printf("Returning %5d \n", particle->GetPdgCode());
 171
 172     return trefs->GetEntries();
 173
 174 }
 175
 176 void AliMCEventHandler::DrawCheck(Int_t i)
 177 {
 178     // Retrieve entry i and draw momentum vector and hits
 179     if (i > -1 && i < fNparticles) {
 180         fTreeTR->GetEntry(fStack->TreeKEntry(i));
 181     } else {
 182         printf("AliMCEventHandler::GetEntry: Index out of range");
 183     }
 184
 185     TParticle* particle = fStack->Particle(i);
 186     Int_t nh =  fTrackReferences->GetEntries();
 187
 188     TH2F*    h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
 189     Float_t x0 = particle->Vx();
 190     Float_t y0 = particle->Vy();
 191
 192     Float_t x1 = particle->Vx() + particle->Px() * 50.;
 193     Float_t y1 = particle->Vy() + particle->Py() * 50.;
 194
 195     TArrow*  a = new TArrow(x0, y0, x1, y1, 0.01);
 196     h->Draw();
 197     a->SetLineColor(2);
 198
 199     a->Draw();
 200
 201     for (Int_t ih = 0; ih < nh; ih++) {
 202         AliTrackReference* ref = (AliTrackReference*) fTrackReferences->At(ih);
 203         TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
 204         m->Draw();
 205         m->SetMarkerSize(0.4);
 206
 207     }
 208 }
 209
 210 Bool_t AliMCEventHandler::Notify()
 211 {
 212     // Notify about directory change
 213     TTree* tree = AliAnalysisManager::GetAnalysisManager()->GetTree();
 214     TFile *curfile = tree->GetCurrentFile();
 215     TString fileName(curfile->GetName());
 216     fileName.ReplaceAll("AliESDs.root", "");
 217     printf("AliMCEventHandler::Notify() file: %s\n", fileName.Data());
 218     fPathName = Form("%s",  fileName.Data());
 219     ResetIO();
 220     InitIO("");
 221     return kTRUE;
 222 }
 223
 224 void AliMCEventHandler::ResetIO()
 225 {
 226     if (fFileE)  delete fFileE;
 227     if (fFileK)  delete fFileK;
 228     if (fFileTR) delete fFileTR;
 229 }
 230
 231
 232 Bool_t AliMCEventHandler::FinishEvent()
 233 {
 234     // Dummy
 235     return kTRUE;
 236 }
 237
 238 Bool_t AliMCEventHandler::Terminate()
 239 {
 240     // Dummy
 241     return kTRUE;
 242 }
 243
 244 Bool_t AliMCEventHandler::TerminateIO()
 245 {
 246     // Dummy
 247     return kTRUE;
 248 }
 249
 250 void AliMCEventHandler::ReorderAndExpandTreeTR()
 251 {
 252 //
 253 //  Reorder and expand the track reference tree in order to match the kinematics tree.
 254 //  Copy the information from different branches into one
 255 //
 256 //  TreeTR
 257     if (fTmpTreeTR) delete fTmpTreeTR;
 258     if (fTmpFileTR) {
 259         fTmpFileTR->Close();
 260         delete fTmpFileTR;
 261     }
 262
 263     fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
 264     fTmpTreeTR = new TTree("TreeTR", "Track References");
 265     if (!fTrackReferences)  fTrackReferences = new TClonesArray("AliTrackReference", 100);
 266     fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 4000);
 267 //
 268     TClonesArray* trefs[7];
 269     for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
 270     if (fTreeTR){
 271         // make branch for central track references
 272         if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
 273         if (fTreeTR->GetBranch("ITS"))    fTreeTR->SetBranchAddress("ITS",    &trefs[1]);
 274         if (fTreeTR->GetBranch("TPC"))    fTreeTR->SetBranchAddress("TPC",    &trefs[2]);
 275         if (fTreeTR->GetBranch("TRD"))    fTreeTR->SetBranchAddress("TRD",    &trefs[3]);
 276         if (fTreeTR->GetBranch("TOF"))    fTreeTR->SetBranchAddress("TOF",    &trefs[4]);
 277         if (fTreeTR->GetBranch("FRAME"))  fTreeTR->SetBranchAddress("FRAME",  &trefs[5]);
 278         if (fTreeTR->GetBranch("MUON"))   fTreeTR->SetBranchAddress("MUON",   &trefs[6]);
 279     }
 280
 281     Int_t np = fStack->GetNprimary();
 282     Int_t nt = fTreeTR->GetEntries();
 283     //
 284     // Loop over tracks and find the secondaries with the help of the kine tree
 285     Int_t ifills = 0;
 286     Int_t it     = 0;
 287     Int_t itlast = 0;
 288
 289     for (Int_t ip = np - 1; ip > -1; ip--) {
 290         TParticle *part = fStack->Particle(ip);
 291 //      printf("Particle %5d %5d %5d %5d %5d %5d \n",
 292 //             ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
 293 //             part->GetLastDaughter(), part->TestBit(kTransportBit));
 294
 295         // Skip primaries that have not been transported
 296         Int_t dau1  = part->GetFirstDaughter();
 297         Int_t dau2  = -1;
 298         // Determine range of secondaries produced by this primary
 299         if (dau1 > -1) {
 300             Int_t inext = ip - 1;
 301             while (dau2 < 0) {
 302                 if (inext >= 0) {
 303                     part = fStack->Particle(inext);
 304                     dau2 =  part->GetFirstDaughter();
 305                     if (dau2 == -1 || dau2 < np) {
 306                         dau2 = -1;
 307                     } else {
 308                         dau2--;
 309                     }
 310                 } else {
 311                     dau2 = fStack->GetNtrack() - 1;
 312                 }
 313                 inext--;
 314             } // find upper bound
 315         }  // dau2 < 0
 316
 317 //      printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
 318
 319         //
 320         // Loop over reference hits and find secondary label
 321         Bool_t hasHits   = kFALSE;
 322         Bool_t isOutside = kFALSE;
 323         it = itlast;
 324         if (dau1 < np) continue;
 325
 326         while (!hasHits && !isOutside && it < nt) {
 327             fTreeTR->GetEntry(it++);
 328             for (Int_t ib = 0; ib < 7; ib++) {
 329                 if (!trefs[ib]) continue;
 330                 Int_t nh = trefs[ib]->GetEntries();
 331                 for (Int_t ih = 0; ih < nh; ih++) {
 332                     AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
 333                     Int_t label = tr->Label();
 334                     if (label >= dau1 && label <= dau2) {
 335                         hasHits = kTRUE;
 336                         itlast = it - 1;
 337                     }
 338                     if (label > dau2 || label < ip) {
 339                         isOutside = kTRUE;
 340                         itlast = it - 1;
 341                     }
 342                 } // hits
 343             } // branches
 344         } // entries
 345
 346         if (!hasHits) {
 347             for (Int_t id = dau1; (id <= dau2); id++) {
 348                 fTmpTreeTR->Fill();
 349                 ifills++;
 350             }
 351         } else {
 352             fTreeTR->GetEntry(itlast);
 353             for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
 354                 for (Int_t ib = 0; ib < 7; ib++) {
 355                     if (!trefs[ib]) continue;
 356                     Int_t nh = trefs[ib]->GetEntries();
 357                     for (Int_t ih = 0; ih < nh; ih++) {
 358                         AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
 359                         Int_t label = tr->Label();
 360                         // Skip primaries
 361                         if (label == ip) continue;
 362                         if (label > dau2 || label < dau1)
 363                             printf("Track Reference Label out of range !: %5d %5d %5d %5d \n", itlast, label, dau1, dau2);
 364                         if (label == id) {
 365                             // secondary found
 366                             tr->SetDetectorId(ib-1);
 367                             Int_t nref =  fTrackReferences->GetEntriesFast();
 368                             TClonesArray &lref = *fTrackReferences;
 369                             new(lref[nref]) AliTrackReference(*tr);
 370                         }
 371                     } // hits
 372                 } // branches
 373                 fTmpTreeTR->Fill();
 374                 fTrackReferences->Clear();
 375                 ifills++;
 376             } // daughters
 377         } // has hits
 378     } // tracks
 379     printf("Secondaries %5d %5d \n", ifills, fStack->GetNtrack() - fStack->GetNprimary());
 380     //
 381     // Now loop again and write the primaries
 382     it = nt - 1;
 383     for (Int_t ip = 0; ip < np; ip++) {
 384         Int_t labmax = -1;
 385         while (labmax < ip && it > -1) {
 386             fTreeTR->GetEntry(it--);
 387             for (Int_t ib = 0; ib < 7; ib++) {
 388                 if (!trefs[ib]) continue;
 389                 Int_t nh = trefs[ib]->GetEntries();
 390                 //
 391                 // Loop over reference hits and find primary labels
 392                 for (Int_t ih = 0; ih < nh; ih++) {
 393                     AliTrackReference* tr = (AliTrackReference*)  trefs[ib]->At(ih);
 394                     Int_t label = tr->Label();
 395                     if (label < np && label > labmax) {
 396                         labmax = label;
 397                     }
 398
 399                     if (label == ip) {
 400                         tr->SetDetectorId(ib-1);
 401                         Int_t nref = fTrackReferences->GetEntriesFast();
 402                         TClonesArray &lref = *fTrackReferences;
 403                         new(lref[nref]) AliTrackReference(*tr);
 404                     }
 405                 } // hits
 406             } // branches
 407         } // entries
 408         it++;
 409         fTmpTreeTR->Fill();
 410         fTrackReferences->Clear();
 411         ifills++;
 412     } // tracks
 413     // Check
 414     printf("All %5d %5d \n", ifills, fStack->GetNtrack());
 415     if (ifills != fStack->GetNtrack())
 416         printf("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, fStack->GetNtrack());
 417     fTreeTR = fTmpTreeTR;
 418
 419 }