Major upgrade of AliRoot code
[u/mrichter/AliRoot.git] / STEER / AliRun.cxx
CommitLineData
99d554c8 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/*
17$Log$
e460afec 18Revision 1.47 2000/12/18 10:44:01 morsch
19Possibility to set field map by passing pointer to objet of type AliMagF via
20SetField().
21Example:
22gAlice->SetField(new AliMagFCM("Map2", "$(ALICE_ROOT)/data/field01.dat",2,1.,10.));
23
d8408e76 24Revision 1.46 2000/12/14 19:29:27 fca
25galice.cuts was not read any more
26
327136d2 27Revision 1.45 2000/11/30 07:12:49 alibrary
28Introducing new Rndm and QA classes
29
65fb704d 30Revision 1.44 2000/10/26 13:58:59 morsch
31Add possibility to choose the lego generator (of type AliGeneratorLego or derived) when running
32RunLego(). Default is the base class AliGeneratorLego.
33
0a520a66 34Revision 1.43 2000/10/09 09:43:17 fca
35Special remapping of hits for TPC and TRD. End-of-primary action introduced
36
24de2263 37Revision 1.42 2000/10/02 21:28:14 fca
38Removal of useless dependecies via forward declarations
39
94de3818 40Revision 1.41 2000/07/13 16:19:09 fca
41Mainly coding conventions + some small bug fixes
42
ef42d733 43Revision 1.40 2000/07/12 08:56:25 fca
44Coding convention correction and warning removal
45
8918e700 46Revision 1.39 2000/07/11 18:24:59 fca
47Coding convention corrections + few minor bug fixes
48
aee8290b 49Revision 1.38 2000/06/20 13:05:45 fca
50Writing down the TREE headers before job starts
51
51e0e89d 52Revision 1.37 2000/06/09 20:05:11 morsch
53Introduce possibility to chose magnetic field version 3: AliMagFDM + field02.dat
54
f1b9d7c3 55Revision 1.36 2000/06/08 14:03:58 hristov
56Only one initializer for a default argument
57
d33c0226 58Revision 1.35 2000/06/07 10:13:14 hristov
59Delete only existent objects.
60
d2ecd553 61Revision 1.34 2000/05/18 10:45:38 fca
62Delete Particle Factory properly
63
c222d2b0 64Revision 1.33 2000/05/16 13:10:40 fca
65New method IsNewTrack and fix for a problem in Father-Daughter relations
66
a01a8b12 67Revision 1.32 2000/04/27 10:38:21 fca
68Correct termination of Lego Run and introduce Lego getter in AliRun
69
838edcaf 70Revision 1.31 2000/04/26 10:17:32 fca
71Changes in Lego for G4 compatibility
72
dffd31ef 73Revision 1.30 2000/04/18 19:11:40 fca
74Introduce variable Config.C function signature
75
45189757 76Revision 1.29 2000/04/07 11:12:34 fca
77G4 compatibility changes
78
875c717b 79Revision 1.28 2000/04/05 06:51:06 fca
80Workaround for an HP compiler problem
81
5eb58812 82Revision 1.27 2000/03/22 18:08:07 fca
83Rationalisation of the virtual MC interfaces
84
80762cb1 85Revision 1.26 2000/03/22 13:42:26 fca
86SetGenerator does not replace an existing generator, ResetGenerator does
87
ee1dd322 88Revision 1.25 2000/02/23 16:25:22 fca
89AliVMC and AliGeant3 classes introduced
90ReadEuclid moved from AliRun to AliModule
91
b13db077 92Revision 1.24 2000/01/19 17:17:20 fca
93Introducing a list of lists of hits -- more hits allowed for detector now
94
1cedd08a 95Revision 1.23 1999/12/03 11:14:31 fca
96Fixing previous wrong checking
97
00719c1b 98Revision 1.21 1999/11/25 10:40:08 fca
99Fixing daughters information also in primary tracks
100
ae23d366 101Revision 1.20 1999/10/04 18:08:49 fca
102Adding protection against inconsistent Euclid files
103
3fcc96a1 104Revision 1.19 1999/09/29 07:50:40 fca
105Introduction of the Copyright and cvs Log
106
99d554c8 107*/
108
fe4da5cc 109///////////////////////////////////////////////////////////////////////////////
110// //
111// Control class for Alice C++ //
112// Only one single instance of this class exists. //
113// The object is created in main program aliroot //
114// and is pointed by the global gAlice. //
115// //
8494b010 116// -Supports the list of all Alice Detectors (fModules). //
fe4da5cc 117// -Supports the list of particles (fParticles). //
118// -Supports the Trees. //
119// -Supports the geometry. //
120// -Supports the event display. //
121//Begin_Html
122/*
1439f98e 123<img src="picts/AliRunClass.gif">
fe4da5cc 124*/
125//End_Html
126//Begin_Html
127/*
1439f98e 128<img src="picts/alirun.gif">
fe4da5cc 129*/
130//End_Html
131// //
132///////////////////////////////////////////////////////////////////////////////
133
94de3818 134#include <stdlib.h>
135#include <stdio.h>
136#include <string.h>
137
fe4da5cc 138#include <TFile.h>
139#include <TRandom.h>
140#include <TBRIK.h>
141#include <TNode.h>
fe4da5cc 142#include <TCint.h>
143#include <TSystem.h>
a8f1fb7c 144#include <TObjectTable.h>
94de3818 145#include <TTree.h>
146#include <TGeometry.h>
147#include <TROOT.h>
148#include "TBrowser.h"
fe4da5cc 149
1578254f 150#include "TParticle.h"
fe4da5cc 151#include "AliRun.h"
fe4da5cc 152#include "AliDisplay.h"
875c717b 153#include "AliMC.h"
dffd31ef 154#include "AliLego.h"
aee8290b 155#include "AliMagFC.h"
156#include "AliMagFCM.h"
157#include "AliMagFDM.h"
94de3818 158#include "AliHit.h"
65fb704d 159#include "TRandom3.h"
160#include "AliMCQA.h"
161#include "AliGenerator.h"
162#include "AliLegoGenerator.h"
94de3818 163
164#include "AliDetector.h"
fe4da5cc 165
fe4da5cc 166AliRun *gAlice;
167
aee8290b 168static AliHeader *gAliHeader;
fe4da5cc 169
fe4da5cc 170ClassImp(AliRun)
171
172//_____________________________________________________________________________
173AliRun::AliRun()
174{
175 //
176 // Default constructor for AliRun
177 //
aee8290b 178 gAliHeader=&fHeader;
fe4da5cc 179 fRun = 0;
180 fEvent = 0;
181 fCurrent = -1;
8494b010 182 fModules = 0;
fe4da5cc 183 fGenerator = 0;
184 fTreeD = 0;
185 fTreeK = 0;
186 fTreeH = 0;
187 fTreeE = 0;
188 fTreeR = 0;
189 fParticles = 0;
190 fGeometry = 0;
191 fDisplay = 0;
192 fField = 0;
875c717b 193 fMC = 0;
fe4da5cc 194 fNdets = 0;
195 fImedia = 0;
196 fTrRmax = 1.e10;
197 fTrZmax = 1.e10;
fe4da5cc 198 fInitDone = kFALSE;
199 fLego = 0;
1578254f 200 fPDGDB = 0; //Particle factory object!
1cedd08a 201 fHitLists = 0;
45189757 202 fConfigFunction = "\0";
65fb704d 203 fRandom = 0;
204 fMCQA = 0;
205 fTransParName = "\0";
fe4da5cc 206}
207
208//_____________________________________________________________________________
209AliRun::AliRun(const char *name, const char *title)
210 : TNamed(name,title)
211{
212 //
213 // Constructor for the main processor.
214 // Creates the geometry
215 // Creates the list of Detectors.
216 // Creates the list of particles.
217 //
218 Int_t i;
219
220 gAlice = this;
221 fTreeD = 0;
222 fTreeK = 0;
223 fTreeH = 0;
224 fTreeE = 0;
225 fTreeR = 0;
226 fTrRmax = 1.e10;
227 fTrZmax = 1.e10;
1141f8e4 228 fGenerator = 0;
fe4da5cc 229 fInitDone = kFALSE;
230 fLego = 0;
231 fField = 0;
45189757 232 fConfigFunction = "Config();";
65fb704d 233
234 // Set random number generator
235 gRandom = fRandom = new TRandom3();
fe4da5cc 236
237 gROOT->GetListOfBrowsables()->Add(this,name);
238 //
239 // create the support list for the various Detectors
8494b010 240 fModules = new TObjArray(77);
fe4da5cc 241 //
242 // Create the TNode geometry for the event display
243
244 BuildSimpleGeometry();
245
246
247 fNtrack=0;
248 fHgwmk=0;
249 fCurrent=-1;
aee8290b 250 gAliHeader=&fHeader;
fe4da5cc 251 fRun = 0;
252 fEvent = 0;
253 //
254 // Create the particle stack
1578254f 255 fParticles = new TClonesArray("TParticle",100);
fe4da5cc 256
257 fDisplay = 0;
258 //
259 // Create default mag field
260 SetField();
261 //
875c717b 262 fMC = gMC;
fe4da5cc 263 //
264 // Prepare the tracking medium lists
265 fImedia = new TArrayI(1000);
266 for(i=0;i<1000;i++) (*fImedia)[i]=-99;
1578254f 267 //
268 // Make particles
269 fPDGDB = TDatabasePDG::Instance(); //Particle factory object!
1cedd08a 270 //
271 // Create HitLists list
272 fHitLists = new TList();
65fb704d 273 //
274 SetTransPar();
fe4da5cc 275}
276
aee8290b 277
fe4da5cc 278//_____________________________________________________________________________
279AliRun::~AliRun()
280{
281 //
282 // Defaullt AliRun destructor
283 //
fe4da5cc 284 delete fImedia;
285 delete fField;
875c717b 286 delete fMC;
fe4da5cc 287 delete fGeometry;
288 delete fDisplay;
289 delete fGenerator;
290 delete fLego;
291 delete fTreeD;
292 delete fTreeK;
293 delete fTreeH;
294 delete fTreeE;
295 delete fTreeR;
8494b010 296 if (fModules) {
297 fModules->Delete();
298 delete fModules;
fe4da5cc 299 }
300 if (fParticles) {
301 fParticles->Delete();
302 delete fParticles;
303 }
1cedd08a 304 delete fHitLists;
c222d2b0 305 delete fPDGDB;
e460afec 306 delete fMCQA;
fe4da5cc 307}
308
309//_____________________________________________________________________________
310void AliRun::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const
311{
312 //
313 // Add a hit to detector id
314 //
8494b010 315 TObjArray &dets = *fModules;
316 if(dets[id]) ((AliModule*) dets[id])->AddHit(track,vol,hits);
fe4da5cc 317}
318
319//_____________________________________________________________________________
320void AliRun::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const
321{
322 //
323 // Add digit to detector id
324 //
8494b010 325 TObjArray &dets = *fModules;
326 if(dets[id]) ((AliModule*) dets[id])->AddDigit(tracks,digits);
fe4da5cc 327}
328
329//_____________________________________________________________________________
330void AliRun::Browse(TBrowser *b)
331{
332 //
333 // Called when the item "Run" is clicked on the left pane
334 // of the Root browser.
335 // It displays the Root Trees and all detectors.
336 //
337 if (fTreeK) b->Add(fTreeK,fTreeK->GetName());
338 if (fTreeH) b->Add(fTreeH,fTreeH->GetName());
339 if (fTreeD) b->Add(fTreeD,fTreeD->GetName());
340 if (fTreeE) b->Add(fTreeE,fTreeE->GetName());
341 if (fTreeR) b->Add(fTreeR,fTreeR->GetName());
342
8494b010 343 TIter next(fModules);
344 AliModule *detector;
345 while((detector = (AliModule*)next())) {
fe4da5cc 346 b->Add(detector,detector->GetName());
347 }
65fb704d 348 b->Add(fMCQA,"AliMCQA");
fe4da5cc 349}
350
351//_____________________________________________________________________________
352void AliRun::Build()
353{
354 //
355 // Initialize Alice geometry
356 // Dummy routine
357 //
358}
359
360//_____________________________________________________________________________
361void AliRun::BuildSimpleGeometry()
362{
363 //
364 // Create a simple TNode geometry used by Root display engine
365 //
366 // Initialise geometry
367 //
368 fGeometry = new TGeometry("AliceGeom","Galice Geometry for Hits");
369 new TMaterial("void","Vacuum",0,0,0); //Everything is void
370 TBRIK *brik = new TBRIK("S_alice","alice volume","void",2000,2000,3000);
371 brik->SetVisibility(0);
372 new TNode("alice","alice","S_alice");
373}
374
375//_____________________________________________________________________________
376void AliRun::CleanDetectors()
377{
378 //
379 // Clean Detectors at the end of event
380 //
8494b010 381 TIter next(fModules);
382 AliModule *detector;
383 while((detector = (AliModule*)next())) {
fe4da5cc 384 detector->FinishEvent();
385 }
386}
387
388//_____________________________________________________________________________
389void AliRun::CleanParents()
390{
391 //
392 // Clean Particles stack.
1578254f 393 // Set parent/daughter relations
fe4da5cc 394 //
395 TClonesArray &particles = *(gAlice->Particles());
1578254f 396 TParticle *part;
fe4da5cc 397 int i;
398 for(i=0; i<fNtrack; i++) {
1578254f 399 part = (TParticle *)particles.UncheckedAt(i);
aee8290b 400 if(!part->TestBit(kDaughtersBit)) {
1578254f 401 part->SetFirstDaughter(-1);
402 part->SetLastDaughter(-1);
fe4da5cc 403 }
404 }
405}
406
407//_____________________________________________________________________________
408Int_t AliRun::DistancetoPrimitive(Int_t, Int_t)
409{
410 //
411 // Return the distance from the mouse to the AliRun object
412 // Dummy routine
413 //
414 return 9999;
415}
416
417//_____________________________________________________________________________
94de3818 418void AliRun::DumpPart (Int_t i) const
fe4da5cc 419{
420 //
421 // Dumps particle i in the stack
422 //
423 TClonesArray &particles = *fParticles;
1578254f 424 ((TParticle*) particles[i])->Print();
fe4da5cc 425}
426
427//_____________________________________________________________________________
94de3818 428void AliRun::DumpPStack () const
fe4da5cc 429{
430 //
431 // Dumps the particle stack
432 //
433 TClonesArray &particles = *fParticles;
434 printf(
435 "\n\n=======================================================================\n");
436 for (Int_t i=0;i<fNtrack;i++)
437 {
1578254f 438 printf("-> %d ",i); ((TParticle*) particles[i])->Print();
fe4da5cc 439 printf("--------------------------------------------------------------\n");
440 }
441 printf(
442 "\n=======================================================================\n\n");
443}
444
d8408e76 445void AliRun::SetField(AliMagF* magField)
446{
447 // Set Magnetic Field Map
448 fField = magField;
449 fField->ReadField();
450}
451
fe4da5cc 452//_____________________________________________________________________________
453void AliRun::SetField(Int_t type, Int_t version, Float_t scale,
454 Float_t maxField, char* filename)
455{
456 //
457 // Set magnetic field parameters
458 // type Magnetic field transport flag 0=no field, 2=helix, 3=Runge Kutta
459 // version Magnetic field map version (only 1 active now)
460 // scale Scale factor for the magnetic field
461 // maxField Maximum value for the magnetic field
462
463 //
464 // --- Sanity check on mag field flags
fe4da5cc 465 if(fField) delete fField;
466 if(version==1) {
d8408e76 467 fField = new AliMagFC("Map1"," ",type,scale,maxField);
f1b9d7c3 468 } else if(version<=2) {
d8408e76 469 fField = new AliMagFCM("Map2-3",filename,type,scale,maxField);
fe4da5cc 470 fField->ReadField();
f1b9d7c3 471 } else if(version==3) {
d8408e76 472 fField = new AliMagFDM("Map4",filename,type,scale,maxField);
f1b9d7c3 473 fField->ReadField();
fe4da5cc 474 } else {
23370b7a 475 Warning("SetField","Invalid map %d\n",version);
fe4da5cc 476 }
477}
478
479//_____________________________________________________________________________
480void AliRun::FillTree()
481{
482 //
483 // Fills all AliRun TTrees
484 //
485 if (fTreeK) fTreeK->Fill();
486 if (fTreeH) fTreeH->Fill();
487 if (fTreeD) fTreeD->Fill();
488 if (fTreeR) fTreeR->Fill();
489}
490
65fb704d 491//_____________________________________________________________________________
492void AliRun::PreTrack()
493{
494 TObjArray &dets = *fModules;
495 AliModule *module;
496
497 for(Int_t i=0; i<=fNdets; i++)
498 if((module = (AliModule*)dets[i]))
499 module->PreTrack();
500
501 fMCQA->PreTrack();
502}
503
504//_____________________________________________________________________________
505void AliRun::PostTrack()
506{
507 TObjArray &dets = *fModules;
508 AliModule *module;
509
510 for(Int_t i=0; i<=fNdets; i++)
511 if((module = (AliModule*)dets[i]))
512 module->PostTrack();
513}
514
fe4da5cc 515//_____________________________________________________________________________
516void AliRun::FinishPrimary()
517{
518 //
519 // Called at the end of each primary track
520 //
521
6c9704e6 522 // static Int_t count=0;
523 // const Int_t times=10;
fe4da5cc 524 // This primary is finished, purify stack
80762cb1 525 PurifyKine();
fe4da5cc 526
24de2263 527 TIter next(fModules);
528 AliModule *detector;
529 while((detector = (AliModule*)next())) {
530 detector->FinishPrimary();
531 }
532
fe4da5cc 533 // Write out hits if any
534 if (gAlice->TreeH()) {
535 gAlice->TreeH()->Fill();
536 }
537
538 // Reset Hits info
539 gAlice->ResetHits();
a8f1fb7c 540
541 //
542 // if(++count%times==1) gObjectTable->Print();
fe4da5cc 543}
544
545//_____________________________________________________________________________
546void AliRun::FinishEvent()
547{
548 //
549 // Called at the end of the event.
550 //
7fb01480 551
dffd31ef 552 //
553 if(fLego) fLego->FinishEvent();
554
fe4da5cc 555 //Update the energy deposit tables
556 Int_t i;
875c717b 557 for(i=0;i<fEventEnergy.GetSize();i++) {
558 fSummEnergy[i]+=fEventEnergy[i];
559 fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
fe4da5cc 560 }
875c717b 561 fEventEnergy.Reset();
fe4da5cc 562
563 // Clean detector information
564 CleanDetectors();
565
566 // Write out the kinematics
567 if (fTreeK) {
568 CleanParents();
569 fTreeK->Fill();
570 }
571
572 // Write out the digits
573 if (fTreeD) {
574 fTreeD->Fill();
575 ResetDigits();
576 }
577
578 // Write out reconstructed clusters
579 if (fTreeR) {
580 fTreeR->Fill();
581 }
582
583 // Write out the event Header information
584 if (fTreeE) fTreeE->Fill();
585
586 // Reset stack info
587 ResetStack();
588
589 // Write Tree headers
59fe9bd2 590 // Int_t ievent = fHeader.GetEvent();
591 // char hname[30];
592 // sprintf(hname,"TreeK%d",ievent);
51e0e89d 593 if (fTreeK) fTreeK->Write(0,TObject::kOverwrite);
59fe9bd2 594 // sprintf(hname,"TreeH%d",ievent);
51e0e89d 595 if (fTreeH) fTreeH->Write(0,TObject::kOverwrite);
59fe9bd2 596 // sprintf(hname,"TreeD%d",ievent);
51e0e89d 597 if (fTreeD) fTreeD->Write(0,TObject::kOverwrite);
59fe9bd2 598 // sprintf(hname,"TreeR%d",ievent);
51e0e89d 599 if (fTreeR) fTreeR->Write(0,TObject::kOverwrite);
875c717b 600
601 ++fEvent;
fe4da5cc 602}
603
604//_____________________________________________________________________________
605void AliRun::FinishRun()
606{
607 //
608 // Called at the end of the run.
609 //
610
dffd31ef 611 //
612 if(fLego) fLego->FinishRun();
613
fe4da5cc 614 // Clean detector information
8494b010 615 TIter next(fModules);
616 AliModule *detector;
617 while((detector = (AliModule*)next())) {
fe4da5cc 618 detector->FinishRun();
619 }
620
621 //Output energy summary tables
622 EnergySummary();
623
624 // file is retrieved from whatever tree
aee8290b 625 TFile *file = 0;
626 if (fTreeK) file = fTreeK->GetCurrentFile();
627 if ((!file) && (fTreeH)) file = fTreeH->GetCurrentFile();
628 if ((!file) && (fTreeD)) file = fTreeD->GetCurrentFile();
629 if ((!file) && (fTreeE)) file = fTreeE->GetCurrentFile();
630 if( NULL==file ) {
fe4da5cc 631 Error("FinishRun","There isn't root file!");
632 exit(1);
633 }
aee8290b 634 file->cd();
51e0e89d 635 fTreeE->Write(0,TObject::kOverwrite);
fe4da5cc 636
637 // Clean tree information
d2ecd553 638 if (fTreeK) {
639 delete fTreeK; fTreeK = 0;
640 }
641 if (fTreeH) {
642 delete fTreeH; fTreeH = 0;
643 }
644 if (fTreeD) {
645 delete fTreeD; fTreeD = 0;
646 }
647 if (fTreeR) {
648 delete fTreeR; fTreeR = 0;
649 }
650 if (fTreeE) {
651 delete fTreeE; fTreeE = 0;
652 }
fe4da5cc 653
654 // Write AliRun info and all detectors parameters
655 Write();
656
657 // Close output file
aee8290b 658 file->Write();
fe4da5cc 659}
660
661//_____________________________________________________________________________
662void AliRun::FlagTrack(Int_t track)
663{
664 //
665 // Flags a track and all its family tree to be kept
666 //
667 int curr;
1578254f 668 TParticle *particle;
fe4da5cc 669
670 curr=track;
671 while(1) {
1578254f 672 particle=(TParticle*)fParticles->UncheckedAt(curr);
fe4da5cc 673
674 // If the particle is flagged the three from here upward is saved already
aee8290b 675 if(particle->TestBit(kKeepBit)) return;
fe4da5cc 676
677 // Save this particle
aee8290b 678 particle->SetBit(kKeepBit);
fe4da5cc 679
680 // Move to father if any
1578254f 681 if((curr=particle->GetFirstMother())==-1) return;
fe4da5cc 682 }
683}
684
685//_____________________________________________________________________________
686void AliRun::EnergySummary()
687{
688 //
689 // Print summary of deposited energy
690 //
691
fe4da5cc 692 Int_t ndep=0;
693 Float_t edtot=0;
694 Float_t ed, ed2;
695 Int_t kn, i, left, j, id;
aee8290b 696 const Float_t kzero=0;
fe4da5cc 697 Int_t ievent=fHeader.GetEvent()+1;
698 //
699 // Energy loss information
700 if(ievent) {
701 printf("***************** Energy Loss Information per event (GEV) *****************\n");
875c717b 702 for(kn=1;kn<fEventEnergy.GetSize();kn++) {
703 ed=fSummEnergy[kn];
fe4da5cc 704 if(ed>0) {
875c717b 705 fEventEnergy[ndep]=kn;
fe4da5cc 706 if(ievent>1) {
707 ed=ed/ievent;
875c717b 708 ed2=fSum2Energy[kn];
fe4da5cc 709 ed2=ed2/ievent;
aee8290b 710 ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
fe4da5cc 711 } else
712 ed2=99;
875c717b 713 fSummEnergy[ndep]=ed;
aee8290b 714 fSum2Energy[ndep]=TMath::Min((Float_t) 99.,TMath::Max(ed2,kzero));
fe4da5cc 715 edtot+=ed;
716 ndep++;
717 }
718 }
719 for(kn=0;kn<(ndep-1)/3+1;kn++) {
720 left=ndep-kn*3;
721 for(i=0;i<(3<left?3:left);i++) {
722 j=kn*3+i;
875c717b 723 id=Int_t (fEventEnergy[j]+0.1);
724 printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
fe4da5cc 725 }
726 printf("\n");
727 }
728 //
729 // Relative energy loss in different detectors
730 printf("******************** Relative Energy Loss per event ********************\n");
731 printf("Total energy loss per event %10.3f GeV\n",edtot);
732 for(kn=0;kn<(ndep-1)/5+1;kn++) {
733 left=ndep-kn*5;
734 for(i=0;i<(5<left?5:left);i++) {
735 j=kn*5+i;
875c717b 736 id=Int_t (fEventEnergy[j]+0.1);
737 printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
fe4da5cc 738 }
739 printf("\n");
740 }
741 for(kn=0;kn<75;kn++) printf("*");
742 printf("\n");
743 }
744 //
745 // Reset the TArray's
875c717b 746 // fEventEnergy.Set(0);
747 // fSummEnergy.Set(0);
748 // fSum2Energy.Set(0);
fe4da5cc 749}
750
751//_____________________________________________________________________________
94de3818 752AliModule *AliRun::GetModule(const char *name) const
fe4da5cc 753{
754 //
755 // Return pointer to detector from name
756 //
8494b010 757 return (AliModule*)fModules->FindObject(name);
fe4da5cc 758}
759
a68348e9 760//_____________________________________________________________________________
94de3818 761AliDetector *AliRun::GetDetector(const char *name) const
a68348e9 762{
763 //
764 // Return pointer to detector from name
765 //
766 return (AliDetector*)fModules->FindObject(name);
767}
768
fe4da5cc 769//_____________________________________________________________________________
94de3818 770Int_t AliRun::GetModuleID(const char *name) const
fe4da5cc 771{
772 //
773 // Return galice internal detector identifier from name
774 //
23370b7a 775 Int_t i=-1;
776 TObject *mod=fModules->FindObject(name);
777 if(mod) i=fModules->IndexOf(mod);
778 return i;
fe4da5cc 779}
780
781//_____________________________________________________________________________
782Int_t AliRun::GetEvent(Int_t event)
783{
784 //
785 // Connect the Trees Kinematics and Hits for event # event
786 // Set branch addresses
787 //
fe4da5cc 788
789 // Reset existing structures
790 ResetStack();
791 ResetHits();
792 ResetDigits();
793
794 // Delete Trees already connected
795 if (fTreeK) delete fTreeK;
796 if (fTreeH) delete fTreeH;
797 if (fTreeD) delete fTreeD;
798 if (fTreeR) delete fTreeR;
59fe9bd2 799
800 // Get header from file
801 if(fTreeE) fTreeE->GetEntry(event);
802 else Error("GetEvent","Cannot file Header Tree\n");
fe4da5cc 803
804 // Get Kine Tree from file
805 char treeName[20];
806 sprintf(treeName,"TreeK%d",event);
807 fTreeK = (TTree*)gDirectory->Get(treeName);
808 if (fTreeK) fTreeK->SetBranchAddress("Particles", &fParticles);
23370b7a 809 else Error("GetEvent","cannot find Kine Tree for event:%d\n",event);
fe4da5cc 810
811 // Get Hits Tree header from file
812 sprintf(treeName,"TreeH%d",event);
813 fTreeH = (TTree*)gDirectory->Get(treeName);
814 if (!fTreeH) {
23370b7a 815 Error("GetEvent","cannot find Hits Tree for event:%d\n",event);
fe4da5cc 816 }
817
818 // Get Digits Tree header from file
819 sprintf(treeName,"TreeD%d",event);
820 fTreeD = (TTree*)gDirectory->Get(treeName);
821 if (!fTreeD) {
07a68c1d 822 Warning("GetEvent","cannot find Digits Tree for event:%d\n",event);
fe4da5cc 823 }
824
825
826 // Get Reconstruct Tree header from file
827 sprintf(treeName,"TreeR%d",event);
828 fTreeR = (TTree*)gDirectory->Get(treeName);
829 if (!fTreeR) {
830 // printf("WARNING: cannot find Reconstructed Tree for event:%d\n",event);
831 }
832
833 // Set Trees branch addresses
8494b010 834 TIter next(fModules);
835 AliModule *detector;
836 while((detector = (AliModule*)next())) {
fe4da5cc 837 detector->SetTreeAddress();
838 }
839
840 if (fTreeK) fTreeK->GetEvent(0);
841 fNtrack = Int_t (fParticles->GetEntries());
842 return fNtrack;
843}
844
845//_____________________________________________________________________________
846TGeometry *AliRun::GetGeometry()
847{
848 //
849 // Import Alice geometry from current file
850 // Return pointer to geometry object
851 //
852 if (!fGeometry) fGeometry = (TGeometry*)gDirectory->Get("AliceGeom");
853 //
854 // Unlink and relink nodes in detectors
855 // This is bad and there must be a better way...
856 //
fe4da5cc 857
8494b010 858 TIter next(fModules);
859 AliModule *detector;
860 while((detector = (AliModule*)next())) {
fe4da5cc 861 detector->SetTreeAddress();
862 TList *dnodes=detector->Nodes();
863 Int_t j;
864 TNode *node, *node1;
865 for ( j=0; j<dnodes->GetSize(); j++) {
866 node = (TNode*) dnodes->At(j);
52d0ab00 867 node1 = fGeometry->GetNode(node->GetName());
fe4da5cc 868 dnodes->Remove(node);
869 dnodes->AddAt(node1,j);
870 }
871 }
872 return fGeometry;
873}
874
875//_____________________________________________________________________________
876void AliRun::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom,
877 Float_t &e, Float_t *vpos, Float_t *polar,
878 Float_t &tof)
879{
880 //
881 // Return next track from stack of particles
882 //
a8f1fb7c 883 TVector3 pol;
fe4da5cc 884 fCurrent=-1;
1578254f 885 TParticle *track;
fe4da5cc 886 for(Int_t i=fNtrack-1; i>=0; i--) {
1578254f 887 track=(TParticle*) fParticles->UncheckedAt(i);
aee8290b 888 if(!track->TestBit(kDoneBit)) {
fe4da5cc 889 //
890 // The track has not yet been processed
891 fCurrent=i;
1578254f 892 ipart=track->GetPdgCode();
893 pmom[0]=track->Px();
894 pmom[1]=track->Py();
895 pmom[2]=track->Pz();
896 e =track->Energy();
897 vpos[0]=track->Vx();
898 vpos[1]=track->Vy();
899 vpos[2]=track->Vz();
a8f1fb7c 900 track->GetPolarisation(pol);
901 polar[0]=pol.X();
902 polar[1]=pol.Y();
903 polar[2]=pol.Z();
1578254f 904 tof=track->T();
aee8290b 905 track->SetBit(kDoneBit);
fe4da5cc 906 break;
907 }
908 }
909 mtrack=fCurrent;
910 //
911 // stop and start timer when we start a primary track
912 Int_t nprimaries = fHeader.GetNprimary();
913 if (fCurrent >= nprimaries) return;
914 if (fCurrent < nprimaries-1) {
915 fTimer.Stop();
1578254f 916 track=(TParticle*) fParticles->UncheckedAt(fCurrent+1);
917 // track->SetProcessTime(fTimer.CpuTime());
fe4da5cc 918 }
919 fTimer.Start();
920}
921
922//_____________________________________________________________________________
94de3818 923Int_t AliRun::GetPrimary(Int_t track) const
fe4da5cc 924{
925 //
926 // return number of primary that has generated track
927 //
928 int current, parent;
1578254f 929 TParticle *part;
fe4da5cc 930 //
931 parent=track;
932 while (1) {
933 current=parent;
1578254f 934 part = (TParticle *)fParticles->UncheckedAt(current);
935 parent=part->GetFirstMother();
fe4da5cc 936 if(parent<0) return current;
937 }
938}
939
940//_____________________________________________________________________________
875c717b 941void AliRun::InitMC(const char *setup)
fe4da5cc 942{
943 //
944 // Initialize the Alice setup
945 //
946
51e0e89d 947 if(fInitDone) {
948 Warning("Init","Cannot initialise AliRun twice!\n");
949 return;
950 }
951
fe4da5cc 952 gROOT->LoadMacro(setup);
45189757 953 gInterpreter->ProcessLine(fConfigFunction.Data());
fe4da5cc 954
cfce8870 955 gMC->DefineParticles(); //Create standard MC particles
fe4da5cc 956
957 TObject *objfirst, *objlast;
958
23370b7a 959 fNdets = fModules->GetLast()+1;
960
fe4da5cc 961 //
875c717b 962 //=================Create Materials and geometry
963 gMC->Init();
964
65fb704d 965 // Added also after in case of interactive initialisation of modules
966 fNdets = fModules->GetLast()+1;
967
8494b010 968 TIter next(fModules);
969 AliModule *detector;
970 while((detector = (AliModule*)next())) {
fe4da5cc 971 detector->SetTreeAddress();
972 objlast = gDirectory->GetList()->Last();
973
fe4da5cc 974 // Add Detector histograms in Detector list of histograms
975 if (objlast) objfirst = gDirectory->GetList()->After(objlast);
976 else objfirst = gDirectory->GetList()->First();
977 while (objfirst) {
978 detector->Histograms()->Add(objfirst);
979 objfirst = gDirectory->GetList()->After(objfirst);
980 }
981 }
327136d2 982 ReadTransPar(); //Read the cuts for all materials
fe4da5cc 983
984 MediaTable(); //Build the special IMEDIA table
985
fe4da5cc 986 //Initialise geometry deposition table
875c717b 987 fEventEnergy.Set(gMC->NofVolumes()+1);
988 fSummEnergy.Set(gMC->NofVolumes()+1);
989 fSum2Energy.Set(gMC->NofVolumes()+1);
fe4da5cc 990
fe4da5cc 991 //Compute cross-sections
875c717b 992 gMC->BuildPhysics();
fe4da5cc 993
994 //Write Geometry object to current file.
995 fGeometry->Write();
996
997 fInitDone = kTRUE;
51e0e89d 998
65fb704d 999 fMCQA = new AliMCQA(fNdets);
1000
51e0e89d 1001 //
1002 // Save stuff at the beginning of the file to avoid file corruption
1003 Write();
fe4da5cc 1004}
1005
1006//_____________________________________________________________________________
1007void AliRun::MediaTable()
1008{
1009 //
1010 // Built media table to get from the media number to
1011 // the detector id
1012 //
ad51aeb0 1013 Int_t kz, nz, idt, lz, i, k, ind;
1014 // Int_t ibeg;
fe4da5cc 1015 TObjArray &dets = *gAlice->Detectors();
8494b010 1016 AliModule *det;
fe4da5cc 1017 //
1018 // For all detectors
1019 for (kz=0;kz<fNdets;kz++) {
1020 // If detector is defined
8494b010 1021 if((det=(AliModule*) dets[kz])) {
ad51aeb0 1022 TArrayI &idtmed = *(det->GetIdtmed());
1023 for(nz=0;nz<100;nz++) {
fe4da5cc 1024 // Find max and min material number
ad51aeb0 1025 if((idt=idtmed[nz])) {
fe4da5cc 1026 det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
1027 det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
1028 }
1029 }
1030 if(det->LoMedium() > det->HiMedium()) {
1031 det->LoMedium() = 0;
1032 det->HiMedium() = 0;
1033 } else {
1034 if(det->HiMedium() > fImedia->GetSize()) {
ad51aeb0 1035 Error("MediaTable","Increase fImedia from %d to %d",
1036 fImedia->GetSize(),det->HiMedium());
fe4da5cc 1037 return;
1038 }
1039 // Tag all materials in rage as belonging to detector kz
1040 for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
1041 (*fImedia)[lz]=kz;
1042 }
1043 }
1044 }
1045 }
1046 //
1047 // Print summary table
1048 printf(" Traking media ranges:\n");
1049 for(i=0;i<(fNdets-1)/6+1;i++) {
1050 for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
1051 ind=i*6+k;
8494b010 1052 det=(AliModule*)dets[ind];
fe4da5cc 1053 if(det)
1054 printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
1055 det->HiMedium());
1056 else
1057 printf(" %6s: %3d -> %3d;","NULL",0,0);
1058 }
1059 printf("\n");
1060 }
1061}
1062
1063//____________________________________________________________________________
1064void AliRun::SetGenerator(AliGenerator *generator)
ee1dd322 1065{
1066 //
1067 // Load the event generator
1068 //
1069 if(!fGenerator) fGenerator = generator;
1070}
1071
1072//____________________________________________________________________________
1073void AliRun::ResetGenerator(AliGenerator *generator)
fe4da5cc 1074{
1075 //
1076 // Load the event generator
1077 //
b13db077 1078 if(fGenerator)
838edcaf 1079 if(generator)
1080 Warning("ResetGenerator","Replacing generator %s with %s\n",
1081 fGenerator->GetName(),generator->GetName());
1082 else
1083 Warning("ResetGenerator","Replacing generator %s with NULL\n",
1084 fGenerator->GetName());
b13db077 1085 fGenerator = generator;
fe4da5cc 1086}
1087
1088//____________________________________________________________________________
65fb704d 1089void AliRun::SetTransPar(char *filename)
1090{
1091 fTransParName = filename;
1092}
1093
1094//____________________________________________________________________________
1095void AliRun::ReadTransPar()
fe4da5cc 1096{
1097 //
1098 // Read filename to set the transport parameters
1099 //
1100
fe4da5cc 1101
aee8290b 1102 const Int_t kncuts=10;
1103 const Int_t knflags=11;
1104 const Int_t knpars=kncuts+knflags;
1105 const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
fe4da5cc 1106 "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
1107 "BREM","COMP","DCAY","DRAY","HADR","LOSS",
1108 "MULS","PAIR","PHOT","RAYL"};
1109 char line[256];
ad51aeb0 1110 char detName[7];
fe4da5cc 1111 char* filtmp;
aee8290b 1112 Float_t cut[kncuts];
1113 Int_t flag[knflags];
fe4da5cc 1114 Int_t i, itmed, iret, ktmed, kz;
1115 FILE *lun;
1116 //
1117 // See whether the file is there
65fb704d 1118 filtmp=gSystem->ExpandPathName(fTransParName.Data());
fe4da5cc 1119 lun=fopen(filtmp,"r");
1120 delete [] filtmp;
1121 if(!lun) {
65fb704d 1122 Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
fe4da5cc 1123 return;
1124 }
1125 //
1126 printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1127 printf(" *%59s\n","*");
1128 printf(" * Please check carefully what you are doing!%10s\n","*");
1129 printf(" *%59s\n","*");
1130 //
1131 while(1) {
1132 // Initialise cuts and flags
aee8290b 1133 for(i=0;i<kncuts;i++) cut[i]=-99;
1134 for(i=0;i<knflags;i++) flag[i]=-99;
fe4da5cc 1135 itmed=0;
1136 for(i=0;i<256;i++) line[i]='\0';
1137 // Read up to the end of line excluded
1138 iret=fscanf(lun,"%[^\n]",line);
1139 if(iret<0) {
1140 //End of file
1141 fclose(lun);
1142 printf(" *%59s\n","*");
1143 printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
1144 return;
1145 }
1146 // Read the end of line
1147 fscanf(lun,"%*c");
1148 if(!iret) continue;
1149 if(line[0]=='*') continue;
1150 // Read the numbers
ad51aeb0 1151 iret=sscanf(line,"%s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d",
1152 detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
1153 &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
1154 &flag[8],&flag[9],&flag[10]);
fe4da5cc 1155 if(!iret) continue;
1156 if(iret<0) {
1157 //reading error
65fb704d 1158 Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
fe4da5cc 1159 continue;
1160 }
ad51aeb0 1161 // Check that the module exist
1162 AliModule *mod = GetModule(detName);
1163 if(mod) {
1164 // Get the array of media numbers
1165 TArrayI &idtmed = *mod->GetIdtmed();
1166 // Check that the tracking medium code is valid
1167 if(0<=itmed && itmed < 100) {
1168 ktmed=idtmed[itmed];
1169 if(!ktmed) {
65fb704d 1170 Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
ad51aeb0 1171 continue;
fe4da5cc 1172 }
ad51aeb0 1173 // Set energy thresholds
aee8290b 1174 for(kz=0;kz<kncuts;kz++) {
ad51aeb0 1175 if(cut[kz]>=0) {
23370b7a 1176 printf(" * %-6s set to %10.3E for tracking medium code %4d for %s\n",
aee8290b 1177 kpars[kz],cut[kz],itmed,mod->GetName());
1178 gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
ad51aeb0 1179 }
fe4da5cc 1180 }
ad51aeb0 1181 // Set transport mechanisms
aee8290b 1182 for(kz=0;kz<knflags;kz++) {
ad51aeb0 1183 if(flag[kz]>=0) {
1184 printf(" * %-6s set to %10d for tracking medium code %4d for %s\n",
aee8290b 1185 kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
1186 gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
ad51aeb0 1187 }
1188 }
1189 } else {
65fb704d 1190 Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
ad51aeb0 1191 continue;
fe4da5cc 1192 }
1193 } else {
65fb704d 1194 Warning("ReadTransPar","Module %s not present\n",detName);
fe4da5cc 1195 continue;
1196 }
1197 }
1198}
1199
1200//_____________________________________________________________________________
1201void AliRun::MakeTree(Option_t *option)
1202{
1203 //
1204 // Create the ROOT trees
1205 // Loop on all detectors to create the Root branch (if any)
1206 //
1207
b13db077 1208 char hname[30];
fe4da5cc 1209 //
1210 // Analyse options
aee8290b 1211 char *oK = strstr(option,"K");
1212 char *oH = strstr(option,"H");
1213 char *oE = strstr(option,"E");
1214 char *oD = strstr(option,"D");
1215 char *oR = strstr(option,"R");
fe4da5cc 1216 //
aee8290b 1217 if (oK && !fTreeK) {
b13db077 1218 sprintf(hname,"TreeK%d",fEvent);
1219 fTreeK = new TTree(hname,"Kinematics");
1220 // Create a branch for particles
1221 fTreeK->Branch("Particles",&fParticles,4000);
51e0e89d 1222 fTreeK->Write();
b13db077 1223 }
aee8290b 1224 if (oH && !fTreeH) {
b13db077 1225 sprintf(hname,"TreeH%d",fEvent);
1226 fTreeH = new TTree(hname,"Hits");
1227 fTreeH->SetAutoSave(1000000000); //no autosave
51e0e89d 1228 fTreeH->Write();
b13db077 1229 }
aee8290b 1230 if (oD && !fTreeD) {
b13db077 1231 sprintf(hname,"TreeD%d",fEvent);
1232 fTreeD = new TTree(hname,"Digits");
51e0e89d 1233 fTreeD->Write();
b13db077 1234 }
aee8290b 1235 if (oR && !fTreeR) {
b13db077 1236 sprintf(hname,"TreeR%d",fEvent);
1237 fTreeR = new TTree(hname,"Reconstruction");
51e0e89d 1238 fTreeR->Write();
b13db077 1239 }
aee8290b 1240 if (oE && !fTreeE) {
b13db077 1241 fTreeE = new TTree("TE","Header");
1242 // Create a branch for Header
aee8290b 1243 fTreeE->Branch("Header","AliHeader",&gAliHeader,4000);
51e0e89d 1244 fTreeE->Write();
b13db077 1245 }
fe4da5cc 1246 //
1247 // Create a branch for hits/digits for each detector
1248 // Each branch is a TClonesArray. Each data member of the Hits classes
1249 // will be in turn a subbranch of the detector master branch
8494b010 1250 TIter next(fModules);
1251 AliModule *detector;
1252 while((detector = (AliModule*)next())) {
aee8290b 1253 if (oH || oD || oR) detector->MakeBranch(option);
fe4da5cc 1254 }
fe4da5cc 1255}
1256
1257//_____________________________________________________________________________
1258Int_t AliRun::PurifyKine(Int_t lastSavedTrack, Int_t nofTracks)
1259{
1260 //
1261 // PurifyKine with external parameters
1262 //
1263 fHgwmk = lastSavedTrack;
1264 fNtrack = nofTracks;
1265 PurifyKine();
1266 return fHgwmk;
1267}
1268
1269//_____________________________________________________________________________
1270void AliRun::PurifyKine()
1271{
1272 //
1273 // Compress kinematic tree keeping only flagged particles
1274 // and renaming the particle id's in all the hits
1275 //
1276 TClonesArray &particles = *fParticles;
1277 int nkeep=fHgwmk+1, parent, i;
1578254f 1278 TParticle *part, *partnew, *father;
fe4da5cc 1279 int *map = new int[particles.GetEntries()];
1280
1281 // Save in Header total number of tracks before compression
1282 fHeader.SetNtrack(fHeader.GetNtrack()+fNtrack-fHgwmk);
1283
a01a8b12 1284 // First pass, invalid Daughter information
fe4da5cc 1285 for(i=0; i<fNtrack; i++) {
a01a8b12 1286 // Preset map, to be removed later
1287 if(i<=fHgwmk) map[i]=i ; else map[i] = -99;
aee8290b 1288 ((TParticle *)particles.UncheckedAt(i))->ResetBit(kDaughtersBit);
a01a8b12 1289 }
fe4da5cc 1290 // Second pass, build map between old and new numbering
1291 for(i=fHgwmk+1; i<fNtrack; i++) {
1578254f 1292 part = (TParticle *)particles.UncheckedAt(i);
aee8290b 1293 if(part->TestBit(kKeepBit)) {
fe4da5cc 1294
1295 // This particle has to be kept
1296 map[i]=nkeep;
1297 if(i!=nkeep) {
1298
1299 // Old and new are different, have to copy
1578254f 1300 partnew = (TParticle *)particles.UncheckedAt(nkeep);
5eb58812 1301 // Change due to a bug in the HP compiler
1302 // *partnew = *part;
1303 memcpy(partnew,part,sizeof(TParticle));
fe4da5cc 1304 } else partnew = part;
1305
1306 // as the parent is always *before*, it must be already
1307 // in place. This is what we are checking anyway!
1578254f 1308 if((parent=partnew->GetFirstMother())>fHgwmk) {
fe4da5cc 1309 if(map[parent]==-99) printf("map[%d] = -99!\n",parent);
1578254f 1310 partnew->SetFirstMother(map[parent]);
fe4da5cc 1311 }
1312 nkeep++;
1313 }
1314 }
1315 fNtrack=nkeep;
1316
1578254f 1317 // Fix daughters information
ae23d366 1318 for (i=0; i<fNtrack; i++) {
1578254f 1319 part = (TParticle *)particles.UncheckedAt(i);
1320 parent = part->GetFirstMother();
ae23d366 1321 if(parent>=0) {
1322 father = (TParticle *)particles.UncheckedAt(parent);
aee8290b 1323 if(father->TestBit(kDaughtersBit)) {
fe4da5cc 1324
ae23d366 1325 if(i<father->GetFirstDaughter()) father->SetFirstDaughter(i);
1326 if(i>father->GetLastDaughter()) father->SetLastDaughter(i);
1327 } else {
1328 // Iitialise daughters info for first pass
1329 father->SetFirstDaughter(i);
1330 father->SetLastDaughter(i);
aee8290b 1331 father->SetBit(kDaughtersBit);
ae23d366 1332 }
fe4da5cc 1333 }
1334 }
1335
1cedd08a 1336#ifdef old
fe4da5cc 1337 // Now loop on all detectors and reset the hits
b13db077 1338 AliHit *OneHit;
8494b010 1339 TIter next(fModules);
1340 AliModule *detector;
1341 while((detector = (AliModule*)next())) {
fe4da5cc 1342 if (!detector->Hits()) continue;
1343 TClonesArray &vHits=*(detector->Hits());
1344 if(vHits.GetEntries() != detector->GetNhits())
1345 printf("vHits.GetEntries()!=detector->GetNhits(): %d != %d\n",
1346 vHits.GetEntries(),detector->GetNhits());
1347 for (i=0; i<detector->GetNhits(); i++) {
1348 OneHit = (AliHit *)vHits.UncheckedAt(i);
1349 OneHit->SetTrack(map[OneHit->GetTrack()]);
1350 }
1351 }
1cedd08a 1352#else
1353
1354 // Now loop on all registered hit lists
1355 TIter next(fHitLists);
1356 TCollection *hitList;
1357 while((hitList = (TCollection*)next())) {
1358 TIter nexthit(hitList);
1359 AliHit *hit;
1360 while((hit = (AliHit*)nexthit())) {
1361 hit->SetTrack(map[hit->GetTrack()]);
1362 }
1363 }
1364#endif
fe4da5cc 1365
24de2263 1366 //
1367 // This for detectors which have a special mapping mechanism
1368 // for hits, such as TPC and TRD
1369 //
1370
1371 TIter nextmod(fModules);
1372 AliModule *detector;
1373 while((detector = (AliModule*)nextmod())) {
1374 detector->RemapTrackHitIDs(map);
1375 }
1376
1377
fe4da5cc 1378 fHgwmk=nkeep-1;
1379 particles.SetLast(fHgwmk);
1380 delete [] map;
1381}
1382
1383//_____________________________________________________________________________
dffd31ef 1384void AliRun::BeginEvent()
fe4da5cc 1385{
1386 //
1387 // Reset all Detectors & kinematics & trees
1388 //
59fe9bd2 1389 char hname[30];
dffd31ef 1390 //
1391
1392 //
1393 if(fLego) {
1394 fLego->BeginEvent();
1395 return;
1396 }
1397
59fe9bd2 1398 //
fe4da5cc 1399 ResetStack();
1400 ResetHits();
1401 ResetDigits();
1402
1403 // Initialise event header
875c717b 1404 fHeader.Reset(fRun,fEvent);
fe4da5cc 1405
59fe9bd2 1406 if(fTreeK) {
1407 fTreeK->Reset();
875c717b 1408 sprintf(hname,"TreeK%d",fEvent);
59fe9bd2 1409 fTreeK->SetName(hname);
1410 }
1411 if(fTreeH) {
1412 fTreeH->Reset();
875c717b 1413 sprintf(hname,"TreeH%d",fEvent);
59fe9bd2 1414 fTreeH->SetName(hname);
1415 }
1416 if(fTreeD) {
1417 fTreeD->Reset();
875c717b 1418 sprintf(hname,"TreeD%d",fEvent);
59fe9bd2 1419 fTreeD->SetName(hname);
1420 }
1421 if(fTreeR) {
1422 fTreeR->Reset();
875c717b 1423 sprintf(hname,"TreeR%d",fEvent);
59fe9bd2 1424 fTreeR->SetName(hname);
1425 }
fe4da5cc 1426}
1427
1428//_____________________________________________________________________________
1429void AliRun::ResetDigits()
1430{
1431 //
1432 // Reset all Detectors digits
1433 //
8494b010 1434 TIter next(fModules);
1435 AliModule *detector;
1436 while((detector = (AliModule*)next())) {
fe4da5cc 1437 detector->ResetDigits();
1438 }
1439}
1440
1441//_____________________________________________________________________________
1442void AliRun::ResetHits()
1443{
1444 //
1445 // Reset all Detectors hits
1446 //
8494b010 1447 TIter next(fModules);
1448 AliModule *detector;
1449 while((detector = (AliModule*)next())) {
fe4da5cc 1450 detector->ResetHits();
1451 }
1452}
1453
1454//_____________________________________________________________________________
1455void AliRun::ResetPoints()
1456{
1457 //
1458 // Reset all Detectors points
1459 //
8494b010 1460 TIter next(fModules);
1461 AliModule *detector;
1462 while((detector = (AliModule*)next())) {
fe4da5cc 1463 detector->ResetPoints();
1464 }
1465}
1466
1467//_____________________________________________________________________________
875c717b 1468void AliRun::RunMC(Int_t nevent, const char *setup)
fe4da5cc 1469{
1470 //
1471 // Main function to be called to process a galice run
1472 // example
1473 // Root > gAlice.Run();
1474 // a positive number of events will cause the finish routine
1475 // to be called
1476 //
1477
fe4da5cc 1478 // check if initialisation has been done
875c717b 1479 if (!fInitDone) InitMC(setup);
fe4da5cc 1480
fe4da5cc 1481 // Create the Root Tree with one branch per detector
80762cb1 1482 MakeTree("KHDER");
fe4da5cc 1483
875c717b 1484 gMC->ProcessRun(nevent);
1485
fe4da5cc 1486 // End of this run, close files
80762cb1 1487 if(nevent>0) FinishRun();
fe4da5cc 1488}
1489
1490//_____________________________________________________________________________
0a520a66 1491void AliRun::RunLego(const char *setup, Int_t nc1, Float_t c1min,
1492 Float_t c1max,Int_t nc2,Float_t c2min,Float_t c2max,
1493 Float_t rmin,Float_t rmax,Float_t zmax, AliLegoGenerator* gener)
fe4da5cc 1494{
1495 //
1496 // Generates lego plots of:
1497 // - radiation length map phi vs theta
1498 // - radiation length map phi vs eta
1499 // - interaction length map
1500 // - g/cm2 length map
1501 //
1502 // ntheta bins in theta, eta
1503 // themin minimum angle in theta (degrees)
1504 // themax maximum angle in theta (degrees)
1505 // nphi bins in phi
1506 // phimin minimum angle in phi (degrees)
1507 // phimax maximum angle in phi (degrees)
1508 // rmin minimum radius
1509 // rmax maximum radius
1510 //
1511 //
1512 // The number of events generated = ntheta*nphi
1513 // run input parameters in macro setup (default="Config.C")
1514 //
1515 // Use macro "lego.C" to visualize the 3 lego plots in spherical coordinates
1516 //Begin_Html
1517 /*
1439f98e 1518 <img src="picts/AliRunLego1.gif">
fe4da5cc 1519 */
1520 //End_Html
1521 //Begin_Html
1522 /*
1439f98e 1523 <img src="picts/AliRunLego2.gif">
fe4da5cc 1524 */
1525 //End_Html
1526 //Begin_Html
1527 /*
1439f98e 1528 <img src="picts/AliRunLego3.gif">
fe4da5cc 1529 */
1530 //End_Html
1531 //
1532
1533 // check if initialisation has been done
875c717b 1534 if (!fInitDone) InitMC(setup);
838edcaf 1535 //Save current generator
1536 AliGenerator *gen=Generator();
1537
0a520a66 1538 // Set new generator
1539 if (!gener) gener = new AliLegoGenerator();
1540 ResetGenerator(gener);
1541 //
1542 // Configure Generator
1543 gener->SetRadiusRange(rmin, rmax);
1544 gener->SetZMax(zmax);
1545 gener->SetCoor1Range(nc1, c1min, c1max);
1546 gener->SetCoor2Range(nc2, c2min, c2max);
1547
1548
b13db077 1549 //Create Lego object
0a520a66 1550 fLego = new AliLego("lego",gener);
b13db077 1551
dffd31ef 1552 //Prepare MC for Lego Run
1553 gMC->InitLego();
1554
b13db077 1555 //Run Lego Object
0a520a66 1556
1557 gMC->ProcessRun(nc1*nc2+1);
fe4da5cc 1558
1559 // Create only the Root event Tree
80762cb1 1560 MakeTree("E");
fe4da5cc 1561
1562 // End of this run, close files
80762cb1 1563 FinishRun();
0a520a66 1564 // Restore current generator
1565 ResetGenerator(gen);
838edcaf 1566 // Delete Lego Object
1567 delete fLego; fLego=0;
fe4da5cc 1568}
1569
94de3818 1570//_____________________________________________________________________________
1571void AliRun::SetConfigFunction(const char * config)
1572{
1573 //
1574 // Set the signature of the function contained in Config.C to configure
1575 // the run
1576 //
1577 fConfigFunction=config;
1578}
1579
fe4da5cc 1580//_____________________________________________________________________________
1581void AliRun::SetCurrentTrack(Int_t track)
1582{
1583 //
1584 // Set current track number
1585 //
1586 fCurrent = track;
1587}
1588
1589//_____________________________________________________________________________
1578254f 1590void AliRun::SetTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom,
fe4da5cc 1591 Float_t *vpos, Float_t *polar, Float_t tof,
65fb704d 1592 AliMCProcess mech, Int_t &ntr, Float_t weight)
fe4da5cc 1593{
1594 //
1595 // Load a track on the stack
1596 //
1597 // done 0 if the track has to be transported
1598 // 1 if not
1599 // parent identifier of the parent track. -1 for a primary
1578254f 1600 // pdg particle code
fe4da5cc 1601 // pmom momentum GeV/c
1602 // vpos position
1603 // polar polarisation
1604 // tof time of flight in seconds
1605 // mecha production mechanism
1606 // ntr on output the number of the track stored
1607 //
1608 TClonesArray &particles = *fParticles;
1578254f 1609 TParticle *particle;
fe4da5cc 1610 Float_t mass;
aee8290b 1611 const Int_t kfirstdaughter=-1;
1612 const Int_t klastdaughter=-1;
1613 const Int_t kS=0;
1578254f 1614 // const Float_t tlife=0;
fe4da5cc 1615
1578254f 1616 //
1617 // Here we get the static mass
1618 // For MC is ok, but a more sophisticated method could be necessary
1619 // if the calculated mass is required
1620 // also, this method is potentially dangerous if the mass
1621 // used in the MC is not the same of the PDG database
1622 //
1623 mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
fe4da5cc 1624 Float_t e=TMath::Sqrt(mass*mass+pmom[0]*pmom[0]+
1625 pmom[1]*pmom[1]+pmom[2]*pmom[2]);
1626
aee8290b 1627 //printf("Loading particle %s mass %f ene %f No %d ip %d pos %f %f %f mom %f %f %f kS %d m %s\n",
1628 //pname,mass,e,fNtrack,pdg,vpos[0],vpos[1],vpos[2],pmom[0],pmom[1],pmom[2],kS,mecha);
fe4da5cc 1629
aee8290b 1630 particle=new(particles[fNtrack]) TParticle(pdg,kS,parent,-1,kfirstdaughter,
1631 klastdaughter,pmom[0],pmom[1],pmom[2],
1578254f 1632 e,vpos[0],vpos[1],vpos[2],tof);
1578254f 1633 ((TParticle*)particles[fNtrack])->SetPolarisation(TVector3(polar[0],polar[1],polar[2]));
1634 ((TParticle*)particles[fNtrack])->SetWeight(weight);
65fb704d 1635 particle->SetUniqueID(mech);
aee8290b 1636 if(!done) particle->SetBit(kDoneBit);
a01a8b12 1637 //Declare that the daughter information is valid
aee8290b 1638 ((TParticle*)particles[fNtrack])->SetBit(kDaughtersBit);
fe4da5cc 1639
1640 if(parent>=0) {
1578254f 1641 particle=(TParticle*) fParticles->UncheckedAt(parent);
1642 particle->SetLastDaughter(fNtrack);
1643 if(particle->GetFirstDaughter()<0) particle->SetFirstDaughter(fNtrack);
fe4da5cc 1644 } else {
1645 //
1646 // This is a primary track. Set high water mark for this event
1647 fHgwmk=fNtrack;
1648 //
1649 // Set also number if primary tracks
1650 fHeader.SetNprimary(fHgwmk+1);
1651 fHeader.SetNtrack(fHgwmk+1);
1652 }
1653 ntr = fNtrack++;
1654}
1655
1656//_____________________________________________________________________________
1657void AliRun::KeepTrack(const Int_t track)
1658{
1659 //
1660 // flags a track to be kept
1661 //
1662 TClonesArray &particles = *fParticles;
aee8290b 1663 ((TParticle*)particles[track])->SetBit(kKeepBit);
fe4da5cc 1664}
1665
1666//_____________________________________________________________________________
875c717b 1667void AliRun::StepManager(Int_t id)
fe4da5cc 1668{
1669 //
1670 // Called at every step during transport
1671 //
1672
fe4da5cc 1673 //
1674 // --- If lego option, do it and leave
b13db077 1675 if (fLego)
fe4da5cc 1676 fLego->StepManager();
b13db077 1677 else {
1678 Int_t copy;
1679 //Update energy deposition tables
875c717b 1680 AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep());
fe4da5cc 1681
b13db077 1682 //Call the appropriate stepping routine;
1683 AliModule *det = (AliModule*)fModules->At(id);
65fb704d 1684 if(det) {
1685 fMCQA->StepManager(id);
1686 det->StepManager();
1687 }
fe4da5cc 1688 }
fe4da5cc 1689}
1690
fe4da5cc 1691//_____________________________________________________________________________
1692void AliRun::Streamer(TBuffer &R__b)
1693{
1694 //
1695 // Stream an object of class AliRun.
1696 //
1697 if (R__b.IsReading()) {
1698 Version_t R__v = R__b.ReadVersion(); if (R__v) { }
1699 TNamed::Streamer(R__b);
1700 if (!gAlice) gAlice = this;
1701 gROOT->GetListOfBrowsables()->Add(this,"Run");
59fe9bd2 1702 fTreeE = (TTree*)gDirectory->Get("TE");
aee8290b 1703 if (fTreeE) fTreeE->SetBranchAddress("Header", &gAliHeader);
59fe9bd2 1704 else Error("Streamer","cannot find Header Tree\n");
fe4da5cc 1705 R__b >> fNtrack;
1706 R__b >> fHgwmk;
1707 R__b >> fDebug;
1708 fHeader.Streamer(R__b);
8494b010 1709 R__b >> fModules;
fe4da5cc 1710 R__b >> fParticles;
1711 R__b >> fField;
875c717b 1712 // R__b >> fMC;
fe4da5cc 1713 R__b >> fNdets;
1714 R__b >> fTrRmax;
1715 R__b >> fTrZmax;
1716 R__b >> fGenerator;
59fe9bd2 1717 if(R__v>1) {
1718 R__b >> fPDGDB; //Particle factory object!
1719 fTreeE->GetEntry(0);
1720 } else {
1721 fHeader.SetEvent(0);
1722 fPDGDB = TDatabasePDG::Instance(); //Particle factory object!
1723 }
45189757 1724 if(R__v>2) {
1725 fConfigFunction.Streamer(R__b);
1726 } else {
1727 fConfigFunction="Config();";
1728 }
65fb704d 1729 if(R__v>3) {
1730 R__b >> fRandom;
1731 gRandom = fRandom;
1732 R__b >> fMCQA;
1733 } else {
1734 fRandom = gRandom = new TRandom3();
1735 fMCQA = new AliMCQA();
1736 }
fe4da5cc 1737 } else {
1738 R__b.WriteVersion(AliRun::IsA());
1739 TNamed::Streamer(R__b);
1740 R__b << fNtrack;
1741 R__b << fHgwmk;
1742 R__b << fDebug;
1743 fHeader.Streamer(R__b);
8494b010 1744 R__b << fModules;
fe4da5cc 1745 R__b << fParticles;
1746 R__b << fField;
875c717b 1747 // R__b << fMC;
fe4da5cc 1748 R__b << fNdets;
1749 R__b << fTrRmax;
1750 R__b << fTrZmax;
1751 R__b << fGenerator;
1578254f 1752 R__b << fPDGDB; //Particle factory object!
45189757 1753 fConfigFunction.Streamer(R__b);
65fb704d 1754 R__b << fRandom;
1755 R__b << fMCQA;
fe4da5cc 1756 }
1757}