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