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