X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliMC.cxx;h=0cca364f14b10177f629b5bb3dc9db985b99752f;hb=6ef9caeb5bb373151828235a601b5ab77f17e1b7;hp=13380f34d50cb181aeb07668110808167942b62e;hpb=504b172db208d0582ee848cf4d9fd438bfe64849;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliMC.cxx b/STEER/AliMC.cxx index 13380f34d50..0cca364f14b 100644 --- a/STEER/AliMC.cxx +++ b/STEER/AliMC.cxx @@ -15,22 +15,41 @@ /* $Id$ */ -#include +// This class is extracted from the AliRun class +// and contains all the MC-related functionality +// The number of dependencies has to be reduced... +// Author: F.Carminati +// Federico.Carminati@cern.ch + +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include +#include "AliCDBEntry.h" +#include "AliCDBManager.h" +#include "AliCDBStorage.h" #include "AliDetector.h" #include "AliGenerator.h" +#include "AliGeomManager.h" #include "AliHeader.h" +#include "AliHit.h" #include "AliLego.h" +#include "AliLog.h" #include "AliMC.h" -#include "AliMCQA.h" +#include "AliMagF.h" #include "AliRun.h" +#include "AliSimulation.h" #include "AliStack.h" #include "AliTrackReference.h" - ClassImp(AliMC) //_______________________________________________________________________ @@ -41,13 +60,20 @@ AliMC::AliMC() : fSum2Energy(0), fTrRmax(1.e10), fTrZmax(1.e10), + fRDecayMax(1.e10), + fRDecayMin(-1.), + fDecayPdg(0), fImedia(0), fTransParName("\0"), - fMCQA(0), fHitLists(0), - fTrackReferences(0) + fTmpTreeTR(0), + fTmpFileTR(0), + fTrackReferences(), + fTmpTrackReferences() { + //default constructor + DecayLimits(); } //_______________________________________________________________________ @@ -59,109 +85,157 @@ AliMC::AliMC(const char *name, const char *title) : fSum2Energy(0), fTrRmax(1.e10), fTrZmax(1.e10), + fRDecayMax(1.e10), + fRDecayMin(-1.), + fDecayPdg(0), fImedia(new TArrayI(1000)), fTransParName("\0"), - fMCQA(0), fHitLists(new TList()), - fTrackReferences(new TClonesArray("AliTrackReference", 100)) + fTmpTreeTR(0), + fTmpFileTR(0), + fTrackReferences("AliTrackReference", 100), + fTmpTrackReferences("AliTrackReference", 100) { - + //constructor // Set transport parameters SetTransPar(); - + DecayLimits(); // Prepare the tracking medium lists for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99; - -} - -//_______________________________________________________________________ -AliMC::AliMC(const AliMC &mc) : - TVirtualMCApplication(mc), - fGenerator(0), - fEventEnergy(0), - fSummEnergy(0), - fSum2Energy(0), - fTrRmax(1.e10), - fTrZmax(1.e10), - fImedia(0), - fTransParName("\0"), - fMCQA(0), - fHitLists(0), - fTrackReferences(0) -{ - // - // Copy constructor for AliRun - // - mc.Copy(*this); } //_______________________________________________________________________ AliMC::~AliMC() { + //destructor delete fGenerator; delete fImedia; - delete fMCQA; delete fHitLists; // Delete track references - if (fTrackReferences) { - fTrackReferences->Delete(); - delete fTrackReferences; - fTrackReferences = 0; - } - -} - -//_______________________________________________________________________ -void AliMC::Copy(AliMC &) const -{ - Fatal("Copy","Not implemented!\n"); } //_______________________________________________________________________ void AliMC::ConstructGeometry() { // - // Create modules, materials, geometry - // - + // Either load geometry from file or create it through usual + // loop on detectors. In the first case the method + // AliModule::CreateMaterials() only builds fIdtmed and is postponed + // at InitGeometry(). + // + + if(gAlice->IsRootGeometry()){ //load geometry either from CDB or from file + if(gAlice->IsGeomFromCDB()){ + AliInfo("Loading geometry from CDB default storage"); + AliCDBPath path("GRP","Geometry","Data"); + AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath()); + if(!entry) AliFatal("Unable to load geometry from CDB!"); + entry->SetOwner(0); + gGeoManager = (TGeoManager*) entry->GetObject(); + if (!gGeoManager) AliFatal("TGeoManager object not found in the specified CDB entry!"); + }else{ + // Load geometry + const char *geomfilename = gAlice->GetGeometryFileName(); + if(gSystem->ExpandPathName(geomfilename)){ + AliInfo(Form("Loading geometry from file:\n %40s",geomfilename)); + TGeoManager::Import(geomfilename); + }else{ + AliInfo(Form("Geometry file %40s not found!\n",geomfilename)); + return; + } + } + }else{ + // Create modules, materials, geometry + if (!gGeoManager) new TGeoManager("ALICE", "ALICE geometry"); TStopwatch stw; TIter next(gAlice->Modules()); AliModule *detector; - if (GetDebug()) Info("ConstructGeometry","Geometry creation:"); + AliDebug(1, "Geometry creation:"); while((detector = dynamic_cast(next()))) { stw.Start(); // Initialise detector materials and geometry detector->CreateMaterials(); detector->CreateGeometry(); - printf("%10s R:%.2fs C:%.2fs\n", - detector->GetName(),stw.RealTime(),stw.CpuTime()); + AliInfo(Form("%10s R:%.2fs C:%.2fs", + detector->GetName(),stw.RealTime(),stw.CpuTime())); } + } + +} + +//_______________________________________________________________________ +Bool_t AliMC::MisalignGeometry() +{ +// Call misalignment code if AliSimulation object was defined. + + if(!gAlice->IsRootGeometry()){ + //Set alignable volumes for the whole geometry + SetAllAlignableVolumes(); + } + // Misalign geometry via AliSimulation instance + if (!AliSimulation::Instance()) return kFALSE; + AliGeomManager::SetGeometry(gGeoManager); + if(!AliGeomManager::CheckSymNamesLUT("ALL")) + AliFatal("Current loaded geometry differs in the definition of symbolic names!"); + + return AliSimulation::Instance()->MisalignGeometry(AliRunLoader::Instance()); +} + +//_______________________________________________________________________ +void AliMC::ConstructOpGeometry() +{ + // + // Loop all detector modules and call DefineOpticalProperties() method + // + + TIter next(gAlice->Modules()); + AliModule *detector; + AliInfo("Optical properties definition"); + while((detector = dynamic_cast(next()))) { + // Initialise detector geometry + if(gAlice->IsRootGeometry()) detector->CreateMaterials(); + // Initialise detector optical properties + detector->DefineOpticalProperties(); + } } //_______________________________________________________________________ void AliMC::InitGeometry() { // - // Initialize detectors and display geometry + // Initialize detectors // - printf("Initialisation:\n"); - TStopwatch stw; - TIter next(gAlice->Modules()); - AliModule *detector; - while((detector = dynamic_cast(next()))) { - stw.Start(); - // Initialise detector and display geometry - detector->Init(); - detector->BuildGeometry(); - printf("%10s R:%.2fs C:%.2fs\n", - detector->GetName(),stw.RealTime(),stw.CpuTime()); - } - + AliInfo("Initialisation:"); + TStopwatch stw; + TIter next(gAlice->Modules()); + AliModule *detector; + while((detector = dynamic_cast(next()))) { + stw.Start(); + detector->Init(); + AliInfo(Form("%10s R:%.2fs C:%.2fs", + detector->GetName(),stw.RealTime(),stw.CpuTime())); + } +} + +//_______________________________________________________________________ +void AliMC::SetAllAlignableVolumes() +{ + // + // Add alignable volumes (TGeoPNEntries) looping on all + // active modules + // + + AliInfo(Form("Setting entries for all alignable volumes of active detectors")); + AliModule *detector; + TIter next(gAlice->Modules()); + while((detector = dynamic_cast(next()))) { + detector->AddAlignableVolumes(); + } } //_______________________________________________________________________ -void AliMC::GeneratePrimaries() +void AliMC::GeneratePrimaries() { // // Generate primary particles and fill them in the stack. @@ -185,13 +259,16 @@ void AliMC::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()); + if(fGenerator) { + if(generator) { + AliWarning(Form("Replacing generator %s with %s", + fGenerator->GetName(),generator->GetName())); + } + else { + AliWarning(Form("Replacing generator %s with NULL", + fGenerator->GetName())); + } + } fGenerator = generator; } @@ -199,13 +276,12 @@ void AliMC::ResetGenerator(AliGenerator *generator) void AliMC::FinishRun() { // Clean generator information - if (GetDebug()) Info("FinishRun"," fGenerator->FinishRun()"); + AliDebug(1, "fGenerator->FinishRun()"); fGenerator->FinishRun(); //Output energy summary tables - if (GetDebug()) Info("FinishRun"," EnergySummary()"); - EnergySummary(); - + AliDebug(1, "EnergySummary()"); + ToAliDebug(1, EnergySummary()); } //_______________________________________________________________________ @@ -218,20 +294,19 @@ void AliMC::BeginPrimary() // Reset Hits info ResetHits(); ResetTrackReferences(); - } //_______________________________________________________________________ void AliMC::PreTrack() { + // Actions before the track's transport + TObjArray &dets = *gAlice->Modules(); AliModule *module; for(Int_t i=0; i<=gAlice->GetNdets(); i++) if((module = dynamic_cast(dets[i]))) module->PreTrack(); - - fMCQA->PreTrack(); } //_______________________________________________________________________ @@ -240,23 +315,42 @@ void AliMC::Stepping() // // Called at every step during transport // - - Int_t id = DetFromMate(gMC->GetMedium()); + Int_t id = DetFromMate(gMC->CurrentMedium()); if (id < 0) return; + + if ( gMC->IsNewTrack() && + gMC->TrackTime() == 0. && + fRDecayMin >= 0. && + fRDecayMax > fRDecayMin && + gMC->TrackPid() == fDecayPdg ) + { + FixParticleDecaytime(); + } + + + // // --- If lego option, do it and leave - if (gAlice->Lego()) - gAlice->Lego()->StepManager(); + if (AliSimulation::Instance()->Lego()) + AliSimulation::Instance()->Lego()->StepManager(); else { Int_t copy; //Update energy deposition tables AddEnergyDeposit(gMC->CurrentVolID(copy),gMC->Edep()); - + // + // write tracke reference for track which is dissapearing - MI + + if (gMC->IsTrackDisappeared() && !(gMC->IsTrackAlive())) { + if (gMC->Etot() > 0.05) AddTrackReference(GetCurrentTrackNumber(), + AliTrackReference::kDisappeared); + + + } + //Call the appropriate stepping routine; AliModule *det = dynamic_cast(gAlice->Modules()->At(id)); if(det && det->StepManagerIsEnabled()) { - fMCQA->StepManager(id); det->StepManager(); } } @@ -265,7 +359,7 @@ void AliMC::Stepping() //_______________________________________________________________________ void AliMC::EnergySummary() { - // + //e // Print summary of deposited energy // @@ -274,7 +368,7 @@ void AliMC::EnergySummary() Float_t ed, ed2; Int_t kn, i, left, j, id; const Float_t kzero=0; - Int_t ievent=gAlice->GetRunLoader()->GetHeader()->GetEvent()+1; + Int_t ievent=AliRunLoader::Instance()->GetHeader()->GetEvent()+1; // // Energy loss information if(ievent) { @@ -334,18 +428,15 @@ void AliMC::BeginEvent() // // Clean-up previous event // Energy scores - if (GetDebug()) - { - Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); - Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); - Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); - Info("BeginEvent"," BEGINNING EVENT "); - Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); - Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); - Info("BeginEvent",">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); - } + AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); + AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); + AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); + AliDebug(1, " BEGINNING EVENT "); + AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); + AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); + AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); - AliRunLoader *runloader=gAlice->GetRunLoader(); + AliRunLoader *runloader=AliRunLoader::Instance(); /*******************************/ /* Clean after eventual */ @@ -356,24 +447,25 @@ void AliMC::BeginEvent() //Set the next event in Run Loader -> Cleans trees (TreeK and all trees in detectors), gAlice->SetEventNrInRun(gAlice->GetEventNrInRun()+1); runloader->SetEventNumber(gAlice->GetEventNrInRun());// sets new files, cleans the previous event stuff, if necessary, etc., - if (GetDebug()) Info("BeginEvent","EventNr is %d",gAlice->GetEventNrInRun()); + AliDebug(1, Form("EventNr is %d",gAlice->GetEventNrInRun())); fEventEnergy.Reset(); // Clean detector information if (runloader->Stack()) - runloader->Stack()->Reset();//clean stack -> tree is unloaded + runloader->Stack()->Reset();//clean stack -> tree is unloaded else - runloader->MakeStack();//or make a new one + runloader->MakeStack();//or make a new one - if(gAlice->Lego() == 0x0) - { - if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTree(K)"); - runloader->MakeTree("K"); - } - - if (GetDebug()) Info("BeginEvent"," gMC->SetStack(fRunLoader->Stack())"); - gMC->SetStack(gAlice->GetRunLoader()->Stack());//Was in InitMC - but was moved here + + if(AliSimulation::Instance()->Lego() == 0x0) + { + AliDebug(1, "fRunLoader->MakeTree(K)"); + runloader->MakeTree("K"); + } + + AliDebug(1, "gMC->SetStack(fRunLoader->Stack())"); + gMC->SetStack(runloader->Stack());//Was in InitMC - but was moved here //because we don't have guarantee that //stack pointer is not going to change from event to event //since it bellobgs to header and is obtained via RunLoader @@ -381,38 +473,33 @@ void AliMC::BeginEvent() // Reset all Detectors & kinematics & make/reset trees // - runloader->GetHeader()->Reset(gAlice->GetRunNumber(),gAlice->GetEvNumber(), - gAlice->GetEventNrInRun()); + runloader->GetHeader()->Reset(AliCDBManager::Instance()->GetRun(),gAlice->GetEvNumber(), + gAlice->GetEventNrInRun()); // fRunLoader->WriteKinematics("OVERWRITE"); is there any reason to rewrite here since MakeTree does so - if(gAlice->Lego()) - { - gAlice->Lego()->BeginEvent(); - return; - } + if(AliSimulation::Instance()->Lego()) + { + AliSimulation::Instance()->Lego()->BeginEvent(); + return; + } - if (GetDebug()) Info("BeginEvent"," ResetHits()"); + AliDebug(1, "ResetHits()"); ResetHits(); - if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTree(H)"); + AliDebug(1, "fRunLoader->MakeTree(H)"); runloader->MakeTree("H"); - if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTrackRefsContainer()"); - runloader->MakeTrackRefsContainer();//for insurance + MakeTmpTrackRefsTree(); //create new branches and SetAdresses TIter next(gAlice->Modules()); AliModule *detector; while((detector = (AliModule*)next())) { - if (GetDebug()) Info("BeginEvent"," %s->MakeBranch(H)",detector->GetName()); - detector->MakeBranch("H"); - if (GetDebug()) Info("BeginEvent"," %s->MakeBranchTR()",detector->GetName()); - detector->MakeBranchTR(); - if (GetDebug()) Info("BeginEvent"," %s->SetTreeAddress()",detector->GetName()); - detector->SetTreeAddress(); + AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName())); + detector->MakeBranch("H"); } } @@ -429,15 +516,43 @@ void AliMC::ResetHits() } } +//_______________________________________________________________________ +void AliMC::ResetDigits() +{ + // + // Reset all Detectors digits + // + TIter next(gAlice->Modules()); + AliModule *detector; + while((detector = dynamic_cast(next()))) { + detector->ResetDigits(); + } +} + +//_______________________________________________________________________ +void AliMC::ResetSDigits() +{ + // + // Reset all Detectors digits + // + TIter next(gAlice->Modules()); + AliModule *detector; + while((detector = dynamic_cast(next()))) { + detector->ResetSDigits(); + } +} + //_______________________________________________________________________ void AliMC::PostTrack() { - TObjArray &dets = *gAlice->Modules(); - AliModule *module; + // Posts tracks for each module - for(Int_t i=0; i<=gAlice->GetNdets(); i++) - if((module = dynamic_cast(dets[i]))) - module->PostTrack(); + TObjArray &dets = *gAlice->Modules(); + AliModule *module; + + for(Int_t i=0; i<=gAlice->GetNdets(); i++) + if((module = dynamic_cast(dets[i]))) + module->PostTrack(); } //_______________________________________________________________________ @@ -446,11 +561,17 @@ void AliMC::FinishPrimary() // // Called at the end of each primary track // - AliRunLoader *runloader=gAlice->GetRunLoader(); + + AliRunLoader *runloader=AliRunLoader::Instance(); // static Int_t count=0; // const Int_t times=10; // This primary is finished, purify stack - runloader->Stack()->PurifyKine(); +#if ROOT_VERSION_CODE > 262152 + if (!(gMC->SecondariesAreOrdered())) { + if (runloader->Stack()->ReorderKine()) RemapHits(); + } +#endif + if (runloader->Stack()->PurifyKine()) RemapHits(); TIter next(gAlice->Modules()); AliModule *detector; @@ -465,7 +586,42 @@ void AliMC::FinishPrimary() } // Write out track references if any - if (runloader->TreeTR()) runloader->TreeTR()->Fill(); + if (fTmpTreeTR) fTmpTreeTR->Fill(); +} + +void AliMC::RemapHits() +{ +// +// Remaps the track labels of the hits + AliRunLoader *runloader=AliRunLoader::Instance(); + AliStack* stack = runloader->Stack(); + TList* hitLists = GetHitLists(); + TIter next(hitLists); + TCollection *hitList; + + while((hitList = dynamic_cast(next()))) { + TIter nexthit(hitList); + AliHit *hit; + while((hit = dynamic_cast(nexthit()))) { + hit->SetTrack(stack->TrackLabel(hit->GetTrack())); + } + } + + // + // This for detectors which have a special mapping mechanism + // for hits, such as TPC and TRD + // + + + TObjArray* modules = gAlice->Modules(); + TIter nextmod(modules); + AliModule *module; + while((module = (AliModule*) nextmod())) { + AliDetector* det = dynamic_cast (module); + if (det) det->RemapTrackHitIDs(stack->TrackLabelMap()); + } + // + RemapTrackReferencesIDs(stack->TrackLabelMap()); } //_______________________________________________________________________ @@ -474,9 +630,8 @@ void AliMC::FinishEvent() // // Called at the end of the event. // - - // - if(gAlice->Lego()) gAlice->Lego()->FinishEvent(); + + if(AliSimulation::Instance()->Lego()) AliSimulation::Instance()->Lego()->FinishEvent(); TIter next(gAlice->Modules()); AliModule *detector; @@ -492,22 +647,24 @@ void AliMC::FinishEvent() fSum2Energy[i]+=fEventEnergy[i]*fEventEnergy[i]; } - AliRunLoader *runloader=gAlice->GetRunLoader(); + AliRunLoader *runloader=AliRunLoader::Instance(); AliHeader* header = runloader->GetHeader(); AliStack* stack = runloader->Stack(); if ( (header == 0x0) || (stack == 0x0) ) {//check if we got header and stack. If not cry and exit aliroot - Fatal("AliRun","Can not get the stack or header from LOADER"); + AliFatal("Can not get the stack or header from LOADER"); return;//never reached } // Update Header information header->SetNprimary(stack->GetNprimary()); header->SetNtrack(stack->GetNtrack()); - // Write out the kinematics - stack->FinishEvent(); + if (!AliSimulation::Instance()->Lego()) stack->FinishEvent(); + + // Synchronize the TreeTR with TreeK + if (fTmpTreeTR) ReorderAndExpandTreeTR(); // Write out the event Header information TTree* treeE = runloader->TreeE(); @@ -518,42 +675,36 @@ void AliMC::FinishEvent() } else { - Error("FinishEvent","Can not get TreeE from RL"); + AliError("Can not get TreeE from RL"); } - if(gAlice->Lego() == 0x0) + if(AliSimulation::Instance()->Lego() == 0x0) { runloader->WriteKinematics("OVERWRITE"); runloader->WriteTrackRefs("OVERWRITE"); runloader->WriteHits("OVERWRITE"); } - if (GetDebug()) - { - Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); - Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); - Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); - Info("FinishEvent"," FINISHING EVENT "); - Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); - Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); - Info("FinishEvent","<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); - } + AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); + AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); + AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); + AliDebug(1, " FINISHING EVENT "); + AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); + AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); + AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); } -//_______________________________________________________________________ -void AliMC::Field(const Double_t* x, Double_t* b) const -{ - gAlice->Field(x,b); -} - //_______________________________________________________________________ void AliMC::Init() { + // MC initialization //=================Create Materials and geometry gMC->Init(); - gMC->DefineParticles(); //Create standard MC particles - + // Set alignable volumes for the whole geometry (with old root) +#if ROOT_VERSION_CODE < 331527 + SetAllAlignableVolumes(); +#endif //Read the cuts for all materials ReadTransPar(); //Build the special IMEDIA table @@ -562,17 +713,11 @@ void AliMC::Init() //Compute cross-sections gMC->BuildPhysics(); - //Write Geometry object to current file. - gAlice->GetRunLoader()->WriteGeometry(); - //Initialise geometry deposition table fEventEnergy.Set(gMC->NofVolumes()+1); fSummEnergy.Set(gMC->NofVolumes()+1); fSum2Energy.Set(gMC->NofVolumes()+1); - // - fMCQA = new AliMCQA(gAlice->GetNdets()); - // Register MC in configuration AliConfig::Instance()->Add(gMC); @@ -598,6 +743,7 @@ void AliMC::MediaTable() if((det=dynamic_cast(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; @@ -609,8 +755,8 @@ void AliMC::MediaTable() det->HiMedium() = 0; } else { if(det->HiMedium() > fImedia->GetSize()) { - Error("MediaTable","Increase fImedia from %d to %d", - fImedia->GetSize(),det->HiMedium()); + AliError(Form("Increase fImedia from %d to %d", + fImedia->GetSize(),det->HiMedium())); return; } // Tag all materials in rage as belonging to detector kz @@ -622,7 +768,8 @@ void AliMC::MediaTable() } // // Print summary table - printf(" Traking media ranges:\n"); + AliInfo("Tracking media ranges:"); + ToAliInfo( for(i=0;i<(ndets-1)/6+1;i++) { for(k=0;k< (6GetName()); + AliWarning(Form("Invalid tracking medium code %d for %s",itmed,mod->GetName())); continue; } // Set energy thresholds for(kz=0;kz=0) { - if(fDebug) printf(" * %-6s set to %10.3E for tracking medium code %4d for %s\n", - kpars[kz],cut[kz],itmed,mod->GetName()); + AliDebug(2, Form("%-6s set to %10.3E for tracking medium code %4d for %s", + kpars[kz],cut[kz],itmed,mod->GetName())); gMC->Gstpar(ktmed,kpars[kz],cut[kz]); } } // Set transport mechanisms for(kz=0;kz=0) { - if(fDebug) printf(" * %-6s set to %10d for tracking medium code %4d for %s\n", - kpars[kncuts+kz],flag[kz],itmed,mod->GetName()); + AliDebug(2, Form("%-6s set to %10d for tracking medium code %4d for %s", + kpars[kncuts+kz],flag[kz],itmed,mod->GetName())); gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz])); } } } else { - Warning("ReadTransPar","Invalid medium code %d *\n",itmed); + AliWarning(Form("Invalid medium code %d",itmed)); continue; } } else { - if(fDebug) printf("%s::ReadTransParModule: %s not present\n",ClassName(),detName); + AliDebug(1, Form("%s not present",detName)); continue; } } @@ -756,19 +893,6 @@ void AliMC::SetTransPar(const char *filename) fTransParName = filename; } -//_______________________________________________________________________ -void AliMC::Browse(TBrowser *b) -{ - // - // Called when the item "Run" is clicked on the left pane - // of the Root browser. - // It displays the Root Trees and all detectors. - // - //detectors are in folders anyway - b->Add(fMCQA,"AliMCQA"); -} - -//PH //_______________________________________________________________________ void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const { @@ -794,7 +918,7 @@ Int_t AliMC::GetCurrentTrackNumber() const { // // Returns current track // - return gAlice->GetRunLoader()->Stack()->GetCurrentTrackNumber(); + return AliRunLoader::Instance()->Stack()->GetCurrentTrackNumber(); } //_______________________________________________________________________ @@ -803,7 +927,7 @@ void AliMC::DumpPart (Int_t i) const // // Dumps particle i in the stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader->Stack()) runloader->Stack()->DumpPart(i); } @@ -814,7 +938,7 @@ void AliMC::DumpPStack () const // // Dumps the particle stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader->Stack()) runloader->Stack()->DumpPStack(); } @@ -825,7 +949,7 @@ Int_t AliMC::GetNtrack() const { // Returns number of tracks in stack // Int_t ntracks = -1; - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader->Stack()) ntracks = runloader->Stack()->GetNtrack(); return ntracks; @@ -838,7 +962,7 @@ Int_t AliMC::GetPrimary(Int_t track) const // return number of primary that has generated track // Int_t nprimary = -999; - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader->Stack()) nprimary = runloader->Stack()->GetPrimary(track); return nprimary; @@ -847,7 +971,9 @@ Int_t AliMC::GetPrimary(Int_t track) const //_______________________________________________________________________ TParticle* AliMC::Particle(Int_t i) const { - AliRunLoader * runloader = gAlice->GetRunLoader(); + // Returns the i-th particle from the stack taking into account + // the remaping done by PurifyKine + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) return runloader->Stack()->Particle(i); @@ -855,11 +981,11 @@ TParticle* AliMC::Particle(Int_t i) const } //_______________________________________________________________________ -TObjArray* AliMC::Particles() const { +const TObjArray* AliMC::Particles() const { // // Returns pointer to Particles array // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) return runloader->Stack()->Particles(); @@ -867,13 +993,13 @@ TObjArray* AliMC::Particles() const { } //_______________________________________________________________________ -void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, Float_t *pmom, - Float_t *vpos, Float_t *polar, Float_t tof, - TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) +void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, const Float_t *pmom, + const Float_t *vpos, const Float_t *polar, Float_t tof, + TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const { // Delegate to stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->PushTrack(done, parent, pdg, pmom, vpos, polar, tof, @@ -885,11 +1011,11 @@ void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, Double_t px, Double_t py, Double_t pz, Double_t e, Double_t vx, Double_t vy, Double_t vz, Double_t tof, Double_t polx, Double_t poly, Double_t polz, - TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) + TMCProcess mech, Int_t &ntr, Float_t weight, Int_t is) const { // Delegate to stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->PushTrack(done, parent, pdg, px, py, pz, e, vx, vy, vz, tof, @@ -897,62 +1023,61 @@ void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, } //_______________________________________________________________________ -void AliMC::SetHighWaterMark(const Int_t nt) +void AliMC::SetHighWaterMark(Int_t nt) const { // // Set high water mark for last track in event - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->SetHighWaterMark(nt); } //_______________________________________________________________________ -void AliMC::KeepTrack(const Int_t track) +void AliMC::KeepTrack(Int_t track) const { // // Delegate to stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->KeepTrack(track); } //_______________________________________________________________________ -void AliMC::FlagTrack(Int_t track) +void AliMC::FlagTrack(Int_t track) const { // Delegate to stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->FlagTrack(track); } //_______________________________________________________________________ -void AliMC::SetCurrentTrack(Int_t track) +void AliMC::SetCurrentTrack(Int_t track) const { // // Set current track number // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) runloader->Stack()->SetCurrentTrack(track); } //_______________________________________________________________________ -void AliMC::AddTrackReference(Int_t label){ +AliTrackReference* AliMC::AddTrackReference(Int_t label, Int_t id) +{ // // add a trackrefernce to the list - if (!fTrackReferences) { - cerr<<"Container trackrefernce not active\n"; - return; - } - Int_t nref = fTrackReferences->GetEntriesFast(); - TClonesArray &lref = *fTrackReferences; - new(lref[nref]) AliTrackReference(label); + Int_t primary = GetPrimary(label); + Particle(primary)->SetBit(kKeepBit); + + Int_t nref = fTmpTrackReferences.GetEntriesFast(); + return new(fTmpTrackReferences[nref]) AliTrackReference(label, id); } @@ -963,33 +1088,206 @@ void AliMC::ResetTrackReferences() // // Reset all references // - if (fTrackReferences) fTrackReferences->Clear(); - - TIter next(gAlice->Modules()); - AliModule *detector; - while((detector = dynamic_cast(next()))) { - detector->ResetTrackReferences(); - } + fTmpTrackReferences.Clear(); } +//_______________________________________________________________________ void AliMC::RemapTrackReferencesIDs(Int_t *map) { // // Remapping track reference // Called at finish primary // - if (!fTrackReferences) return; - for (Int_t i=0;iGetEntries();i++){ - AliTrackReference * ref = dynamic_cast(fTrackReferences->UncheckedAt(i)); - if (ref) { - Int_t newID = map[ref->GetTrack()]; - if (newID>=0) ref->SetTrack(newID); - else { - //ref->SetTrack(-1); - ref->SetBit(kNotDeleted,kFALSE); - fTrackReferences->RemoveAt(i); - } - } + + Int_t nEntries = fTmpTrackReferences.GetEntries(); + for (Int_t i=0; i < nEntries; i++){ + AliTrackReference * ref = dynamic_cast(fTmpTrackReferences.UncheckedAt(i)); + if (ref) { + Int_t newID = map[ref->GetTrack()]; + if (newID>=0) ref->SetTrack(newID); + else { + ref->SetBit(kNotDeleted,kFALSE); + fTmpTrackReferences.RemoveAt(i); + } + } // if ref } - fTrackReferences->Compress(); + fTmpTrackReferences.Compress(); +} + +//_______________________________________________________________________ +void AliMC::FixParticleDecaytime() +{ + // + // Fix the particle decay time according to rmin and rmax for decays + // + + TLorentzVector p; + gMC->TrackMomentum(p); + Double_t tmin, tmax; + Double_t b; + + // Transverse velocity + Double_t vt = p.Pt() / p.E(); + + if ((b = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField()) > 0.) { // [kG] + + // Radius of helix + + Double_t rho = p.Pt() / 0.0003 / b; // [cm] + + // Revolution frequency + + Double_t omega = vt / rho; + + // Maximum and minimum decay time + // + // Check for curlers first + const Double_t kOvRhoSqr2 = 1./(rho*TMath::Sqrt(2.)); + if (fRDecayMax * kOvRhoSqr2 > 1.) return; + + // + + tmax = TMath::ACos((1.-fRDecayMax*kOvRhoSqr2)*(1.+fRDecayMax*kOvRhoSqr2)) / omega; // [ct] + tmin = TMath::ACos((1.-fRDecayMin*kOvRhoSqr2)*(1.+fRDecayMin*kOvRhoSqr2)) / omega; // [ct] + } else { + tmax = fRDecayMax / vt; // [ct] + tmin = fRDecayMin / vt; // [ct] + } + + // + // Dial t using the two limits + Double_t t = tmin + (tmax - tmin) * gRandom->Rndm(); // [ct] + // + // + // Force decay time in transport code + // + gMC->ForceDecayTime(t / 2.99792458e10); +} + +void AliMC::MakeTmpTrackRefsTree() +{ + // Make the temporary track reference tree + fTmpFileTR = new TFile("TrackRefsTmp.root", "recreate"); + fTmpTreeTR = new TTree("TreeTR", "Track References"); + TClonesArray* pRef = &fTmpTrackReferences; + fTmpTreeTR->Branch("TrackReferences", "TClonesArray", &pRef, 4000); +} + +//_______________________________________________________________________ +void AliMC::ReorderAndExpandTreeTR() +{ +// +// Reorder and expand the temporary track reference tree in order to match the kinematics tree +// + + AliRunLoader *rl = AliRunLoader::Instance(); +// +// TreeTR + AliDebug(1, "fRunLoader->MakeTrackRefsContainer()"); + rl->MakeTrackRefsContainer(); + TTree * treeTR = rl->TreeTR(); + if (treeTR){ + // make branch for central track references + TBranch *branch; + TClonesArray* pRef = &fTrackReferences; + branch = treeTR->Branch("TrackReferences", &pRef); + branch->SetAddress(&pRef); + } + + AliStack* stack = rl->Stack(); + Int_t np = stack->GetNprimary(); + Int_t nt = fTmpTreeTR->GetEntries(); + // + // Loop over tracks and find the secondaries with the help of the kine tree + Int_t ifills = 0; + Int_t it = 0; + for (Int_t ip = np - 1; ip > -1; ip--) { + TParticle *part = stack->Particle(ip); + //printf("Particle %5d %5d %5d %5d %5d \n", ip, part->GetPdgCode(), part->GetFirstMother(), part->GetFirstDaughter(), part->GetLastDaughter()); + + // Skip primaries that have not been transported + Int_t dau1 = part->GetFirstDaughter(); + Int_t dau2 = -1; + if (!part->TestBit(kTransportBit)) continue; + // + fTmpTreeTR->GetEntry(it++); + Int_t nh = fTmpTrackReferences.GetEntries(); + // Determine range of secondaries produced by this primary + if (dau1 > -1) { + Int_t inext = ip - 1; + while (dau2 < 0) { + if (inext >= 0) { + part = stack->Particle(inext); + dau2 = part->GetFirstDaughter(); + if (!(part->TestBit(kTransportBit)) || dau2 == -1 || dau2 < np) { +// if (dau2 == -1 || dau2 < np) { + dau2 = -1; + } else { + dau2--; + } + } else { + dau2 = stack->GetNtrack() - 1; + } + inext--; + } // find upper bound + } // dau2 < 0 +// printf("Check (1) %5d %5d %5d %5d %5d \n", ip, np, it, dau1, dau2); + // + // Loop over reference hits and find secondary label + for (Int_t id = dau1; (id <= dau2) && (dau1 > -1); id++) { + for (Int_t ih = 0; ih < nh; ih++) { + AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences.At(ih); + Int_t label = tr->Label(); + // Skip primaries + if (label == ip) continue; + if (label > dau2 || label < dau1) + AliWarning(Form("Track Reference Label out of range !: %5d %5d %5d \n", label, dau1, dau2)); + if (label == id) { + // secondary found + Int_t nref = fTrackReferences.GetEntriesFast(); + new(fTrackReferences[nref]) AliTrackReference(*tr); + } + } // hits + treeTR->Fill(); + fTrackReferences.Clear(); + ifills++; + } // daughters + } // tracks + // + // Now loop again and write the primaries + it = nt - 1; + for (Int_t ip = 0; ip < np; ip++) { + TParticle* part = stack->Particle(ip); +// if ((part->GetFirstDaughter() == -1 && part->GetStatusCode() <= 1) || part->GetFirstDaughter() >= np) + if (part->TestBit(kTransportBit)) + { + // Skip particles that have not been transported + fTmpTreeTR->GetEntry(it--); + Int_t nh = fTmpTrackReferences.GetEntries(); + // + // Loop over reference hits and find primary labels + for (Int_t ih = 0; ih < nh; ih++) { + AliTrackReference* tr = (AliTrackReference*) fTmpTrackReferences.At(ih); + Int_t label = tr->Label(); + if (label == ip) { + Int_t nref = fTrackReferences.GetEntriesFast(); + new(fTrackReferences[nref]) AliTrackReference(*tr); + } + } + } + treeTR->Fill(); + fTrackReferences.Clear(); + ifills++; + } // tracks + // Check + if (ifills != stack->GetNtrack()) + AliWarning(Form("Number of entries in TreeTR (%5d) unequal to TreeK (%5d) \n", ifills, stack->GetNtrack())); +// +// Clean-up + delete fTmpTreeTR; + fTmpFileTR->Close(); + delete fTmpFileTR; + fTmpTrackReferences.Clear(); + gSystem->Exec("rm -rf TrackRefsTmp.root"); } +