X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliMC.cxx;h=f18473ffadeefdc00988a945c8575dbc9472cc77;hb=29d9710ea9413e67b964c0cb2f7d7bc32f211304;hp=531be4140ec0da60de9221ec7ff2340d369d3a89;hpb=ced249e6cb0ac19ef0ef8d5dadbc6e1a30f43922;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliMC.cxx b/STEER/AliMC.cxx index 531be4140ec..f18473ffade 100644 --- a/STEER/AliMC.cxx +++ b/STEER/AliMC.cxx @@ -21,32 +21,36 @@ // Author: F.Carminati // Federico.Carminati@cern.ch +#include + #include -#include +#include #include +#include +#include #include +#include +#include #include #include #include -#include -#include "AliLog.h" +#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 "AliMagF.h" #include "AliTrackReference.h" -#include "AliSimulation.h" -#include "AliGeomManager.h" -#include "AliCDBManager.h" -#include "AliCDBStorage.h" -#include "AliCDBEntry.h" - ClassImp(AliMC) @@ -63,9 +67,11 @@ AliMC::AliMC() : fDecayPdg(0), fImedia(0), fTransParName("\0"), - fMCQA(0), fHitLists(0), - fTrackReferences(0) + fTmpTreeTR(0), + fTmpFileTR(0), + fTrackReferences(), + fTmpTrackReferences() { //default constructor @@ -86,9 +92,11 @@ AliMC::AliMC(const char *name, const char *title) : 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 @@ -98,52 +106,14 @@ AliMC::AliMC(const char *name, const char *title) : 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), - fRDecayMax(1.e10), - fRDecayMin(-1.), - fDecayPdg(0), - fImedia(0), - fTransParName("\0"), - fMCQA(0), - fHitLists(0), - fTrackReferences(0) -{ - // - // Copy constructor for AliMC - // - 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(TObject &) const -{ - //dummy Copy function - AliFatal("Not implemented!"); } //_______________________________________________________________________ @@ -156,8 +126,8 @@ void AliMC::ConstructGeometry() // at InitGeometry(). // - if(gAlice->IsRootGeometry()){ //load geometry either from CDB or from file - if(gAlice->IsGeomFromCDB()){ + if(AliSimulation::Instance()->IsGeometryFromFile()){ //load geometry either from CDB or from file + if(IsGeometryFromCDB()){ AliInfo("Loading geometry from CDB default storage"); AliCDBPath path("GRP","Geometry","Data"); AliCDBEntry *entry=AliCDBManager::Instance()->Get(path.GetPath()); @@ -167,7 +137,7 @@ void AliMC::ConstructGeometry() if (!gGeoManager) AliFatal("TGeoManager object not found in the specified CDB entry!"); }else{ // Load geometry - const char *geomfilename = gAlice->GetGeometryFileName(); + const char *geomfilename = AliSimulation::Instance()->GetGeometryFile(); if(gSystem->ExpandPathName(geomfilename)){ AliInfo(Form("Loading geometry from file:\n %40s",geomfilename)); TGeoManager::Import(geomfilename); @@ -176,8 +146,10 @@ void AliMC::ConstructGeometry() return; } } + gMC->SetRootGeometry(); }else{ // Create modules, materials, geometry + if (!gGeoManager) new TGeoManager("ALICE", "ALICE geometry"); TStopwatch stw; TIter next(gAlice->Modules()); AliModule *detector; @@ -197,16 +169,19 @@ void AliMC::ConstructGeometry() //_______________________________________________________________________ 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::GetInstance()) return kFALSE; - AliGeomManager::SetGeometry(gGeoManager); - return AliSimulation::GetInstance()->MisalignGeometry(gAlice->GetRunLoader()); + // Call misalignment code if AliSimulation object was defined. + + if(!AliSimulation::Instance()->IsGeometryFromFile()){ + //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()); } //_______________________________________________________________________ @@ -220,6 +195,8 @@ void AliMC::ConstructOpGeometry() AliModule *detector; AliInfo("Optical properties definition"); while((detector = dynamic_cast(next()))) { + // Initialise detector geometry + if(AliSimulation::Instance()->IsGeometryFromFile()) detector->CreateMaterials(); // Initialise detector optical properties detector->DefineOpticalProperties(); } @@ -238,14 +215,30 @@ void AliMC::InitGeometry() AliModule *detector; while((detector = dynamic_cast(next()))) { stw.Start(); - // Initialise detector geometry - if(gAlice->IsRootGeometry()) detector->CreateMaterials(); detector->Init(); AliInfo(Form("%10s R:%.2fs C:%.2fs", detector->GetName(),stw.RealTime(),stw.CpuTime())); } } +//_______________________________________________________________________ +void AliMC::SetGeometryFromCDB() +{ + // Set the loading of geometry from cdb instead of creating it + // A default CDB storage needs to be set before this method is called + if(AliCDBManager::Instance()->IsDefaultStorageSet() && + AliCDBManager::Instance()->GetRun() >= 0) + AliSimulation::Instance()->SetGeometryFile("*OCDB*"); + else + AliError("Loading of geometry from CDB ignored. First set a default CDB storage!"); +} + +//_______________________________________________________________________ +Bool_t AliMC::IsGeometryFromCDB() const +{ + return (strcmp(AliSimulation::Instance()->GetGeometryFile(),"*OCDB*")==0); +} + //_______________________________________________________________________ void AliMC::SetAllAlignableVolumes() { @@ -322,21 +315,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(); } //_______________________________________________________________________ @@ -345,7 +336,6 @@ void AliMC::Stepping() // // Called at every step during transport // - Int_t id = DetFromMate(gMC->CurrentMedium()); if (id < 0) return; @@ -363,22 +353,25 @@ void AliMC::Stepping() // // --- 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()) { - if (gMC->Etot()>0.05) AddTrackReference(GetCurrentTrackNumber()); + + 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()) { - if(AliLog::GetGlobalDebugLevel()>0) fMCQA->StepManager(id); det->StepManager(); } } @@ -387,7 +380,7 @@ void AliMC::Stepping() //_______________________________________________________________________ void AliMC::EnergySummary() { - // + //e // Print summary of deposited energy // @@ -396,7 +389,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) { @@ -464,7 +457,7 @@ void AliMC::BeginEvent() AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); AliDebug(1, ">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"); - AliRunLoader *runloader=gAlice->GetRunLoader(); + AliRunLoader *runloader=AliRunLoader::Instance(); /*******************************/ /* Clean after eventual */ @@ -481,18 +474,19 @@ void AliMC::BeginEvent() // 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(AliSimulation::Instance()->Lego() == 0x0) + { + AliDebug(1, "fRunLoader->MakeTree(K)"); + runloader->MakeTree("K"); + } - if(gAlice->Lego() == 0x0) - { - AliDebug(1, "fRunLoader->MakeTree(K)"); - runloader->MakeTree("K"); - } - AliDebug(1, "gMC->SetStack(fRunLoader->Stack())"); - gMC->SetStack(gAlice->GetRunLoader()->Stack());//Was in InitMC - but was moved here + 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 @@ -500,15 +494,15 @@ 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; + } AliDebug(1, "ResetHits()"); @@ -517,32 +511,17 @@ void AliMC::BeginEvent() AliDebug(1, "fRunLoader->MakeTree(H)"); runloader->MakeTree("H"); - AliDebug(1, "fRunLoader->MakeTrackRefsContainer()"); - runloader->MakeTrackRefsContainer();//for insurance + MakeTmpTrackRefsTree(); //create new branches and SetAdresses TIter next(gAlice->Modules()); AliModule *detector; while((detector = (AliModule*)next())) { - AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName())); - detector->MakeBranch("H"); - AliDebug(2, Form("%s->MakeBranchTR()",detector->GetName())); - detector->MakeBranchTR(); - AliDebug(2, Form("%s->SetTreeAddress()",detector->GetName())); - detector->SetTreeAddress(); + AliDebug(2, Form("%s->MakeBranch(H)",detector->GetName())); + detector->MakeBranch("H"); } - // make branch for AliRun track References - TTree * treeTR = runloader->TreeTR(); - if (treeTR){ - // make branch for central track references - if (!fTrackReferences) fTrackReferences = new TClonesArray("AliTrackReference",0); - TBranch *branch; - branch = treeTR->Branch("AliRun",&fTrackReferences); - branch->SetAddress(&fTrackReferences); - } - // } //_______________________________________________________________________ @@ -558,10 +537,37 @@ 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() { // Posts tracks for each module + TObjArray &dets = *gAlice->Modules(); AliModule *module; @@ -576,15 +582,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 #if ROOT_VERSION_CODE > 262152 - if (!(gMC->SecondariesAreOrdered())) - runloader->Stack()->ReorderKine(); + if (!(gMC->SecondariesAreOrdered())) { + if (runloader->Stack()->ReorderKine()) RemapHits(); + } #endif - runloader->Stack()->PurifyKine(); + if (runloader->Stack()->PurifyKine()) RemapHits(); TIter next(gAlice->Modules()); AliModule *detector; @@ -599,7 +607,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()); } //_______________________________________________________________________ @@ -608,10 +651,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; @@ -627,7 +668,7 @@ 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(); @@ -640,9 +681,11 @@ void AliMC::FinishEvent() header->SetNprimary(stack->GetNprimary()); header->SetNtrack(stack->GetNtrack()); - // Write out the kinematics - if (!gAlice->Lego()) 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(); @@ -656,7 +699,7 @@ void AliMC::FinishEvent() AliError("Can not get TreeE from RL"); } - if(gAlice->Lego() == 0x0) + if(AliSimulation::Instance()->Lego() == 0x0) { runloader->WriteKinematics("OVERWRITE"); runloader->WriteTrackRefs("OVERWRITE"); @@ -672,13 +715,6 @@ void AliMC::FinishEvent() AliDebug(1, "<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<"); } -//_______________________________________________________________________ -void AliMC::Field(const Double_t* x, Double_t* b) const -{ - // Calculates field "b" at point "x" - gAlice->Field(x,b); -} - //_______________________________________________________________________ void AliMC::Init() { @@ -703,9 +739,6 @@ void AliMC::Init() fSummEnergy.Set(gMC->NofVolumes()+1); fSum2Energy.Set(gMC->NofVolumes()+1); - // - fMCQA = new AliMCQA(gAlice->GetNdets()); - // Register MC in configuration AliConfig::Instance()->Add(gMC); @@ -782,12 +815,12 @@ void AliMC::ReadTransPar() const Int_t kncuts=10; - const Int_t knflags=11; + const Int_t knflags=12; const Int_t knpars=kncuts+knflags; const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO", "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI", "BREM","COMP","DCAY","DRAY","HADR","LOSS", - "MULS","PAIR","PHOT","RAYL"}; + "MULS","PAIR","PHOT","RAYL","STRA"}; char line[256]; char detName[7]; char* filtmp; @@ -810,9 +843,9 @@ void AliMC::ReadTransPar() for(i=0;iAdd(fMCQA,"AliMCQA"); -} - //_______________________________________________________________________ void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const { @@ -900,7 +921,7 @@ void AliMC::AddHit(Int_t id, Int_t track, Int_t *vol, Float_t *hits) const // Add a hit to detector id // TObjArray &dets = *gAlice->Modules(); - if(dets[id]) dynamic_cast(dets[id])->AddHit(track,vol,hits); + if(dets[id]) static_cast(dets[id])->AddHit(track,vol,hits); } //_______________________________________________________________________ @@ -910,7 +931,7 @@ void AliMC::AddDigit(Int_t id, Int_t *tracks, Int_t *digits) const // Add digit to detector id // TObjArray &dets = *gAlice->Modules(); - if(dets[id]) dynamic_cast(dets[id])->AddDigit(tracks,digits); + if(dets[id]) static_cast(dets[id])->AddDigit(tracks,digits); } //_______________________________________________________________________ @@ -918,7 +939,7 @@ Int_t AliMC::GetCurrentTrackNumber() const { // // Returns current track // - return gAlice->GetRunLoader()->Stack()->GetCurrentTrackNumber(); + return AliRunLoader::Instance()->Stack()->GetCurrentTrackNumber(); } //_______________________________________________________________________ @@ -927,7 +948,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); } @@ -938,7 +959,7 @@ void AliMC::DumpPStack () const // // Dumps the particle stack // - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader->Stack()) runloader->Stack()->DumpPStack(); } @@ -949,7 +970,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; @@ -962,7 +983,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; @@ -973,7 +994,7 @@ TParticle* AliMC::Particle(Int_t i) const { // Returns the i-th particle from the stack taking into account // the remaping done by PurifyKine - AliRunLoader * runloader = gAlice->GetRunLoader(); + AliRunLoader * runloader = AliRunLoader::Instance(); if (runloader) if (runloader->Stack()) return runloader->Stack()->Particle(i); @@ -981,11 +1002,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(); @@ -993,13 +1014,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, +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, @@ -1015,7 +1036,7 @@ void AliMC::PushTrack(Int_t done, Int_t parent, Int_t pdg, { // 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, @@ -1027,7 +1048,7 @@ 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); @@ -1039,7 +1060,7 @@ 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); @@ -1050,7 +1071,7 @@ 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); @@ -1062,23 +1083,22 @@ 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) { - AliError("Container trackrefernce not active"); - 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); } @@ -1089,37 +1109,33 @@ 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) +//_______________________________________________________________________ +void AliMC::RemapTrackReferencesIDs(const Int_t *map) { // // Remapping track reference // Called at finish primary // - if (!fTrackReferences) return; - Int_t nEntries = fTrackReferences->GetEntries(); + + Int_t nEntries = fTmpTrackReferences.GetEntries(); for (Int_t i=0; i < nEntries; i++){ - AliTrackReference * ref = dynamic_cast(fTrackReferences->UncheckedAt(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); - fTrackReferences->RemoveAt(i); + fTmpTrackReferences.RemoveAt(i); } } // if ref } - fTrackReferences->Compress(); + fTmpTrackReferences.Compress(); } +//_______________________________________________________________________ void AliMC::FixParticleDecaytime() { // @@ -1134,7 +1150,7 @@ void AliMC::FixParticleDecaytime() // Transverse velocity Double_t vt = p.Pt() / p.E(); - if ((b = gAlice->Field()->SolenoidField()) > 0.) { // [kG] + if ((b = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->SolenoidField()) > 0.) { // [kG] // Radius of helix @@ -1147,12 +1163,13 @@ void AliMC::FixParticleDecaytime() // Maximum and minimum decay time // // Check for curlers first - if (fRDecayMax * fRDecayMax / rho / rho / 2. > 1.) return; + const Double_t kOvRhoSqr2 = 1./(rho*TMath::Sqrt(2.)); + if (fRDecayMax * kOvRhoSqr2 > 1.) return; // - tmax = TMath::ACos(1. - fRDecayMax * fRDecayMax / rho / rho / 2.) / omega; // [ct] - tmin = TMath::ACos(1. - fRDecayMin * fRDecayMin / rho / rho / 2.) / omega; // [ct] + 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] @@ -1167,3 +1184,127 @@ void AliMC::FixParticleDecaytime() // 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", &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(); + // make branch for central track references + TClonesArray* pRef = &fTrackReferences; + treeTR->Branch("TrackReferences", &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"); +} +