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