X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=STEER%2FAliMC.cxx;h=16a642ff6fb24f59d6d93427d12ecab0175e5efa;hb=8229c6583569f40405a7bccf29d69166205b12d4;hp=3fb4102314b027f07f432fb8ace0e4514b41a75f;hpb=2057aecb08b8cfad37f9af46c11de054cd39dd1a;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliMC.cxx b/STEER/AliMC.cxx index 3fb4102314b..16a642ff6fb 100644 --- a/STEER/AliMC.cxx +++ b/STEER/AliMC.cxx @@ -21,11 +21,16 @@ // Author: F.Carminati // Federico.Carminati@cern.ch +#include #include #include #include #include +#include +#include "TGeant3.h" + +#include "AliLog.h" #include "AliDetector.h" #include "AliGenerator.h" #include "AliHeader.h" @@ -34,6 +39,7 @@ #include "AliMCQA.h" #include "AliRun.h" #include "AliStack.h" +#include "AliMagF.h" #include "AliTrackReference.h" @@ -47,7 +53,9 @@ AliMC::AliMC() : fSum2Energy(0), fTrRmax(1.e10), fTrZmax(1.e10), - fDebug(0), + fRDecayMax(1.e10), + fRDecayMin(0), + fDecayPdg(0), fImedia(0), fTransParName("\0"), fMCQA(0), @@ -56,6 +64,7 @@ AliMC::AliMC() : { //default constructor + DecayLimits(); } //_______________________________________________________________________ @@ -67,7 +76,9 @@ AliMC::AliMC(const char *name, const char *title) : fSum2Energy(0), fTrRmax(1.e10), fTrZmax(1.e10), - fDebug(0), + fRDecayMax(1.e10), + fRDecayMin(0), + fDecayPdg(0), fImedia(new TArrayI(1000)), fTransParName("\0"), fMCQA(0), @@ -77,10 +88,9 @@ AliMC::AliMC(const char *name, const char *title) : //constructor // Set transport parameters SetTransPar(); - + DecayLimits(); // Prepare the tracking medium lists for(Int_t i=0;i<1000;i++) (*fImedia)[i]=-99; - } //_______________________________________________________________________ @@ -92,7 +102,9 @@ AliMC::AliMC(const AliMC &mc) : fSum2Energy(0), fTrRmax(1.e10), fTrZmax(1.e10), - fDebug(0), + fRDecayMax(1.e10), + fRDecayMin(0), + fDecayPdg(0), fImedia(0), fTransParName("\0"), fMCQA(0), @@ -100,7 +112,7 @@ AliMC::AliMC(const AliMC &mc) : fTrackReferences(0) { // - // Copy constructor for AliRun + // Copy constructor for AliMC // mc.Copy(*this); } @@ -126,28 +138,45 @@ AliMC::~AliMC() void AliMC::Copy(TObject &) const { //dummy Copy function - Fatal("Copy","Not implemented!\n"); + AliFatal("Not implemented!"); } //_______________________________________________________________________ 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 + const char *geomfilename = gAlice->GetGeometryFileName(); + if(gSystem->ExpandPathName(geomfilename)){ + AliInfo(Form("Loading geometry from file:\n %40s\n\n",geomfilename)); + TGeoManager::Import(geomfilename); + }else{ + AliInfo(Form("Geometry file %40s not found!\n",geomfilename)); + return; + } + }else{ + // Create modules, materials, 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())); } + } + } //_______________________________________________________________________ @@ -157,19 +186,20 @@ void AliMC::InitGeometry() // Initialize detectors and display geometry // - 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(); + // Initialise detector and display geometry + if(gAlice->IsRootGeometry()) detector->CreateMaterials(); + detector->Init(); + detector->BuildGeometry(); + AliInfo(Form("%10s R:%.2fs C:%.2fs", + detector->GetName(),stw.RealTime(),stw.CpuTime())); + } + } //_______________________________________________________________________ @@ -199,11 +229,11 @@ void AliMC::ResetGenerator(AliGenerator *generator) // if(fGenerator) if(generator) - Warning("ResetGenerator","Replacing generator %s with %s\n", - fGenerator->GetName(),generator->GetName()); + AliWarning(Form("Replacing generator %s with %s", + fGenerator->GetName(),generator->GetName())) else - Warning("ResetGenerator","Replacing generator %s with NULL\n", - fGenerator->GetName()); + AliWarning(Form("Replacing generator %s with NULL", + fGenerator->GetName())); fGenerator = generator; } @@ -211,15 +241,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()); } //_______________________________________________________________________ @@ -255,10 +282,22 @@ void AliMC::Stepping() // // Called at every step during transport // - + Int_t id = DetFromMate(gMC->GetMedium()); 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()) @@ -270,7 +309,7 @@ void AliMC::Stepping() // // write tracke reference for track which is dissapearing - MI if (gMC->IsTrackDisappeared()) { - if (gMC->Etot()>0.05) AddTrackReference(GetCurrentTrackNumber()); + if (gMC->Etot()>0.05) AddTrackReference(GetCurrentTrackNumber()); } //Call the appropriate stepping routine; @@ -354,16 +393,13 @@ 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(); @@ -376,7 +412,7 @@ 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 @@ -388,11 +424,11 @@ void AliMC::BeginEvent() if(gAlice->Lego() == 0x0) { - if (GetDebug()) Info("BeginEvent"," fRunLoader->MakeTree(K)"); + AliDebug(1, "fRunLoader->MakeTree(K)"); runloader->MakeTree("K"); } - if (GetDebug()) Info("BeginEvent"," gMC->SetStack(fRunLoader->Stack())"); + AliDebug(1, "gMC->SetStack(fRunLoader->Stack())"); gMC->SetStack(gAlice->GetRunLoader()->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 @@ -412,13 +448,13 @@ void AliMC::BeginEvent() } - 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()"); + AliDebug(1, "fRunLoader->MakeTrackRefsContainer()"); runloader->MakeTrackRefsContainer();//for insurance @@ -427,11 +463,11 @@ void AliMC::BeginEvent() AliModule *detector; while((detector = (AliModule*)next())) { - if (GetDebug()) Info("BeginEvent"," %s->MakeBranch(H)",detector->GetName()); + AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName())); detector->MakeBranch("H"); - if (GetDebug()) Info("BeginEvent"," %s->MakeBranchTR()",detector->GetName()); + AliDebug(2, Form("%s->MakeBranchTR()",detector->GetName())); detector->MakeBranchTR(); - if (GetDebug()) Info("BeginEvent"," %s->SetTreeAddress()",detector->GetName()); + AliDebug(2, Form("%s->SetTreeAddress()",detector->GetName())); detector->SetTreeAddress(); } // make branch for AliRun track References @@ -463,12 +499,12 @@ void AliMC::ResetHits() void AliMC::PostTrack() { // Posts tracks for each module - TObjArray &dets = *gAlice->Modules(); - AliModule *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(); } //_______________________________________________________________________ @@ -481,6 +517,10 @@ void AliMC::FinishPrimary() // static Int_t count=0; // const Int_t times=10; // This primary is finished, purify stack +#if ROOT_VERSION_CODE > 262152 + if (!(gMC->SecondariesAreOrdered())) + runloader->Stack()->ReorderKine(); +#endif runloader->Stack()->PurifyKine(); TIter next(gAlice->Modules()); @@ -529,7 +569,7 @@ void AliMC::FinishEvent() 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 @@ -549,7 +589,7 @@ void AliMC::FinishEvent() } else { - Error("FinishEvent","Can not get TreeE from RL"); + AliError("Can not get TreeE from RL"); } if(gAlice->Lego() == 0x0) @@ -559,16 +599,13 @@ void AliMC::FinishEvent() 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, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); } //_______________________________________________________________________ @@ -629,6 +666,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; @@ -640,8 +678,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 @@ -653,7 +691,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; } } @@ -788,7 +817,7 @@ void AliMC::SetTransPar(const char *filename) } //_______________________________________________________________________ -void AliMC::Browse(TBrowser *b) const +void AliMC::Browse(TBrowser *b) { // // Called when the item "Run" is clicked on the left pane @@ -980,7 +1009,7 @@ void AliMC::AddTrackReference(Int_t label){ // // add a trackrefernce to the list if (!fTrackReferences) { - cerr<<"Container trackrefernce not active\n"; + AliError("Container trackrefernce not active"); return; } Int_t nref = fTrackReferences->GetEntriesFast(); @@ -1026,3 +1055,52 @@ void AliMC::RemapTrackReferencesIDs(Int_t *map) } fTrackReferences->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 = gAlice->Field()->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 + if (fRDecayMax * fRDecayMax / rho / rho / 2. > 1.) return; + + // + + tmax = TMath::ACos(1. - fRDecayMax * fRDecayMax / rho / rho / 2.) / omega; // [ct] + tmin = TMath::ACos(1. - fRDecayMin * fRDecayMin / rho / rho / 2.) / 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 + // + TGeant3 * geant = (TGeant3*) gMC; + geant->Gcphys()->sumlif = t / p.Beta() / p.Gamma(); +}