X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliLego.cxx;h=2474f35dabd70e0c5bf802916b7cd9e36989d493;hb=69f2a1049cad647f6bed7c5f2c1833e1a79a16b0;hp=9a43a5b5ca71ac367fb3cb0402f293074ed13af9;hpb=65fb704d5a3c40dcad260b5c6ceb87aaa8697172;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliLego.cxx b/STEER/AliLego.cxx index 9a43a5b5ca7..2474f35dabd 100644 --- a/STEER/AliLego.cxx +++ b/STEER/AliLego.cxx @@ -13,50 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.19 2000/10/26 14:13:05 morsch -- Change from coordinates theta, phi to general coordinates Coor1 and Coor2. -- Lego generator instance can be passed in constructor. - -Revision 1.18 2000/10/02 21:28:14 fca -Removal of useless dependecies via forward declarations - -Revision 1.17 2000/07/12 08:56:25 fca -Coding convention correction and warning removal - -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$ */ ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// @@ -88,74 +45,157 @@ Introduction of the Copyright and cvs Log // ////////////////////////////////////////////////////////////// +#include "TClonesArray.h" +#include "TH2.h" #include "TMath.h" +#include "TString.h" +#include "TVirtualMC.h" +#include "AliLog.h" +#include "AliDebugVolume.h" #include "AliLego.h" #include "AliLegoGenerator.h" -#include "AliConst.h" #include "AliMC.h" -#include "TH2.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), + fStopped(0) { // // Default constructor // - fHistRadl = 0; - fHistAbso = 0; - fHistGcm2 = 0; - fHistReta = 0; } -//___________________________________________ +//_______________________________________________________________________ +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), + fStopped(0) +{ + // + // Copy constructor + // + lego.Copy(*this); +} + + +//_______________________________________________________________________ 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) + 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), + fStopped(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", ntheta,thetamin,thetamax,nphi,phimin,phimax); fHistAbso = new TH2F("habso","Interaction length map", ntheta,thetamin,thetamax,nphi,phimin,phimax); fHistGcm2 = new TH2F("hgcm2","g/cm2 length map", ntheta,thetamin,thetamax,nphi,phimin,phimax); +// + fVolumesFwd = new TClonesArray("AliDebugVolume",1000); + fVolumesBwd = new TClonesArray("AliDebugVolume",1000); + fDebug = AliDebugLevel(); } -AliLego::AliLego(const char *title, AliLegoGenerator* generator) - : TNamed("Lego Generator",title) +//_______________________________________________________________________ +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), + fStopped(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); + 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); + // + // - fHistRadl = new TH2F("hradl","Radiation length map", - n1, c1min, c1max, n2, c2min, c2max); - fHistAbso = new TH2F("habso","Interaction length map", - n1, c1min, c1max, n2, c2min, c2max); - fHistGcm2 = new TH2F("hgcm2","g/cm2 length map", - n1, c1min, c1max, n2, c2min, c2max); + fVolumesFwd = new TClonesArray("AliDebugVolume",1000); + fVolumesBwd = new TClonesArray("AliDebugVolume",1000); + fDebug = AliDebugLevel(); } -//___________________________________________ +//_______________________________________________________________________ AliLego::~AliLego() { // @@ -165,9 +205,11 @@ AliLego::~AliLego() delete fHistAbso; delete fHistGcm2; delete fGener; + delete fVolumesFwd; + delete fVolumesBwd; } -//___________________________________________ +//_______________________________________________________________________ void AliLego::BeginEvent() { // @@ -176,9 +218,20 @@ void AliLego::BeginEvent() fTotRadl = 0; fTotAbso = 0; fTotGcm2 = 0; + fStopped = 0; + + if (fDebug) { + if (fErrorCondition) ToAliDebug(1, DumpVolumes()); + fVolumesFwd->Delete(); + fVolumesBwd->Delete(); + fStepsForward = 0; + fStepsBackward = 0; + fErrorCondition = 0; + if (gAlice->GetMCApp()->GetCurrentTrackNumber() == 0) fStepBack = 0; + } } -//___________________________________________ +//_______________________________________________________________________ void AliLego::FinishEvent() { // @@ -187,12 +240,12 @@ void AliLego::FinishEvent() Double_t c1, c2; c1 = fGener->CurCoor1(); c2 = fGener->CurCoor2(); - fHistRadl->Fill(c1,c2,fTotRadl); - fHistAbso->Fill(c1,c2,fTotAbso); - fHistGcm2->Fill(c1,c2,fTotGcm2); + fHistRadl->Fill(c2,c1,fTotRadl); + fHistAbso->Fill(c2,c1,fTotAbso); + fHistGcm2->Fill(c2,c1,fTotGcm2); } -//___________________________________________ +//_______________________________________________________________________ void AliLego::FinishRun() { // @@ -206,77 +259,201 @@ void AliLego::FinishRun() fHistRadl->Delete(); fHistRadl=0; fHistAbso->Delete(); fHistAbso=0; fHistGcm2->Delete(); fHistGcm2=0; + // + if (fErrorCondition) ToAliError(DumpVolumes()); } -//___________________________________________ -void AliLego::Copy(AliLego &lego) const +//_______________________________________________________________________ +void AliLego::Copy(TObject&) const { // // Copy *this onto lego -- not implemented // - Fatal("Copy","Not implemented!\n"); + AliFatal("Not implemented!"); } -//___________________________________________ -void AliLego::StepManager() +//_______________________________________________________________________ +void AliLego::StepManager() { -// 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; - - 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; + // + // 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]; - 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; -// Int_t copy; -// Int_t id = gMC->CurrentVolID(copy); -// char* vol = gMC->VolName(id); - -// printf("\n %f %f %f %f %s ", fTotRadl, vect[0], vect[1], vect[2], vol); - - } - 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; -// Int_t copy; -// Int_t id = gMC->CurrentVolID(copy); -// char* vol = gMC->VolName(id); -// printf("\n %f %f %f %f %s %f ", fTotRadl, vect[0], vect[1], vect[2], vol, t); - } + TString tmp1, tmp2; + copy = 1; + id = gMC->CurrentVolID(copy); + vol = gMC->VolName(id); + Float_t step = gMC->TrackStep(); + + 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 + if (!fStopped) { + if (absl) fTotAbso += t/absl; + if (radl) fTotRadl += t/radl; + fTotGcm2 += t*dens; + } + +// printf("We will stop now %5d %13.3f !\n", fStopped, t); +// printf("%13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %s %13.3f\n", +// pos[2], TMath::Sqrt(pos[0] * pos[0] + pos[1] * pos[1]), step, a, z, radl, absl, gMC->CurrentVolName(), fTotRadl); + 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; + fStopped = kTRUE; + 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()); + + if(step) { + + if (absl) fTotAbso += step/absl; + if (radl) fTotRadl += step/radl; + fTotGcm2 += step*dens; +// printf("%13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %13.3f %s %13.3f\n", +// pos[2], TMath::Sqrt(pos[0] * pos[0] + pos[1] * pos[1]), step, a, z, radl, absl, gMC->CurrentVolName(), fTotRadl); + } + + } 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)) + { + AliWarning(Form("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"); +}