- Warning("SetField","Invalid map %d\n",version);
- }
-}
-
-//_____________________________________________________________________________
-void AliRun::PreTrack()
-{
- TObjArray &dets = *fModules;
- AliModule *module;
-
- for(Int_t i=0; i<=fNdets; i++)
- if((module = (AliModule*)dets[i]))
- module->PreTrack();
-
- fMCQA->PreTrack();
-}
-
-//_____________________________________________________________________________
-void AliRun::PostTrack()
-{
- TObjArray &dets = *fModules;
- AliModule *module;
-
- for(Int_t i=0; i<=fNdets; i++)
- if((module = (AliModule*)dets[i]))
- module->PostTrack();
-}
-
-//_____________________________________________________________________________
-void AliRun::FinishPrimary()
-{
- //
- // Called at the end of each primary track
- //
-
- // static Int_t count=0;
- // const Int_t times=10;
- // This primary is finished, purify stack
- fStack->PurifyKine();
-
- TIter next(fModules);
- AliModule *detector;
- while((detector = (AliModule*)next())) {
- detector->FinishPrimary();
- }
-
- // Write out hits if any
- if (gAlice->TreeH()) {
- gAlice->TreeH()->Fill();
- }
-
- //
- // if(++count%times==1) gObjectTable->Print();
-}
-
-//_____________________________________________________________________________
-void AliRun::BeginPrimary()
-{
- //
- // Called at the beginning of each primary track
- //
-
- // Reset Hits info
- gAlice->ResetHits();
-
-}
-
-//_____________________________________________________________________________
-void AliRun::FinishEvent()
-{
- //
- // Called at the end of the event.
- //
-
- //
- if(fLego) fLego->FinishEvent();
-
- //Update the energy deposit tables
- Int_t i;
- for(i=0;i<fEventEnergy.GetSize();i++) {
- fSummEnergy[i]+=fEventEnergy[i];
- fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i];
- }
-
-
-
- // Update Header information
-
- fHeader->SetNprimary(fStack->GetNprimary());
- fHeader->SetNtrack(fStack->GetNtrack());
-
-
- // Write out the kinematics
- fStack->FinishEvent();
-
- // Write out the event Header information
- if (fTreeE) {
- fHeader->SetStack(fStack);
- fTreeE->Fill();
- }
-
-
- // Write Tree headers
- TTree* pTreeK = fStack->TreeK();
- if (pTreeK) pTreeK->Write(0,TObject::kOverwrite);
- if (fTreeH) fTreeH->Write(0,TObject::kOverwrite);
-
- ++fEvent;
-}
-
-//_____________________________________________________________________________
-void AliRun::FinishRun()
-{
- //
- // Called at the end of the run.
- //
-
- //
- if(fLego) fLego->FinishRun();
-
- // Clean detector information
- TIter next(fModules);
- AliModule *detector;
- while((detector = (AliModule*)next())) {
- detector->FinishRun();
- }
-
- //Output energy summary tables
- EnergySummary();
-
- TFile *file = fTreeE->GetCurrentFile();
-
- file->cd();
-
- fTreeE->Write(0,TObject::kOverwrite);
-
- // Write AliRun info and all detectors parameters
- Write(0,TObject::kOverwrite);
-
- // Clean tree information
-
- fStack->FinishRun();
-
- if (fTreeH) {
- delete fTreeH; fTreeH = 0;
- }
- if (fTreeD) {
- delete fTreeD; fTreeD = 0;
- }
- if (fTreeR) {
- delete fTreeR; fTreeR = 0;
- }
- if (fTreeE) {
- delete fTreeE; fTreeE = 0;
- }
- if (fTreeS) {
- delete fTreeS; fTreeS = 0;
- }
-
- // Close output file
- file->Write();
-}
-
-//_____________________________________________________________________________
-void AliRun::FlagTrack(Int_t track)
-{
- // Delegate to stack
- //
- fStack->FlagTrack(track);
-}
-
-//_____________________________________________________________________________
-void AliRun::EnergySummary()
-{
- //
- // Print summary of deposited energy
- //
-
- Int_t ndep=0;
- Float_t edtot=0;
- Float_t ed, ed2;
- Int_t kn, i, left, j, id;
- const Float_t kzero=0;
- Int_t ievent=fHeader->GetEvent()+1;
- //
- // Energy loss information
- if(ievent) {
- printf("***************** Energy Loss Information per event (GEV) *****************\n");
- for(kn=1;kn<fEventEnergy.GetSize();kn++) {
- ed=fSummEnergy[kn];
- if(ed>0) {
- fEventEnergy[ndep]=kn;
- if(ievent>1) {
- ed=ed/ievent;
- ed2=fSum2Energy[kn];
- ed2=ed2/ievent;
- ed2=100*TMath::Sqrt(TMath::Max(ed2-ed*ed,kzero))/ed;
- } else
- ed2=99;
- fSummEnergy[ndep]=ed;
- fSum2Energy[ndep]=TMath::Min((Float_t) 99.,TMath::Max(ed2,kzero));
- edtot+=ed;
- ndep++;
- }
- }
- for(kn=0;kn<(ndep-1)/3+1;kn++) {
- left=ndep-kn*3;
- for(i=0;i<(3<left?3:left);i++) {
- j=kn*3+i;
- id=Int_t (fEventEnergy[j]+0.1);
- printf(" %s %10.3f +- %10.3f%%;",gMC->VolName(id),fSummEnergy[j],fSum2Energy[j]);
- }
- printf("\n");
- }
- //
- // Relative energy loss in different detectors
- printf("******************** Relative Energy Loss per event ********************\n");
- printf("Total energy loss per event %10.3f GeV\n",edtot);
- for(kn=0;kn<(ndep-1)/5+1;kn++) {
- left=ndep-kn*5;
- for(i=0;i<(5<left?5:left);i++) {
- j=kn*5+i;
- id=Int_t (fEventEnergy[j]+0.1);
- printf(" %s %10.3f%%;",gMC->VolName(id),100*fSummEnergy[j]/edtot);
- }
- printf("\n");
- }
- for(kn=0;kn<75;kn++) printf("*");
- printf("\n");
- }
- //
- // Reset the TArray's
- // fEventEnergy.Set(0);
- // fSummEnergy.Set(0);
- // fSum2Energy.Set(0);
-}
-
-//_____________________________________________________________________________
-AliModule *AliRun::GetModule(const char *name) const
-{
- //
- // Return pointer to detector from name
- //
- return (AliModule*)fModules->FindObject(name);
-}
-
-//_____________________________________________________________________________
-AliDetector *AliRun::GetDetector(const char *name) const
-{
- //
- // Return pointer to detector from name
- //
- return (AliDetector*)fModules->FindObject(name);
-}
-
-//_____________________________________________________________________________
-Int_t AliRun::GetModuleID(const char *name) const
-{
- //
- // Return galice internal detector identifier from name
- //
- Int_t i=-1;
- TObject *mod=fModules->FindObject(name);
- if(mod) i=fModules->IndexOf(mod);
- return i;
-}
-
-//_____________________________________________________________________________
-Int_t AliRun::GetEvent(Int_t event)
-{
- //
- // Connect the Trees Kinematics and Hits for event # event
- // Set branch addresses
- //
-
- // Reset existing structures
- ResetHits();
- ResetDigits();
- ResetSDigits();
-
- // Delete Trees already connected
- if (fTreeH) delete fTreeH;
- if (fTreeD) delete fTreeD;
- if (fTreeR) delete fTreeR;
- if (fTreeS) delete fTreeS;
-
- // Create the particle stack
- if (fHeader) delete fHeader;
- fHeader = 0;
-
- // Get header from file
- if(fTreeE) {
- fTreeE->SetBranchAddress("Header", &fHeader);
- fTreeE->GetEntry(event);
- }
- else
- Error("GetEvent","Cannot find Header Tree (TE)\n");
-
- // Get the stack from the header
- if (fStack) delete fStack;
- fStack = fHeader->Stack();
- fStack->GetEvent(event);
- //
- TFile *file = fTreeE->GetCurrentFile();
- char treeName[20];
- file->cd();
-
- // Get Hits Tree header from file
- sprintf(treeName,"TreeH%d",event);
- fTreeH = (TTree*)gDirectory->Get(treeName);
- if (!fTreeH) {
- Error("GetEvent","cannot find Hits Tree for event:%d\n",event);
- }
-
- // Get Digits Tree header from file
- sprintf(treeName,"TreeD%d",event);
- fTreeD = (TTree*)gDirectory->Get(treeName);
- if (!fTreeD) {
- // Warning("GetEvent","cannot find Digits Tree for event:%d\n",event);
- }
-
- file->cd();
-
- // Get SDigits Tree header from file
- sprintf(treeName,"TreeS%d",event);
- fTreeS = (TTree*)gDirectory->Get(treeName);
- if (!fTreeS) {
- // Warning("GetEvent","cannot find SDigits Tree for event:%d\n",event);
- }
-
- file->cd();
-
- // Get Reconstruct Tree header from file
- sprintf(treeName,"TreeR%d",event);
- fTreeR = (TTree*)gDirectory->Get(treeName);
- if (!fTreeR) {
- // printf("WARNING: cannot find Reconstructed Tree for event:%d\n",event);
- }
-
- file->cd();
-
- // Set Trees branch addresses
- TIter next(fModules);
- AliModule *detector;
- while((detector = (AliModule*)next())) {
- detector->SetTreeAddress();
- }
-
- return fStack->GetNtrack();
-}
-
-//_____________________________________________________________________________
-TGeometry *AliRun::GetGeometry()
-{
- //
- // Import Alice geometry from current file
- // Return pointer to geometry object
- //
- if (!fGeometry) fGeometry = (TGeometry*)gDirectory->Get("AliceGeom");
- //
- // Unlink and relink nodes in detectors
- // This is bad and there must be a better way...
- //
-
- TIter next(fModules);
- AliModule *detector;
- while((detector = (AliModule*)next())) {
- TList *dnodes=detector->Nodes();
- Int_t j;
- TNode *node, *node1;
- for ( j=0; j<dnodes->GetSize(); j++) {
- node = (TNode*) dnodes->At(j);
- node1 = fGeometry->GetNode(node->GetName());
- dnodes->Remove(node);
- dnodes->AddAt(node1,j);
- }
- }
- return fGeometry;
-}
-
-//_____________________________________________________________________________
-void AliRun::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom,
- Float_t &e, Float_t *vpos, Float_t *polar,
- Float_t &tof)
-{
- // Delegate to stack
- //
- fStack->GetNextTrack(mtrack, ipart, pmom, e, vpos, polar, tof);
-}
-
-//_____________________________________________________________________________
-Int_t AliRun::GetPrimary(Int_t track) const
-{
- //
- // return number of primary that has generated track
- //
- return fStack->GetPrimary(track);
-}
-
-//_____________________________________________________________________________
-void AliRun::InitMC(const char *setup)
-{
- //
- // Initialize the Alice setup
- //
-
- if(fInitDone) {
- Warning("Init","Cannot initialise AliRun twice!\n");
- return;
- }
-
- gROOT->LoadMacro(setup);
- gInterpreter->ProcessLine(fConfigFunction.Data());
-
-
- gMC->DefineParticles(); //Create standard MC particles
-
- TObject *objfirst, *objlast;
-
- fNdets = fModules->GetLast()+1;
-
- //
- //=================Create Materials and geometry
- gMC->Init();
-
- // Added also after in case of interactive initialisation of modules
- fNdets = fModules->GetLast()+1;
-
- TIter next(fModules);
- AliModule *detector;
- while((detector = (AliModule*)next())) {
- detector->SetTreeAddress();
- objlast = gDirectory->GetList()->Last();
-
- // Add Detector histograms in Detector list of histograms
- if (objlast) objfirst = gDirectory->GetList()->After(objlast);
- else objfirst = gDirectory->GetList()->First();
- while (objfirst) {
- detector->Histograms()->Add(objfirst);
- objfirst = gDirectory->GetList()->After(objfirst);
- }
- }
- ReadTransPar(); //Read the cuts for all materials
-
- MediaTable(); //Build the special IMEDIA table
-
- //Initialise geometry deposition table
- fEventEnergy.Set(gMC->NofVolumes()+1);
- fSummEnergy.Set(gMC->NofVolumes()+1);
- fSum2Energy.Set(gMC->NofVolumes()+1);
-
- //Compute cross-sections
- gMC->BuildPhysics();
-
- //Write Geometry object to current file.
- fGeometry->Write();
-
- fInitDone = kTRUE;
-
- fMCQA = new AliMCQA(fNdets);
-
- AliConfig *config = AliConfig::Instance();
- //
- // Save stuff at the beginning of the file to avoid file corruption
- Write();
-}
-
-//_____________________________________________________________________________
-void AliRun::MediaTable()
-{
- //
- // Built media table to get from the media number to
- // the detector id
- //
- Int_t kz, nz, idt, lz, i, k, ind;
- // Int_t ibeg;
- TObjArray &dets = *gAlice->Detectors();
- AliModule *det;
- //
- // For all detectors
- for (kz=0;kz<fNdets;kz++) {
- // If detector is defined
- if((det=(AliModule*) dets[kz])) {
- TArrayI &idtmed = *(det->GetIdtmed());
- for(nz=0;nz<100;nz++) {
- // Find max and min material number
- if((idt=idtmed[nz])) {
- det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
- det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
- }
- }
- if(det->LoMedium() > det->HiMedium()) {
- det->LoMedium() = 0;
- det->HiMedium() = 0;
- } else {
- if(det->HiMedium() > fImedia->GetSize()) {
- Error("MediaTable","Increase fImedia from %d to %d",
- fImedia->GetSize(),det->HiMedium());
- return;
- }
- // Tag all materials in rage as belonging to detector kz
- for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
- (*fImedia)[lz]=kz;
- }
- }
- }
- }
- //
- // Print summary table
- printf(" Traking media ranges:\n");
- for(i=0;i<(fNdets-1)/6+1;i++) {
- for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
- ind=i*6+k;
- det=(AliModule*)dets[ind];
- if(det)
- printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
- det->HiMedium());
- else
- printf(" %6s: %3d -> %3d;","NULL",0,0);
- }
- printf("\n");
- }
-}
-
-//____________________________________________________________________________
-void AliRun::SetGenerator(AliGenerator *generator)
-{
- //
- // Load the event generator
- //
- if(!fGenerator) fGenerator = generator;
-}
-
-//____________________________________________________________________________
-void AliRun::ResetGenerator(AliGenerator *generator)
-{
- //
- // Load the event generator
- //
- if(fGenerator)
- if(generator)
- Warning("ResetGenerator","Replacing generator %s with %s\n",
- fGenerator->GetName(),generator->GetName());
- else
- Warning("ResetGenerator","Replacing generator %s with NULL\n",
- fGenerator->GetName());
- fGenerator = generator;