STEER/AliMCEvent.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2007, 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 //-------------------------------------------------------------------------
  19 //     Class for Kinematic Events
  20 //     Author: Andreas Morsch, CERN
  21 //-------------------------------------------------------------------------
  22 #include <TArrow.h>
  23 #include <TMarker.h>
  24 #include <TH2F.h>
  25 #include <TTree.h>
  26 #include <TFile.h>
  27 #include <TParticle.h>
  28 #include <TClonesArray.h>
  29 #include <TRefArray.h>
  30
  31 #include "AliLog.h"
  32 #include "AliMCEvent.h"
  33 #include "AliStack.h"
  34 #include "AliTrackReference.h"
  35 #include "AliHeader.h"
  36 #include "AliGenEventHeader.h"
  37
  38
  39 AliMCEvent::AliMCEvent():
  40     AliVEvent(),
  41     fStack(0),
  42     fMCParticles(new TClonesArray("AliMCParticle",1000)),
  43     fMCParticleMap(0),
  44     fHeader(new AliHeader()),
  45     fTRBuffer(0),
  46     fTrackReferences(new TClonesArray("AliTrackReference", 1000)),
  47     fTreeTR(0),
  48     fTmpTreeTR(0),
  49     fTmpFileTR(0),
  50     fNprimaries(-1),
  51     fNparticles(-1)
  52 {
  53     // Default constructor
  54 }
  55
  56 AliMCEvent::AliMCEvent(const AliMCEvent& mcEvnt) :
  57     AliVEvent(mcEvnt),
  58     fStack(0),
  59     fMCParticles(0),
  60     fMCParticleMap(0),
  61     fHeader(new AliHeader()),
  62     fTRBuffer(0),
  63     fTrackReferences(0),
  64     fTreeTR(0),
  65     fTmpTreeTR(0),
  66     fTmpFileTR(0),
  67     fNprimaries(-1),
  68     fNparticles(-1)
  69 {
  70 // Copy constructor
  71 }
  72
  73
  74 AliMCEvent& AliMCEvent::operator=(const AliMCEvent& mcEvnt)
  75 {
  76     // assignment operator
  77     if (this!=&mcEvnt) {
  78         AliVEvent::operator=(mcEvnt);
  79     }
  80
  81     return *this;
  82 }
  83
  84 void AliMCEvent::ConnectTreeE (TTree* tree)
  85 {
  86     // Connect the event header tree
  87     tree->SetBranchAddress("Header", &fHeader);
  88 }
  89
  90 void AliMCEvent::ConnectTreeK (TTree* tree)
  91 {
  92     // Connect the kinematics tree to the stack
  93     fStack = fHeader->Stack();
  94     fStack->ConnectTree(tree);
  95     //
  96     // Load the event
  97     fStack->GetEvent();
  98     fNparticles = fStack->GetNtrack();
  99     fNprimaries = fStack->GetNprimary();
 100     Int_t iev  = fHeader->GetEvent();
 101     Int_t ievr = fHeader->GetEventNrInRun();
 102
 103     AliInfo(Form("AliMCEvent# %5d %5d: Number of particles: %5d (all) %5d (primaries)\n",
 104                  iev, ievr, fNparticles, fNprimaries));
 105
 106     // This is a cache for the TParticles converted to MCParticles on user request
 107     if (fMCParticleMap) {
 108         fMCParticleMap->Clear();
 109         fMCParticles->Delete();
 110         if (fNparticles>0) fMCParticleMap->Expand(fNparticles);
 111     }
 112     else
 113         fMCParticleMap = new TRefArray(fNparticles);
 114 }
 115
 116 void AliMCEvent::ConnectTreeTR (TTree* tree)
 117 {
 118     // Connect the track reference tree
 119     fTreeTR = tree;
 120
 121     if (fTreeTR->GetBranch("AliRun")) {
 122         if (fTmpFileTR) {
 123             fTmpFileTR->Close();
 124             delete fTmpFileTR;
 125         }
 126         // This is an old format with one branch per detector not in synch with TreeK
 127         ReorderAndExpandTreeTR();
 128     } else {
 129         // New format
 130         fTreeTR->SetBranchAddress("TrackReferences", &fTRBuffer);
 131     }
 132 }
 133
 134 Int_t AliMCEvent::GetParticleAndTR(Int_t i, TParticle*& particle, TClonesArray*& trefs)
 135 {
 136     // Retrieve entry i
 137     if (i < 0 || i >= fNparticles) {
 138         AliWarning(Form("AliMCEventHandler::GetEntry: Index out of range"));
 139         particle = 0;
 140         trefs    = 0;
 141         return (-1);
 142     }
 143     particle = fStack->Particle(i);
 144     if (fTreeTR) {
 145         fTreeTR->GetEntry(fStack->TreeKEntry(i));
 146         trefs    = fTRBuffer;
 147         return trefs->GetEntries();
 148     } else {
 149         trefs = 0;
 150         return -1;
 151     }
 152 }
 153
 154
 155 void AliMCEvent::Clean()
 156 {
 157     // Clean-up before new trees are connected
 158 //    if (fHeader) {
 159 //      delete fHeader;
 160 //      fHeader = 0;
 161 //    }
 162
 163     delete fStack; fStack = 0;
 164
 165     // Clear TR
 166     if (fTRBuffer) {
 167         fTRBuffer->Delete();
 168         delete fTRBuffer;
 169         fTRBuffer = 0;
 170     }
 171 }
 172
 173 void AliMCEvent::FinishEvent()
 174 {
 175   // Clean-up after event
 176   //
 177     fStack->Reset(0);
 178     fMCParticles->Delete();
 179     fMCParticleMap->Clear();
 180     if (fTRBuffer)
 181       fTRBuffer->Delete();
 182     fTrackReferences->Delete();
 183     fNparticles = -1;
 184     fNprimaries = -1;
 185     fStack      =  0;
 186
 187 }
 188
 189
 190
 191 void AliMCEvent::DrawCheck(Int_t i, Int_t search)
 192 {
 193     //
 194     // Simple event display for debugging
 195     if (!fTreeTR) {
 196         AliWarning("No Track Reference information available");
 197         return;
 198     }
 199
 200     if (i > -1 && i < fNparticles) {
 201         fTreeTR->GetEntry(fStack->TreeKEntry(i));
 202     } else {
 203         AliWarning("AliMCEvent::GetEntry: Index out of range");
 204     }
 205
 206     Int_t nh = fTRBuffer->GetEntries();
 207
 208
 209     if (search) {
 210         while(nh <= search && i < fNparticles - 1) {
 211             i++;
 212             fTreeTR->GetEntry(fStack->TreeKEntry(i));
 213             nh =  fTRBuffer->GetEntries();
 214         }
 215         printf("Found Hits at %5d\n", i);
 216     }
 217     TParticle* particle = fStack->Particle(i);
 218
 219     TH2F*    h = new TH2F("", "", 100, -500, 500, 100, -500, 500);
 220     Float_t x0 = particle->Vx();
 221     Float_t y0 = particle->Vy();
 222
 223     Float_t x1 = particle->Vx() + particle->Px() * 50.;
 224     Float_t y1 = particle->Vy() + particle->Py() * 50.;
 225
 226     TArrow*  a = new TArrow(x0, y0, x1, y1, 0.01);
 227     h->Draw();
 228     a->SetLineColor(2);
 229
 230     a->Draw();
 231
 232     for (Int_t ih = 0; ih < nh; ih++) {
 233         AliTrackReference* ref = (AliTrackReference*) fTRBuffer->At(ih);
 234         TMarker* m = new TMarker(ref->X(), ref->Y(), 20);
 235         m->Draw();
 236         m->SetMarkerSize(0.4);
 237
 238     }
 239 }
 240
 241
 242 void AliMCEvent::ReorderAndExpandTreeTR()
 243 {
 244 //
 245 //  Reorder and expand the track reference tree in order to match the kinematics tree.
 246 //  Copy the information from different branches into one
 247 //
 248 //  TreeTR
 249
 250     fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate");
 251     fTmpTreeTR = new TTree("TreeTR", "TrackReferences");
 252     if (!fTRBuffer)  fTRBuffer = new TClonesArray("AliTrackReference", 100);
 253     fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &fTRBuffer, 32000, 0);
 254
 255
 256 //
 257 //  Activate the used branches only. Otherwisw we get a bad memory leak.
 258     fTreeTR->SetBranchStatus("*",        0);
 259     fTreeTR->SetBranchStatus("AliRun.*", 1);
 260     fTreeTR->SetBranchStatus("ITS.*",    1);
 261     fTreeTR->SetBranchStatus("TPC.*",    1);
 262     fTreeTR->SetBranchStatus("TRD.*",    1);
 263     fTreeTR->SetBranchStatus("TOF.*",    1);
 264     fTreeTR->SetBranchStatus("FRAME.*",  1);
 265     fTreeTR->SetBranchStatus("MUON.*",   1);
 266 //
 267 //  Connect the active branches
 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     //
 285     // Loop over tracks and find the secondaries with the help of the kine tree
 286     Int_t ifills = 0;
 287     Int_t it     = 0;
 288     Int_t itlast = 0;
 289     TParticle* part;
 290
 291     for (Int_t ip = np - 1; ip > -1; ip--) {
 292         part = fStack->Particle(ip);
 293 //      printf("Particle %5d %5d %5d %5d %5d %5d \n",
 294 //             ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(),
 295 //             part->GetLastDaughter(), part->TestBit(kTransportBit));
 296
 297         // Determine range of secondaries produced by this primary during transport
 298         Int_t dau1  = part->GetFirstDaughter();
 299         if (dau1 < np) continue;  // This particle has no secondaries produced during transport
 300         Int_t dau2  = -1;
 301         if (dau1 > -1) {
 302             Int_t inext = ip - 1;
 303             while (dau2 < 0) {
 304                 if (inext >= 0) {
 305                     part = fStack->Particle(inext);
 306                     dau2 =  part->GetFirstDaughter();
 307                     if (dau2 == -1 || dau2 < np) {
 308                         dau2 = -1;
 309                     } else {
 310                         dau2--;
 311                     }
 312                 } else {
 313                     dau2 = fStack->GetNtrack() - 1;
 314                 }
 315                 inext--;
 316             } // find upper bound
 317         }  // dau2 < 0
 318
 319
 320 //      printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2);
 321 //
 322 // Loop over reference hits and find secondary label
 323 // First the tricky part: find the entry in treeTR than contains the hits or
 324 // make sure that no hits exist.
 325 //
 326         Bool_t hasHits   = kFALSE;
 327         Bool_t isOutside = kFALSE;
 328
 329         it = itlast;
 330         while (!hasHits && !isOutside && it < nt) {
 331             fTreeTR->GetEntry(it++);
 332             for (Int_t ib = 0; ib < 7; ib++) {
 333                 if (!trefs[ib]) continue;
 334                 Int_t nh = trefs[ib]->GetEntries();
 335                 for (Int_t ih = 0; ih < nh; ih++) {
 336                     AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
 337                     Int_t label = tr->Label();
 338                     if (label >= dau1 && label <= dau2) {
 339                         hasHits = kTRUE;
 340                         itlast = it - 1;
 341                         break;
 342                     }
 343                     if (label > dau2 || label < ip) {
 344                         isOutside = kTRUE;
 345                         itlast = it - 1;
 346                         break;
 347                     }
 348                 } // hits
 349                 if (hasHits || isOutside) break;
 350             } // branches
 351         } // entries
 352
 353         if (!hasHits) {
 354             // Write empty entries
 355             for (Int_t id = dau1; (id <= dau2); id++) {
 356                 fTmpTreeTR->Fill();
 357                 ifills++;
 358             }
 359         } else {
 360             // Collect all hits
 361             fTreeTR->GetEntry(itlast);
 362             for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) {
 363                 for (Int_t ib = 0; ib < 7; ib++) {
 364                     if (!trefs[ib]) continue;
 365                     Int_t nh = trefs[ib]->GetEntries();
 366                     for (Int_t ih = 0; ih < nh; ih++) {
 367                         AliTrackReference* tr = (AliTrackReference*) trefs[ib]->At(ih);
 368                         Int_t label = tr->Label();
 369                         // Skip primaries
 370                         if (label == ip) continue;
 371                         if (label > dau2 || label < dau1)
 372                             printf("AliMCEventHandler::Track Reference Label out of range !: %5d %5d %5d %5d \n",
 373                                    itlast, label, dau1, dau2);
 374                         if (label == id) {
 375                             // secondary found
 376                             tr->SetDetectorId(ib-1);
 377                             Int_t nref =  fTRBuffer->GetEntriesFast();
 378                             TClonesArray &lref = *fTRBuffer;
 379                             new(lref[nref]) AliTrackReference(*tr);
 380                         }
 381                     } // hits
 382                 } // branches
 383                 fTmpTreeTR->Fill();
 384                 fTRBuffer->Delete();
 385                 ifills++;
 386             } // daughters
 387         } // has hits
 388     } // tracks
 389
 390     //
 391     // Now loop again and write the primaries
 392     //
 393     it = nt - 1;
 394     for (Int_t ip = 0; ip < np; ip++) {
 395         Int_t labmax = -1;
 396         while (labmax < ip && it > -1) {
 397             fTreeTR->GetEntry(it--);
 398             for (Int_t ib = 0; ib < 7; ib++) {
 399                 if (!trefs[ib]) continue;
 400                 Int_t nh = trefs[ib]->GetEntries();
 401                 //
 402                 // Loop over reference hits and find primary labels
 403                 for (Int_t ih = 0; ih < nh; ih++) {
 404                     AliTrackReference* tr = (AliTrackReference*)  trefs[ib]->At(ih);
 405                     Int_t label = tr->Label();
 406                     if (label < np && label > labmax) {
 407                         labmax = label;
 408                     }
 409
 410                     if (label == ip) {
 411                         tr->SetDetectorId(ib-1);
 412                         Int_t nref = fTRBuffer->GetEntriesFast();
 413                         TClonesArray &lref = *fTRBuffer;
 414                         new(lref[nref]) AliTrackReference(*tr);
 415                     }
 416                 } // hits
 417             } // branches
 418         } // entries
 419         it++;
 420         fTmpTreeTR->Fill();
 421         fTRBuffer->Delete();
 422         ifills++;
 423     } // tracks
 424     // Check
 425
 426
 427     // Clean-up
 428     delete fTreeTR; fTreeTR = 0;
 429
 430     for (Int_t ib = 0; ib < 7; ib++) {
 431         if (trefs[ib]) {
 432             trefs[ib]->Clear();
 433             delete trefs[ib];
 434             trefs[ib] = 0;
 435         }
 436     }
 437
 438     if (ifills != fStack->GetNtrack())
 439         printf("AliMCEvent:Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n",
 440                ifills, fStack->GetNtrack());
 441
 442     fTmpTreeTR->Write();
 443     fTreeTR = fTmpTreeTR;
 444 }
 445
 446 AliMCParticle* AliMCEvent::GetTrack(Int_t i) const
 447 {
 448     // Get MC Particle i
 449     AliMCParticle *mcParticle = 0;
 450     TParticle     *particle   = 0;
 451     TClonesArray  *trefs      = 0;
 452     Int_t          ntref      = 0;
 453     TRefArray     *rarray     = 0;
 454     // Out of range check
 455     if (i < 0 || i >= fNparticles) {
 456         AliWarning(Form("AliMCEvent::GetEntry: Index out of range"));
 457         mcParticle = 0;
 458         return (mcParticle);
 459     }
 460
 461
 462     //
 463     // First check of the MC Particle has been already cached
 464     if(!fMCParticleMap->At(i)) {
 465         // Get particle from the stack
 466         particle   = fStack->Particle(i);
 467         // Get track references from Tree TR
 468         if (fTreeTR) {
 469             fTreeTR->GetEntry(fStack->TreeKEntry(i));
 470             trefs     = fTRBuffer;
 471             ntref     = trefs->GetEntriesFast();
 472             rarray    = new TRefArray(ntref);
 473             Int_t nen = fTrackReferences->GetEntriesFast();
 474             for (Int_t j = 0; j < ntref; j++) {
 475                 // Save the track references in a TClonesArray
 476                 AliTrackReference* ref = dynamic_cast<AliTrackReference*>((*fTRBuffer)[j]);
 477                 // Save the pointer in a TRefArray
 478                 new ((*fTrackReferences)[nen]) AliTrackReference(*ref);
 479                 rarray->AddAt((*fTrackReferences)[nen], j);
 480                 nen++;
 481             } // loop over track references for entry i
 482         } // if TreeTR available
 483         Int_t nentries = fMCParticles->GetEntriesFast();
 484         new ((*fMCParticles)[nentries]) AliMCParticle(particle, rarray, i);
 485         mcParticle = dynamic_cast<AliMCParticle*>((*fMCParticles)[nentries]);
 486         fMCParticleMap->AddAt(mcParticle, i);
 487     } else {
 488         mcParticle = dynamic_cast<AliMCParticle*>(fMCParticleMap->At(i));
 489     }
 490     return mcParticle;
 491 }
 492
 493  AliGenEventHeader* AliMCEvent::GenEventHeader() {return (fHeader->GenEventHeader());}
 494
 495 ClassImp(AliMCEvent)