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