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 //
 362 //    fTreeTR->SetBranchStatus("*",      0);
 363     fTreeTR->SetBranchStatus("ABSO",  0);
 364     fTreeTR->SetBranchStatus("BODY",  0);
 365     fTreeTR->SetBranchStatus("DIPO",  0);
 366     fTreeTR->SetBranchStatus("EMCAL", 0);
 367     fTreeTR->SetBranchStatus("FMD",   0);
 368     fTreeTR->SetBranchStatus("HALL",  0);
 369     fTreeTR->SetBranchStatus("MAG",   0);
 370     fTreeTR->SetBranchStatus("PHOS",  0);
 371     fTreeTR->SetBranchStatus("PIPE",  0);
 372     fTreeTR->SetBranchStatus("PMD",   0);
 373     fTreeTR->SetBranchStatus("RICH",  0);
 374     fTreeTR->SetBranchStatus("SHIL",  0);
 375     fTreeTR->SetBranchStatus("START", 0);
 376     fTreeTR->SetBranchStatus("VZERO", 0);
 377     fTreeTR->SetBranchStatus("ZDC",   0);
 378
 379     fTreeTR->SetBranchStatus("AliRun", 1);
 380     fTreeTR->SetBranchStatus("ITS",    1);
 381     fTreeTR->SetBranchStatus("TPC",    1);
 382     fTreeTR->SetBranchStatus("TRD",    1);
 383     fTreeTR->SetBranchStatus("TOF",    1);
 384     fTreeTR->SetBranchStatus("FRAME",  1);
 385     fTreeTR->SetBranchStatus("MUON",   1);
 386
 387     TClonesArray* trefs[7];
 388     for (Int_t i = 0; i < 7; i++) trefs[i] = 0;
 389     if (fTreeTR){
 390         // make branch for central track references
 391         if (fTreeTR->GetBranch("AliRun")) fTreeTR->SetBranchAddress("AliRun", &trefs[0]);
 392         if (fTreeTR->GetBranch("ITS"))    fTreeTR->SetBranchAddress("ITS",    &trefs[1]);
 393         if (fTreeTR->GetBranch("TPC"))    fTreeTR->SetBranchAddress("TPC",    &trefs[2]);
 394         if (fTreeTR->GetBranch("TRD"))    fTreeTR->SetBranchAddress("TRD",    &trefs[3]);
 395         if (fTreeTR->GetBranch("TOF"))    fTreeTR->SetBranchAddress("TOF",    &trefs[4]);
 396         if (fTreeTR->GetBranch("FRAME"))  fTreeTR->SetBranchAddress("FRAME",  &trefs[5]);
 397         if (fTreeTR->GetBranch("MUON"))   fTreeTR->SetBranchAddress("MUON",   &trefs[6]);
 398     }
 399
 400     Int_t np = fStack->GetNprimary();
 401     Int_t nt = fTreeTR->GetEntries();
 402     //
 403     // Loop over tracks and find the secondaries with the help of the kine tree
 404     Int_t ifills = 0;
 405     Int_t it     = 0;
 406     Int_t itlast = 0;
 407
 408     for (Int_t ip = np - 1; ip > -1; ip--) {
 409         TParticle *part = fStack->Particle(ip);
 410 //      printf("Particle %5d %5d %5d %5d %5d %5d \n",
 411 //             ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
 412 //             part->GetLastDaughter(), part->TestBit(kTransportBit));
 413
 414         // Determine range of secondaries produced by this primary during transport
 415         Int_t dau1  = part->GetFirstDaughter();
 416         if (dau1 < np) continue;  // This particle has no secondaries produced during transport
 417         Int_t dau2  = -1;
 418         if (dau1 > -1) {
 419             Int_t inext = ip - 1;
 420             while (dau2 < 0) {
 421                 if (inext >= 0) {
 422                     part = fStack->Particle(inext);
 423                     dau2 =  part->GetFirstDaughter();
 424                     if (dau2 == -1 || dau2 < np) {
 425                         dau2 = -1;
 426                     } else {
 427                         dau2--;
 428                     }
 429                 } else {
 430                     dau2 = fStack->GetNtrack() - 1;
 431                 }
 432                 inext--;
 433             } // find upper bound
 434         }  // dau2 < 0
 435
 436
 437 //      printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
 438 //
 439 // Loop over reference hits and find secondary label
 440 // First the tricky part: find the entry in treeTR than contains the hits or
 441 // make sure that no hits exist.
 442 //
 443         Bool_t hasHits   = kFALSE;
 444         Bool_t isOutside = kFALSE;
 445
 446         it = itlast;
 447         while (!hasHits && !isOutside && it < nt) {
 448             fTreeTR->GetEntry(it++);
 449             for (Int_t ib = 0; ib < 7; ib++) {
 450                 if (!trefs[ib]) continue;
 451                 Int_t nh = trefs[ib]->GetEntries();
 452                 for (Int_t ih = 0; ih < nh; ih++) {
 453                     AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
 454                     Int_t label = tr->Label();
 455                     if (label >= dau1 && label <= dau2) {
 456                         hasHits = kTRUE;
 457                         itlast = it - 1;
 458                         break;
 459                     }
 460                     if (label > dau2 || label < ip) {
 461                         isOutside = kTRUE;
 462                         itlast = it - 1;
 463                         break;
 464                     }
 465                 } // hits
 466                 if (hasHits || isOutside) break;
 467             } // branches
 468         } // entries
 469
 470         if (!hasHits) {
 471             // Write empty entries
 472             for (Int_t id = dau1; (id <= dau2); id++) {
 473                 fTmpTreeTR->Fill();
 474                 ifills++;
 475             }
 476         } else {
 477             // Collect all hits
 478             fTreeTR->GetEntry(itlast);
 479             for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
 480                 for (Int_t ib = 0; ib < 7; ib++) {
 481                     if (!trefs[ib]) continue;
 482                     Int_t nh = trefs[ib]->GetEntries();
 483                     for (Int_t ih = 0; ih < nh; ih++) {
 484                         AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
 485                         Int_t label = tr->Label();
 486                         // Skip primaries
 487                         if (label == ip) continue;
 488                         if (label > dau2 || label < dau1)
 489                             printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
 490                                    itlast, label, dau1, dau2);
 491                         if (label == id) {
 492                             // secondary found
 493                             tr->SetDetectorId(ib-1);
 494                             Int_t nref =  fTrackReferences->GetEntriesFast();
 495                             TClonesArray &lref = *fTrackReferences;
 496                             new(lref[nref]) AliTrackReference(*tr);
 497                         }
 498                     } // hits
 499                 } // branches
 500                 fTmpTreeTR->Fill();
 501                 fTrackReferences->Clear();
 502                 ifills++;
 503             } // daughters
 504         } // has hits
 505     } // tracks
 506     //
 507     // Now loop again and write the primaries
 508     //
 509     it = nt - 1;
 510     for (Int_t ip = 0; ip < np; ip++) {
 511         Int_t labmax = -1;
 512         while (labmax < ip && it > -1) {
 513             fTreeTR->GetEntry(it--);
 514             for (Int_t ib = 0; ib < 7; ib++) {
 515                 if (!trefs[ib]) continue;
 516                 Int_t nh = trefs[ib]->GetEntries();
 517                 //
 518                 // Loop over reference hits and find primary labels
 519                 for (Int_t ih = 0; ih < nh; ih++) {
 520                     AliTrackReference* tr = (AliTrackReference*)  trefs[ib]->At(ih);
 521                     Int_t label = tr->Label();
 522                     if (label < np && label > labmax) {
 523                         labmax = label;
 524                     }
 525
 526                     if (label == ip) {
 527                         tr->SetDetectorId(ib-1);
 528                         Int_t nref = fTrackReferences->GetEntriesFast();
 529                         TClonesArray &lref = *fTrackReferences;
 530                         new(lref[nref]) AliTrackReference(*tr);
 531                     }
 532                 } // hits
 533             } // branches
 534         } // entries
 535         it++;
 536         fTmpTreeTR->Fill();
 537         fTrackReferences->Clear();
 538         ifills++;
 539     } // tracks
 540     // Check
 541     delete fTreeTR;
 542     for (Int_t ib = 0; ib < 7; ib++) {
 543         if (trefs[ib]) delete trefs[ib];
 544     }
 545
 546     if (ifills != fStack->GetNtrack())
 547         printf("AliMCEventHandler:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
 548                ifills, fStack->GetNtrack());
 549     fTreeTR = fTmpTreeTR;
 550 }
 551
 552 void AliMCEventHandler::SetInputPath(char* fname)
 553 {
 554     // Set the input path name
 555     delete fPathName;
 556     fPathName = new TString(fname);
 557 }