X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliLego.cxx;h=81e9ccd70c858c39aec154dcf94d68d6ebd76d4a;hb=6667b60242f1be6df722f3d4d5efb64cc50fae75;hp=254a74c2a1322d95fc3b68d84e9ba82d2c56a2d4;hpb=8918e70084c0bd41164dfeb3795567b013b2693e;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliLego.cxx b/STEER/AliLego.cxx index 254a74c2a13..81e9ccd70c8 100644 --- a/STEER/AliLego.cxx +++ b/STEER/AliLego.cxx @@ -13,40 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.16 2000/05/26 08:35:03 fca -Move the check on z after z has been retrieved - -Revision 1.15 2000/05/16 13:10:40 fca -New method IsNewTrack and fix for a problem in Father-Daughter relations - -Revision 1.14 2000/04/27 10:38:21 fca -Correct termination of Lego Run and introduce Lego getter in AliRun - -Revision 1.13 2000/04/26 10:17:31 fca -Changes in Lego for G4 compatibility - -Revision 1.12 2000/04/07 11:12:33 fca -G4 compatibility changes - -Revision 1.11 2000/03/22 13:42:26 fca -SetGenerator does not replace an existing generator, ResetGenerator does - -Revision 1.10 2000/02/23 16:25:22 fca -AliVMC and AliGeant3 classes introduced -ReadEuclid moved from AliRun to AliModule - -Revision 1.9 1999/12/03 10:54:01 fca -Fix lego summary - -Revision 1.8 1999/10/01 09:54:33 fca -Correct logics for Lego StepManager - -Revision 1.7 1999/09/29 09:24:29 fca -Introduction of the Copyright and cvs Log - -*/ +/* $Id$ */ ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// @@ -78,57 +45,152 @@ Introduction of the Copyright and cvs Log // ////////////////////////////////////////////////////////////// +#include "TClonesArray.h" +#include "TH2.h" #include "TMath.h" +#include "TString.h" +#include "TVirtualMC.h" + +#include "AliDebugVolume.h" #include "AliLego.h" #include "AliLegoGenerator.h" -#include "AliRun.h" -#include "AliConst.h" #include "AliMC.h" +#include "AliRun.h" ClassImp(AliLego) - -//___________________________________________ -AliLego::AliLego() +//_______________________________________________________________________ +AliLego::AliLego(): + fGener(0), + fTotRadl(0), + fTotAbso(0), + fTotGcm2(0), + fHistRadl(0), + fHistAbso(0), + fHistGcm2(0), + fHistReta(0), + fVolumesFwd(0), + fVolumesBwd(0), + fStepBack(0), + fStepsBackward(0), + fStepsForward(0), + fErrorCondition(0), + fDebug(0) { // // Default constructor // - fHistRadl = 0; - fHistAbso = 0; - fHistGcm2 = 0; - fHistReta = 0; } -//___________________________________________ -AliLego::AliLego(const char *title, Int_t ntheta, Float_t themin, - Float_t themax, Int_t nphi, Float_t phimin, Float_t phimax, - Float_t rmin, Float_t rmax, Float_t zmax) - : TNamed("Lego Generator",title) +//_______________________________________________________________________ +AliLego::AliLego(const AliLego &lego): + TNamed(lego), + fGener(0), + fTotRadl(0), + fTotAbso(0), + fTotGcm2(0), + fHistRadl(0), + fHistAbso(0), + fHistGcm2(0), + fHistReta(0), + fVolumesFwd(0), + fVolumesBwd(0), + fStepBack(0), + fStepsBackward(0), + fStepsForward(0), + fErrorCondition(0), + fDebug(0) { // - // specify the angular limits and the size of the rectangular box + // Copy constructor // - fGener = new AliLegoGenerator(ntheta, themin, themax, - nphi, phimin, phimax, rmin, rmax, zmax); - - gAlice->ResetGenerator(fGener); + lego.Copy(*this); +} - Float_t etamin = -TMath::Log(TMath::Tan(TMath::Min((Double_t)themax*kDegrad/2,TMath::Pi()/2-1.e-10))); - Float_t etamax = -TMath::Log(TMath::Tan(TMath::Max((Double_t)themin*kDegrad/2, 1.e-10))); +//_______________________________________________________________________ +AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin, + Float_t thetamax, Int_t nphi, Float_t phimin, Float_t phimax, + Float_t rmin, Float_t rmax, Float_t zmax): + TNamed("Lego Generator",title), + fGener(0), + fTotRadl(0), + fTotAbso(0), + fTotGcm2(0), + fHistRadl(0), + fHistAbso(0), + fHistGcm2(0), + fHistReta(0), + fVolumesFwd(0), + fVolumesBwd(0), + fStepBack(0), + fStepsBackward(0), + fStepsForward(0), + fErrorCondition(0), + fDebug(0) +{ + // + // specify the angular limits and the size of the rectangular box + // + fGener = new AliLegoGenerator(ntheta, thetamin, thetamax, + nphi, phimin, phimax, rmin, rmax, zmax); fHistRadl = new TH2F("hradl","Radiation length map", - nphi,phimin,phimax,ntheta,themin,themax); + ntheta,thetamin,thetamax,nphi,phimin,phimax); fHistAbso = new TH2F("habso","Interaction length map", - nphi,phimin,phimax,ntheta,themin,themax); + ntheta,thetamin,thetamax,nphi,phimin,phimax); fHistGcm2 = new TH2F("hgcm2","g/cm2 length map", - nphi,phimin,phimax,ntheta,themin,themax); - fHistReta = new TH2F("hetar","Radiation length vs. eta", - nphi,phimin,phimax,ntheta,etamin,etamax); + ntheta,thetamin,thetamax,nphi,phimin,phimax); +// + fVolumesFwd = new TClonesArray("AliDebugVolume",1000); + fVolumesBwd = new TClonesArray("AliDebugVolume",1000); + fDebug = gAlice->GetDebug(); +} + +//_______________________________________________________________________ +AliLego::AliLego(const char *title, AliLegoGenerator* generator): + TNamed("Lego Generator",title), + fGener(0), + fTotRadl(0), + fTotAbso(0), + fTotGcm2(0), + fHistRadl(0), + fHistAbso(0), + fHistGcm2(0), + fHistReta(0), + fVolumesFwd(0), + fVolumesBwd(0), + fStepBack(0), + fStepsBackward(0), + fStepsForward(0), + fErrorCondition(0), + fDebug(0) +{ + // + // specify the angular limits and the size of the rectangular box + // + fGener = generator; + Float_t c1min, c1max, c2min, c2max; + Int_t n1 = fGener->NCoor1(); + Int_t n2 = fGener->NCoor2(); + + fGener->Coor1Range(c1min, c1max); + fGener->Coor2Range(c2min, c2max); + + fHistRadl = new TH2F("hradl","Radiation length map", + n2, c2min, c2max, n1, c1min, c1max); + fHistAbso = new TH2F("habso","Interaction length map", + n2, c2min, c2max, n1, c1min, c1max); + fHistGcm2 = new TH2F("hgcm2","g/cm2 length map", + n2, c2min, c2max, n1, c1min, c1max); + // + // + fVolumesFwd = new TClonesArray("AliDebugVolume",1000); + fVolumesBwd = new TClonesArray("AliDebugVolume",1000); + fDebug = gAlice->GetDebug(); } -//___________________________________________ +//_______________________________________________________________________ AliLego::~AliLego() { // @@ -137,12 +199,12 @@ AliLego::~AliLego() delete fHistRadl; delete fHistAbso; delete fHistGcm2; - delete fHistReta; - gAlice->ResetGenerator(0); delete fGener; + delete fVolumesFwd; + delete fVolumesBwd; } -//___________________________________________ +//_______________________________________________________________________ void AliLego::BeginEvent() { // @@ -151,28 +213,33 @@ void AliLego::BeginEvent() fTotRadl = 0; fTotAbso = 0; fTotGcm2 = 0; + + if (fDebug) { + if (fErrorCondition) DumpVolumes(); + fVolumesFwd->Delete(); + fVolumesBwd->Delete(); + fStepsForward = 0; + fStepsBackward = 0; + fErrorCondition = 0; + if (gAlice->GetMCApp()->GetCurrentTrackNumber() == 0) fStepBack = 0; + } } -//___________________________________________ +//_______________________________________________________________________ void AliLego::FinishEvent() { // // Finish the event and update the histos // - Double_t thed, phid, eta; - thed = fGener->CurTheta()*kRaddeg; - phid = fGener->CurPhi()*kRaddeg; - eta = -TMath::Log(TMath::Tan(TMath::Max( - TMath::Min((Double_t)(fGener->CurTheta())/2, - TMath::Pi()/2-1.e-10),1.e-10))); - - fHistRadl->Fill(phid,thed,fTotRadl); - fHistAbso->Fill(phid,thed,fTotAbso); - fHistGcm2->Fill(phid,thed,fTotGcm2); - fHistReta->Fill(phid,eta,fTotRadl); + Double_t c1, c2; + c1 = fGener->CurCoor1(); + c2 = fGener->CurCoor2(); + fHistRadl->Fill(c2,c1,fTotRadl); + fHistAbso->Fill(c2,c1,fTotAbso); + fHistGcm2->Fill(c2,c1,fTotGcm2); } -//___________________________________________ +//_______________________________________________________________________ void AliLego::FinishRun() { // @@ -181,18 +248,17 @@ void AliLego::FinishRun() fHistRadl->Write(); fHistAbso->Write(); fHistGcm2->Write(); - fHistReta->Write(); // Delete histograms from memory fHistRadl->Delete(); fHistRadl=0; fHistAbso->Delete(); fHistAbso=0; fHistGcm2->Delete(); fHistGcm2=0; - fHistReta->Delete(); fHistReta=0; - + // + if (fErrorCondition) DumpVolumes(); } -//___________________________________________ -void AliLego::Copy(AliLego &lego) const +//_______________________________________________________________________ +void AliLego::Copy(TObject&) const { // // Copy *this onto lego -- not implemented @@ -200,52 +266,178 @@ void AliLego::Copy(AliLego &lego) const Fatal("Copy","Not implemented!\n"); } -//___________________________________________ -void AliLego::StepManager() +//_______________________________________________________________________ +void AliLego::StepManager() { -// called from AliRun::Stepmanager from gustep. -// Accumulate the 3 parameters step by step + // + // called from AliRun::Stepmanager from gustep. + // Accumulate the 3 parameters step by step + // + static Float_t t; + Float_t a,z,dens,radl,absl; + Int_t i, id, copy; + const char* vol; + static Float_t vect[3], dir[3]; + + TString tmp1, tmp2; + copy = 1; + id = gMC->CurrentVolID(copy); + vol = gMC->VolName(id); + Float_t step = gMC->TrackStep(); - static Float_t t; - Float_t a,z,dens,radl,absl; - Int_t i; - - Float_t step = gMC->TrackStep(); - - Float_t vect[3], dir[3]; - TLorentzVector pos, mom; - - gMC->CurrentMaterial(a,z,dens,radl,absl); - - if (z < 1) return; + TLorentzVector pos, mom; + gMC->TrackPosition(pos); + gMC->TrackMomentum(mom); + + Int_t status = 0; + if (gMC->IsTrackEntering()) status = 1; + if (gMC->IsTrackExiting()) status = 2; + + if (! fStepBack) { + // printf("\n volume %s %d", vol, status); + // + // Normal Forward stepping + // + if (fDebug) { + // printf("\n steps fwd %d %s %d %f", fStepsForward, vol, fErrorCondition, step); + + // + // store volume if different from previous + // + + TClonesArray &lvols = *fVolumesFwd; + if (fStepsForward > 0) { + AliDebugVolume* tmp = dynamic_cast((*fVolumesFwd)[fStepsForward-1]); + if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) { + fStepsForward -= 2; + fVolumesFwd->RemoveAt(fStepsForward); + fVolumesFwd->RemoveAt(fStepsForward+1); + } + } + + new(lvols[fStepsForward++]) + AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status); + + } // Debug + // + // Get current material properties + + gMC->CurrentMaterial(a,z,dens,radl,absl); + + if (z < 1) return; + + // --- See if we have to stop now + if (TMath::Abs(pos[2]) > fGener->ZMax() || + pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) { + if (!gMC->IsNewTrack()) { + // Not the first step, add past contribution + fTotAbso += t/absl; + fTotRadl += t/radl; + fTotGcm2 += t*dens; + if (fDebug) { + // + // generate "mirror" particle flying back + // + fStepsBackward = fStepsForward; + + Float_t pmom[3], orig[3]; + Float_t polar[3] = {0.,0.,0.}; + Int_t ntr; + pmom[0] = -dir[0]; + pmom[1] = -dir[1]; + pmom[2] = -dir[2]; + orig[0] = vect[0]; + orig[1] = vect[1]; + orig[2] = vect[2]; + + gAlice->GetMCApp()->PushTrack(1, gAlice->GetMCApp()->GetCurrentTrackNumber(), + 0, pmom, orig, polar, 0., kPNoProcess, ntr); + } // debug + + } // not a new track ! + + if (fDebug) fStepBack = 1; + gMC->StopTrack(); + return; + } // outside scoring region ? + + // --- See how long we have to go + for(i=0;i<3;++i) { + vect[i]=pos[i]; + dir[i]=mom[i]; + } + + t = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax()); - gMC->TrackPosition(pos); - gMC->TrackMomentum(mom); -// --- See if we have to stop now - if (TMath::Abs(pos[2]) > fGener->ZMax() || - pos[0]*pos[0] +pos[1]*pos[1] > fGener->RadMax()*fGener->RadMax()) { - if (!gMC->IsNewTrack()) { - // Not the first step, add past contribution - fTotAbso += t/absl; - fTotRadl += t/radl; - fTotGcm2 += t*dens; - } - gMC->StopTrack(); - return; - } - -// --- See how long we have to go - for(i=0;i<3;++i) { - vect[i]=pos[i]; - dir[i]=mom[i]; - } - - t = fGener->PropagateCylinder(vect,dir,fGener->RadMax(),fGener->ZMax()); - - if(step) { - fTotAbso += step/absl; - fTotRadl += step/radl; - fTotGcm2 += step*dens; - } + if(step) { + + fTotAbso += step/absl; + fTotRadl += step/radl; + fTotGcm2 += step*dens; + } + + } else { + if (fDebug) { + // + // Geometry debugging + // Fly back and compare volume sequence + // + TClonesArray &lvols = *fVolumesBwd; + if (fStepsBackward < fStepsForward) { + AliDebugVolume* tmp = dynamic_cast((*fVolumesBwd)[fStepsBackward]); + if (tmp->IsVEqual(vol, copy) && gMC->IsTrackEntering()) { + fStepsBackward += 2; + fVolumesBwd->RemoveAt(fStepsBackward-1); + fVolumesBwd->RemoveAt(fStepsBackward-2); + } + } + + fStepsBackward--; + // printf("\n steps %d %s %d", fStepsBackward, vol, fErrorCondition); + if (fStepsBackward < 0) { + gMC->StopTrack(); + fStepBack = 0; + return; + } + + new(lvols[fStepsBackward]) AliDebugVolume(vol,copy,step,pos[0], pos[1], pos[2], status); + + AliDebugVolume* tmp = dynamic_cast((*fVolumesFwd)[fStepsBackward]); + if (! (tmp->IsVEqual(vol, copy)) && (!fErrorCondition)) + { + printf("\n Problem at (x,y,z): %d %f %f %f, volumes: %s %s step: %f\n", + fStepsBackward, pos[0], pos[1], pos[2], tmp->GetName(), vol, step); + fErrorCondition = 1; + } + } // Debug + } // bwd/fwd } +//_______________________________________________________________________ +void AliLego::DumpVolumes() +{ + // + // Dump volume sequence in case of error + // + printf("\n Dumping Volume Sequence:"); + printf("\n =============================================="); + + for (Int_t i = fStepsForward-1; i>=0; i--) + { + AliDebugVolume* tmp1 = dynamic_cast((*fVolumesFwd)[i]); + AliDebugVolume* tmp2 = dynamic_cast((*fVolumesBwd)[i]); + if (tmp1) + printf("\n Volume Fwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n" + , i, + tmp1->GetName(), tmp1->CopyNumber(), tmp1->Step(), + tmp1->X(), tmp1->Y(), tmp1->Z(), tmp1->Status()); + if (tmp2 && i>= fStepsBackward) + printf("\n Volume Bwd: %3d: %5s (%3d) step: %12.5e (x,y,z) (%12.5e %12.5e %12.5e) status: %9s \n" + , i, + tmp2->GetName(), tmp2->CopyNumber(), tmp2->Step(), + tmp2->X(), tmp2->Y(), tmp2->Z(), tmp2->Status()); + + printf("\n ............................................................................\n"); + } + printf("\n ==============================================\n"); +}