X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliLego.cxx;h=2474f35dabd70e0c5bf802916b7cd9e36989d493;hb=63a1afff7e98dc012a9aad8c431d820506a1a14d;hp=30076cd50c223bbc43f4b70a66358fb16723b227;hpb=19237f5fa90f951ab222b61d4f3ad8be3796f9d5;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliLego.cxx b/STEER/AliLego.cxx index 30076cd50c2..2474f35dabd 100644 --- a/STEER/AliLego.cxx +++ b/STEER/AliLego.cxx @@ -13,64 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.23 2001/05/20 10:10:39 morsch -- Debug output at the beginning of new event and end of run. -- Filter out boundary loops. - -Revision 1.22 2001/05/11 13:22:40 morsch -If run with debug option (from gAlice) geantinos are sent back and volume sequence forward/backward is compared. -Can be very verbous in some cases. - -Revision 1.21 2000/12/15 10:33:59 morsch -Invert coordinates to make meaningful zylindrical plots. - -Revision 1.20 2000/11/30 07:12:48 alibrary -Introducing new Rndm and QA classes - -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$ */ ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// @@ -102,49 +45,99 @@ Introduction of the Copyright and cvs Log // ////////////////////////////////////////////////////////////// +#include "TClonesArray.h" +#include "TH2.h" #include "TMath.h" +#include "TString.h" +#include "TVirtualMC.h" -#include "AliLego.h" - +#include "AliLog.h" #include "AliDebugVolume.h" -#include "AliRun.h" +#include "AliLego.h" #include "AliLegoGenerator.h" -#include "AliConst.h" #include "AliMC.h" -#include "TH2.h" -#include "../TGeant3/TGeant3.h" -#include "TString.h" -#include "TClonesArray.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", @@ -152,49 +145,57 @@ AliLego::AliLego(const char *title, Int_t ntheta, Float_t thetamin, fHistGcm2 = new TH2F("hgcm2","g/cm2 length map", ntheta,thetamin,thetamax,nphi,phimin,phimax); // - fStepBack = 0; - fStepsBackward = 0; - fStepsForward = 0; - fStepBack = 0; fVolumesFwd = new TClonesArray("AliDebugVolume",1000); fVolumesBwd = new TClonesArray("AliDebugVolume",1000); - fDebug = gAlice->GetDebug(); - fErrorCondition = 0; + 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", - 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); -// -// - fStepBack = 0; - fStepsBackward = 0; - fStepsForward = 0; - fStepBack = 0; - fVolumesFwd = new TClonesArray("AliDebugVolume",1000); - fVolumesBwd = new TClonesArray("AliDebugVolume",1000); - fDebug = gAlice->GetDebug(); - fErrorCondition =0; + fVolumesFwd = new TClonesArray("AliDebugVolume",1000); + fVolumesBwd = new TClonesArray("AliDebugVolume",1000); + fDebug = AliDebugLevel(); } -//___________________________________________ +//_______________________________________________________________________ AliLego::~AliLego() { // @@ -204,11 +205,11 @@ AliLego::~AliLego() delete fHistAbso; delete fHistGcm2; delete fGener; - if (fVolumesFwd) delete fVolumesFwd; - if (fVolumesBwd) delete fVolumesBwd; + delete fVolumesFwd; + delete fVolumesBwd; } -//___________________________________________ +//_______________________________________________________________________ void AliLego::BeginEvent() { // @@ -217,19 +218,20 @@ void AliLego::BeginEvent() fTotRadl = 0; fTotAbso = 0; fTotGcm2 = 0; - + fStopped = 0; + if (fDebug) { - if (fErrorCondition) DumpVolumes(); - fVolumesFwd->Delete(); - fVolumesBwd->Delete(); - fStepsForward = 0; - fStepsBackward = 0; - fErrorCondition = 0; - if (gAlice->CurrentTrack() == 0) fStepBack = 0; + 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() { // @@ -243,7 +245,7 @@ void AliLego::FinishEvent() fHistGcm2->Fill(c2,c1,fTotGcm2); } -//___________________________________________ +//_______________________________________________________________________ void AliLego::FinishRun() { // @@ -258,23 +260,25 @@ void AliLego::FinishRun() fHistAbso->Delete(); fHistAbso=0; fHistGcm2->Delete(); fHistGcm2=0; // - if (fErrorCondition) DumpVolumes(); + 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() { -// called from AliRun::Stepmanager from gustep. -// Accumulate the 3 parameters step by step - +//_______________________________________________________________________ +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, id, copy; @@ -286,171 +290,170 @@ void AliLego::StepManager() { 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 = (AliDebugVolume*) (*fVolumesFwd)[fStepsForward-1]; - if (tmp->IsEqual(vol, copy) && gMC->IsTrackEntering()) { - step += ((AliDebugVolume*) lvols[fStepsForward])->Step(); - fStepsForward -= 2; - delete ((AliDebugVolume*) lvols[fStepsForward]); - delete ((AliDebugVolume*) lvols[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; + + 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); -// --- 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->SetTrack(1, gAlice->CurrentTrack(), - 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()); - - 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 = (AliDebugVolume*) (*fVolumesBwd)[fStepsBackward]; - if (tmp->IsEqual(vol, copy) && gMC->IsTrackEntering()) { - fStepsBackward += 2; - delete ((AliDebugVolume*) lvols[fStepsBackward-1]); - delete ((AliDebugVolume*) lvols[fStepsBackward-2]); - step += ((AliDebugVolume*) lvols[fStepsBackward])->Step(); - } - } - - 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); + // + // store volume if different from previous + // - AliDebugVolume* tmp = (AliDebugVolume*) (*fVolumesFwd)[fStepsBackward]; - if (! (tmp->IsEqual(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 + 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--) + // + // 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 = (AliDebugVolume*) (*fVolumesFwd)[i]; - AliDebugVolume* tmp2 = (AliDebugVolume*) (*fVolumesBwd)[i]; - 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()); - - 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"); + 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"); + printf("\n ==============================================\n"); } - - - - - - -