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 "AliLog.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     AliVEventHandler(),
  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     fExtension(""),
  65     fFileNumber(0),
  66     fEventsPerFile(0)
  67 {
  68     // Default constructor
  69 }
  70
  71 AliMCEventHandler::AliMCEventHandler(const char* name, const char* title) :
  72     AliVEventHandler(name, title),
  73     fFileE(0),
  74     fFileK(0),
  75     fFileTR(0),
  76     fTmpFileTR(0),
  77     fTreeE(0),
  78     fTreeK(0),
  79     fTreeTR(0),
  80     fTmpTreeTR(0),
  81     fStack(0),
  82     fHeader(new AliHeader()),
  83     fTrackReferences(new TClonesArray("AliTrackReference", 200)),
  84     fNEvent(-1),
  85     fEvent(-1),
  86     fNprimaries(-1),
  87     fNparticles(-1),
  88     fPathName("./"),
  89     fExtension(""),
  90     fFileNumber(0),
  91     fEventsPerFile(0)
  92 {
  93     // Constructor
  94 }
  95 AliMCEventHandler::~AliMCEventHandler()
  96 {
  97     // Destructor
  98     delete fFileE;
  99     delete fFileK;
 100     delete fFileTR;
 101 }
 102
 103 Bool_t AliMCEventHandler::InitIO(Option_t* /*opt*/)
 104 {
 105     // Initialize input
 106     //
 107     fFileE = TFile::Open(Form("%sgalice.root", fPathName));
 108     if (!fFileE) AliFatal(Form("AliMCEventHandler:galice.root not found in directory %s ! \n", fPathName));
 109
 110     fFileE->GetObject("TE", fTreeE);
 111     fTreeE->SetBranchAddress("Header", &fHeader);
 112     fNEvent = fTreeE->GetEntries();
 113     //
 114     // Tree K
 115     fFileK = TFile::Open(Form("%sKinematics%s.root", fPathName, fExtension));
 116     if (!fFileK) AliFatal(Form("AliMCEventHandler:Kinematics.root not found in directory %s ! \n", fPathName));
 117     fEventsPerFile = fFileK->GetNkeys() - fFileK->GetNProcessIDs();
 118     //
 119     // Tree TR
 120     fFileTR = TFile::Open(Form("%sTrackRefs%s.root", fPathName, fExtension));
 121     if (!fFileTR) AliWarning(Form("AliMCEventHandler:TrackRefs.root not found in directory %s ! \n", fPathName));
 122     //
 123     // Reset the event number
 124     fEvent      = -1;
 125     fFileNumber =  0;
 126
 127     AliInfo(Form("AliMCEventHandler:Number of events in this directory %5d \n", fNEvent));
 128     return kTRUE;
 129
 130 }
 131
 132 Bool_t AliMCEventHandler::GetEvent(Int_t iev)
 133 {
 134     // Load the event number iev
 135     //
 136     // Calculate the file number
 137     Int_t inew  = iev/fEventsPerFile;
 138     if (inew != fFileNumber) {
 139         fFileNumber = inew;
 140         if (!OpenFile(fFileNumber)){
 141             return kFALSE;
 142         }
 143     }
 144     // Folder name
 145     char folder[20];
 146     sprintf(folder, "Event%d", iev);
 147
 148     // TreeE
 149     fTreeE->GetEntry(iev);
 150     fStack = fHeader->Stack();
 151     // Tree K
 152     TDirectoryFile* dirK  = 0;
 153     fFileK->GetObject(folder, dirK);
 154     if (!dirK) {
 155         AliWarning(Form("AliMCEventHandler: Event #%5d not found\n", iev));
 156         return kFALSE;
 157     }
 158     dirK->GetObject("TreeK", fTreeK);
 159     fStack->ConnectTree(fTreeK);
 160     fStack->GetEvent();
 161     //Tree TR
 162     if (fFileTR) {
 163         TDirectoryFile* dirTR = 0;
 164         fFileTR->GetObject(folder, dirTR);
 165         dirTR->GetObject("TreeTR", fTreeTR);
 166         if (fTreeTR->GetBranch("AliRun")) {
 167             // This is an old format with one branch per detector not in synch with TreeK
 168             ReorderAndExpandTreeTR();
 169         } else {
 170             // New format
 171             fTreeTR->SetBranchAddress("TrackReferences", &fTrackReferences);
 172         }
 173     }
 174
 175     //
 176     fNparticles = fStack->GetNtrack();
 177     fNprimaries = fStack->GetNprimary();
 178     AliInfo(Form("AliMCEventHandler: Event#: %5d Number of particles %5d \n", fEvent, fNparticles));
 179
 180     return kTRUE;
 181 }
 182
 183 Bool_t AliMCEventHandler::OpenFile(Int_t i)
 184 {
 185     // Open file i
 186     Bool_t ok = kTRUE;
 187     if (i > 0) {
 188         fExtension = Form("%d", i);
 189     } else {
 190         fExtension = "";
 191     }
 192
 193
 194     delete fFileK;
 195     fFileK = new TFile(Form("%sKinematics%s.root", fPathName, fExtension));
 196     if (!fFileK) {
 197         AliFatal(Form("AliMCEventHandler:Kinematics%s.root not found in directory %s ! \n", fExtension, fPathName));
 198         ok = kFALSE;
 199     }
 200
 201     delete fFileTR;
 202     fFileTR = new TFile(Form("%sTrackRefs%s.root", fPathName, fExtension));
 203     if (!fFileTR) {
 204         AliWarning(Form("AliMCEventHandler:TrackRefs%s.root not found in directory %s ! \n", fExtension, fPathName));
 205         ok = kFALSE;
 206     }
 207
 208     return ok;
 209 }
 210
 211 Bool_t AliMCEventHandler::BeginEvent()
 212 {
 213     // Read the next event
 214     fEvent++;
 215     if (fEvent >= fNEvent) {
 216         AliWarning(Form("AliMCEventHandler: Event number out of range %5d\n", fEvent));
 217         return kFALSE;
 218     }
 219     return GetEvent(fEvent);
 220 }
 221
 222 Int_t AliMCEventHandler::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
 223 {
 224     // Retrieve entry i
 225     if (i > -1 && i < fNparticles) {
 226         AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
 227         particle = 0;
 228         trefs    = 0;
 229         return (-1);
 230     }
 231     particle = fStack->Particle(i);
 232     if (fFileTR) {
 233         fTreeTR->GetEntry(fStack->TreeKEntry(i));
 234         trefs    = fTrackReferences;
 235         return trefs->GetEntries();
 236     } else {
 237         trefs = 0;
 238         return -1;
 239     }
 240 }
 241
 242 void AliMCEventHandler::DrawCheck(Int_t i, Bool_t search)
 243 {
 244     // Retrieve entry i and draw momentum vector and hits
 245     if (!fFileTR) {
 246         AliWarning("No Track Reference information available");
 247         return;
 248     }
 249
 250     if (i > -1 && i < fNparticles) {
 251         fTreeTR->GetEntry(fStack->TreeKEntry(i));
 252     } else {
 253         AliWarning("AliMCEventHandler::GetEntry: Index out of range");
 254     }
 255
 256     Int_t nh = fTrackReferences->GetEntries();
 257
 258
 259     if (search) {
 260         while(nh == 0 && i < fNparticles - 1) {
 261             i++;
 262             fTreeTR->GetEntry(fStack->TreeKEntry(i));
 263             nh =  fTrackReferences->GetEntries();
 264         }
 265         printf("Found Hits at %5d\n", i);
 266     }
 267     TParticle* particle = fStack->Particle(i);
 268
 269     TH2F*    h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
 270     Float_t x0 = particle->Vx();
 271     Float_t y0 = particle->Vy();
 272
 273     Float_t x1 = particle->Vx() + particle->Px() * 50.;
 274     Float_t y1 = particle->Vy() + particle->Py() * 50.;
 275
 276     TArrow*  a = new TArrow(x0, y0, x1, y1, 0.01);
 277     h->Draw();
 278     a->SetLineColor(2);
 279
 280     a->Draw();
 281
 282     for (Int_t ih = 0; ih < nh; ih++) {
 283         AliTrackReference* ref = (AliTrackReference*) fTrackReferences->At(ih);
 284         TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
 285         m->Draw();
 286         m->SetMarkerSize(0.4);
 287
 288     }
 289 }
 290
 291 Bool_t AliMCEventHandler::Notify(const char *path)
 292 {
 293     // Notify about directory change
 294     // The directory is taken from the 'path' argument
 295     // Reconnect trees
 296
 297     printf("AliMCEventHandler::Notify() file: %s\n", path);
 298     fPathName = Form("%s",  path);
 299     ResetIO();
 300     InitIO("");
 301     return kTRUE;
 302 }
 303
 304 void AliMCEventHandler::ResetIO()
 305 {
 306     // Reset files
 307     if (fFileE)  delete fFileE;
 308     if (fFileK)  delete fFileK;
 309     if (fFileTR) delete fFileTR;
 310 }
 311
 312
 313 Bool_t AliMCEventHandler::FinishEvent()
 314 {
 315     // Dummy
 316     return kTRUE;
 317 }
 318
 319 Bool_t AliMCEventHandler::Terminate()
 320 {
 321     // Dummy
 322     return kTRUE;
 323 }
 324
 325 Bool_t AliMCEventHandler::TerminateIO()
 326 {
 327     // Dummy
 328     return kTRUE;
 329 }
 330
 331 void AliMCEventHandler::ReorderAndExpandTreeTR()
 332 {
 333 //
 334 //  Reorder and expand the track reference tree in order to match the kinematics tree.
 335 //  Copy the information from different branches into one
 336 //
 337 //  TreeTR
 338     if (fTmpTreeTR) delete fTmpTreeTR;
 339     if (fTmpFileTR) {
 340         fTmpFileTR->Close();
 341         delete fTmpFileTR;
 342     }
 343
 344     fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
 345     fTmpTreeTR = new TTree("TreeTR", "Track References");
 346     if (!fTrackReferences)  fTrackReferences = new TClonesArray("AliTrackReference", 100);
 347     fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTrackReferences, 4000);
 348 //
 349     TClonesArray* trefs[7];
 350     for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
 351     if (fTreeTR){
 352         // make branch for central track references
 353         if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
 354         if (fTreeTR->GetBranch("ITS"))    fTreeTR->SetBranchAddress("ITS",    &trefs[1]);
 355         if (fTreeTR->GetBranch("TPC"))    fTreeTR->SetBranchAddress("TPC",    &trefs[2]);
 356         if (fTreeTR->GetBranch("TRD"))    fTreeTR->SetBranchAddress("TRD",    &trefs[3]);
 357         if (fTreeTR->GetBranch("TOF"))    fTreeTR->SetBranchAddress("TOF",    &trefs[4]);
 358         if (fTreeTR->GetBranch("FRAME"))  fTreeTR->SetBranchAddress("FRAME",  &trefs[5]);
 359         if (fTreeTR->GetBranch("MUON"))   fTreeTR->SetBranchAddress("MUON",   &trefs[6]);
 360     }
 361
 362     Int_t np = fStack->GetNprimary();
 363     Int_t nt = fTreeTR->GetEntries();
 364     //
 365     // Loop over tracks and find the secondaries with the help of the kine tree
 366     Int_t ifills = 0;
 367     Int_t it     = 0;
 368     Int_t itlast = 0;
 369
 370     for (Int_t ip = np - 1; ip > -1; ip--) {
 371         TParticle *part = fStack->Particle(ip);
 372 //      printf("Particle %5d %5d %5d %5d %5d %5d \n",
 373 //             ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
 374 //             part->GetLastDaughter(), part->TestBit(kTransportBit));
 375
 376         // Determine range of secondaries produced by this primary during transport
 377         Int_t dau1  = part->GetFirstDaughter();
 378         if (dau1 < np) continue;  // This particle has no secondaries produced during transport
 379         Int_t dau2  = -1;
 380         if (dau1 > -1) {
 381             Int_t inext = ip - 1;
 382             while (dau2 < 0) {
 383                 if (inext >= 0) {
 384                     part = fStack->Particle(inext);
 385                     dau2 =  part->GetFirstDaughter();
 386                     if (dau2 == -1 || dau2 < np) {
 387                         dau2 = -1;
 388                     } else {
 389                         dau2--;
 390                     }
 391                 } else {
 392                     dau2 = fStack->GetNtrack() - 1;
 393                 }
 394                 inext--;
 395             } // find upper bound
 396         }  // dau2 < 0
 397
 398
 399 //      printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
 400 //
 401 // Loop over reference hits and find secondary label
 402 // First the tricky part: find the entry in treeTR than contains the hits or
 403 // make sure that no hits exist.
 404 //
 405         Bool_t hasHits   = kFALSE;
 406         Bool_t isOutside = kFALSE;
 407
 408         it = itlast;
 409         while (!hasHits && !isOutside && it < nt) {
 410             fTreeTR->GetEntry(it++);
 411             for (Int_t ib = 0; ib < 7; ib++) {
 412                 if (!trefs[ib]) continue;
 413                 Int_t nh = trefs[ib]->GetEntries();
 414                 for (Int_t ih = 0; ih < nh; ih++) {
 415                     AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
 416                     Int_t label = tr->Label();
 417                     if (label >= dau1 && label <= dau2) {
 418                         hasHits = kTRUE;
 419                         itlast = it - 1;
 420                         break;
 421                     }
 422                     if (label > dau2 || label < ip) {
 423                         isOutside = kTRUE;
 424                         itlast = it - 1;
 425                         break;
 426                     }
 427                 } // hits
 428                 if (hasHits || isOutside) break;
 429             } // branches
 430         } // entries
 431
 432         if (!hasHits) {
 433             // Write empty entries
 434             for (Int_t id = dau1; (id <= dau2); id++) {
 435                 fTmpTreeTR->Fill();
 436                 ifills++;
 437             }
 438         } else {
 439             // Collect all hits
 440             fTreeTR->GetEntry(itlast);
 441             for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
 442                 for (Int_t ib = 0; ib < 7; ib++) {
 443                     if (!trefs[ib]) continue;
 444                     Int_t nh = trefs[ib]->GetEntries();
 445                     for (Int_t ih = 0; ih < nh; ih++) {
 446                         AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
 447                         Int_t label = tr->Label();
 448                         // Skip primaries
 449                         if (label == ip) continue;
 450                         if (label > dau2 || label < dau1)
 451                             printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
 452                                    itlast, label, dau1, dau2);
 453                         if (label == id) {
 454                             // secondary found
 455                             tr->SetDetectorId(ib-1);
 456                             Int_t nref =  fTrackReferences->GetEntriesFast();
 457                             TClonesArray &lref = *fTrackReferences;
 458                             new(lref[nref]) AliTrackReference(*tr);
 459                         }
 460                     } // hits
 461                 } // branches
 462                 fTmpTreeTR->Fill();
 463                 fTrackReferences->Clear();
 464                 ifills++;
 465             } // daughters
 466         } // has hits
 467     } // tracks
 468     //
 469     // Now loop again and write the primaries
 470     //
 471     it = nt - 1;
 472     for (Int_t ip = 0; ip < np; ip++) {
 473         Int_t labmax = -1;
 474         while (labmax < ip && it > -1) {
 475             fTreeTR->GetEntry(it--);
 476             for (Int_t ib = 0; ib < 7; ib++) {
 477                 if (!trefs[ib]) continue;
 478                 Int_t nh = trefs[ib]->GetEntries();
 479                 //
 480                 // Loop over reference hits and find primary labels
 481                 for (Int_t ih = 0; ih < nh; ih++) {
 482                     AliTrackReference* tr = (AliTrackReference*)  trefs[ib]->At(ih);
 483                     Int_t label = tr->Label();
 484                     if (label < np && label > labmax) {
 485                         labmax = label;
 486                     }
 487
 488                     if (label == ip) {
 489                         tr->SetDetectorId(ib-1);
 490                         Int_t nref = fTrackReferences->GetEntriesFast();
 491                         TClonesArray &lref = *fTrackReferences;
 492                         new(lref[nref]) AliTrackReference(*tr);
 493                     }
 494                 } // hits
 495             } // branches
 496         } // entries
 497         it++;
 498         fTmpTreeTR->Fill();
 499         fTrackReferences->Clear();
 500         ifills++;
 501     } // tracks
 502     // Check
 503     if (ifills != fStack->GetNtrack())
 504         printf("AliMCEventHandler:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
 505                ifills, fStack->GetNtrack());
 506     fTreeTR = fTmpTreeTR;
 507 }