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