X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliLego.cxx;h=5e5e53563ec6418279bf06ddaf13217498cf8920;hb=daa897f53dc01128202a7e3d9fff23fcedae5db7;hp=870eb0615a4fcdccf5cb4d97653e6e878c993cea;hpb=875c717b027620ba3ebc16f35c0fccb2398e8ee4;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliLego.cxx b/STEER/AliLego.cxx index 870eb0615a4..5e5e53563ec 100644 --- a/STEER/AliLego.cxx +++ b/STEER/AliLego.cxx @@ -13,25 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -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$ */ ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// @@ -63,255 +45,415 @@ 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 "AliRun.h" -#include "AliConst.h" +#include "AliLegoGenerator.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), + fStopped(0) { - fHistRadl = 0; - fHistAbso = 0; - fHistGcm2 = 0; - fHistReta = 0; + // + // Default constructor + // } -//___________________________________________ -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), + fStopped(0) { -// specify the angular limits and the size of the rectangular box - - fGener = new AliLegoGenerator(ntheta, themin, themax, - nphi, phimin, phimax, rmin, rmax, zmax); - - gAlice->ResetGenerator(fGener); + // + // Copy constructor + // + 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), + 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", - 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 = AliDebugLevel(); } -//___________________________________________ -AliLego::~AliLego() +//_______________________________________________________________________ +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) { - delete fHistRadl; - delete fHistAbso; - delete fHistGcm2; - delete fHistReta; -} + // + // 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 = AliDebugLevel(); +} -//___________________________________________ -void AliLego::Run() +//_______________________________________________________________________ +AliLego::~AliLego() { - // loop on phi,theta bins - gMC->InitLego(); - Float_t thed, phid, eta; - for (Int_t i=0; i<=fGener->Nphi()*fGener->Ntheta(); ++i) { -// --- Set to 0 radiation length, absorption length and g/cm2 --- - fTotRadl = 0; - fTotAbso = 0; - fTotGcm2 = 0; - - gMC->ProcessEvent(); - - 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); - gAlice->FinishEvent(); - } - // store histograms in current Root file - fHistRadl->Write(); - fHistAbso->Write(); - fHistGcm2->Write(); - fHistReta->Write(); + // + // Destructor + // + delete fHistRadl; + delete fHistAbso; + delete fHistGcm2; + delete fGener; + delete fVolumesFwd; + delete fVolumesBwd; } -//___________________________________________ -void AliLego::StepManager() +//_______________________________________________________________________ +void AliLego::BeginEvent() { -// called from AliRun::Stepmanager from gustep. -// Accumulate the 3 parameters step by step + // + // --- Set to 0 radiation length, absorption length and g/cm2 --- + // + fTotRadl = 0; + fTotAbso = 0; + fTotGcm2 = 0; + fStopped = 0; - 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->TrackPosition(pos); - gMC->TrackMomentum(mom); - 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->TrackLength()) { - // 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 (fDebug) { + if (fErrorCondition) ToAliDebug(1, DumpVolumes()); + fVolumesFwd->Delete(); + fVolumesBwd->Delete(); + fStepsForward = 0; + fStepsBackward = 0; + fErrorCondition = 0; + if (gAlice->GetMCApp()->GetCurrentTrackNumber() == 0) fStepBack = 0; + } } -ClassImp(AliLegoGenerator) - -//___________________________________________ -AliLegoGenerator::AliLegoGenerator(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) : - AliGenerator(0), fRadMin(rmin), fRadMax(rmax), fZMax(zmax), fNtheta(ntheta), - fNphi(nphi), fThetaBin(ntheta), fPhiBin(-1), fCurTheta(0), fCurPhi(0) - +//_______________________________________________________________________ +void AliLego::FinishEvent() { - SetPhiRange(phimin,phimax); - SetThetaRange(themin,themax); - SetName("Lego"); + // + // Finish the event and update the histos + // + 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 AliLegoGenerator::Generate() +//_______________________________________________________________________ +void AliLego::FinishRun() { -// Create a geantino with kinematics corresponding to the current -// bins in theta and phi. - // - // Rootinos are 0 - const Int_t mpart = 0; - Float_t orig[3], pmom[3]; - Float_t t, cost, sint, cosp, sinp; - - // Prepare for next step - if(fThetaBin>=fNtheta-1) - if(fPhiBin>=fNphi-1) { - Warning("Generate","End of Lego Generation"); - return; - } else { - fPhiBin++; - printf("Generating rays in phi bin:%d\n",fPhiBin); - fThetaBin=0; - } else fThetaBin++; - - fCurTheta = (fThetaMin+(fThetaBin+0.5)*(fThetaMax-fThetaMin)/fNtheta); - fCurPhi = (fPhiMin+(fPhiBin+0.5)*(fPhiMax-fPhiMin)/fNphi); - cost = TMath::Cos(fCurTheta); - sint = TMath::Sin(fCurTheta); - cosp = TMath::Cos(fCurPhi); - sinp = TMath::Sin(fCurPhi); - - pmom[0] = cosp*sint; - pmom[1] = sinp*sint; - pmom[2] = cost; - - // --- Where to start - orig[0] = orig[1] = orig[2] = 0; - Float_t dalicz = 3000; - if (fRadMin > 0) { - t = PropagateCylinder(orig,pmom,fRadMin,dalicz); - orig[0] = pmom[0]*t; - orig[1] = pmom[1]*t; - orig[2] = pmom[2]*t; - if (TMath::Abs(orig[2]) > fZMax) return; - } - - Float_t polar[3]={0.,0.,0.}; - Int_t ntr; - gAlice->SetTrack(1, 0, mpart, pmom, orig, polar, 0, "LEGO ray", ntr); - + // Store histograms in current Root file + // + fHistRadl->Write(); + fHistAbso->Write(); + fHistGcm2->Write(); + + // Delete histograms from memory + fHistRadl->Delete(); fHistRadl=0; + fHistAbso->Delete(); fHistAbso=0; + fHistGcm2->Delete(); fHistGcm2=0; + // + if (fErrorCondition) ToAliError(DumpVolumes()); } -//___________________________________________ -Float_t AliLegoGenerator::PropagateCylinder(Float_t *x, Float_t *v, Float_t r, Float_t z) +//_______________________________________________________________________ +void AliLego::Copy(TObject&) const { -// Propagate to cylinder from inside - - Double_t hnorm, sz, t, t1, t2, t3, sr; - Double_t d[3]; - const Float_t kSmall = 1e-8; - const Float_t kSmall2 = kSmall*kSmall; - -// ---> Find intesection with Z planes - d[0] = v[0]; - d[1] = v[1]; - d[2] = v[2]; - hnorm = TMath::Sqrt(1/(d[0]*d[0]+d[1]*d[1]+d[2]*d[2])); - d[0] *= hnorm; - d[1] *= hnorm; - d[2] *= hnorm; - if (d[2] > kSmall) sz = (z-x[2])/d[2]; - else if (d[2] < -kSmall) sz = -(z+x[2])/d[2]; - else sz = 1.e10; // ---> Direction parallel to X-Y, no intersection + // + // Copy *this onto lego -- not implemented + // + AliFatal("Not implemented!"); +} -// ---> Intersection with cylinders -// Intersection point (x,y,z) -// (x,y,z) is on track : x=X(1)+t*D(1) -// y=X(2)+t*D(2) -// z=X(3)+t*D(3) -// (x,y,z) is on cylinder : x**2 + y**2 = R**2 -// -// (D(1)**2+D(2)**2)*t**2 -// +2.*(X(1)*D(1)+X(2)*D(2))*t -// +X(1)**2+X(2)**2-R**2=0 -// ---> Solve second degree equation - t1 = d[0]*d[0] + d[1]*d[1]; - if (t1 <= kSmall2) { - t = sz; // ---> Track parallel to the z-axis, take distance to planes - } else { - t2 = x[0]*d[0] + x[1]*d[1]; - t3 = x[0]*x[0] + x[1]*x[1]; - // ---> It should be positive, but there may be numerical problems - sr = (t2 +TMath::Sqrt(TMath::Max(t2*t2-(t3-r*r)*t1,0.)))/t1; - // ---> Find minimum distance between planes and cylinder - t = TMath::Min(sz,sr); - } - return t; +//_______________________________________________________________________ +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; + 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(); + + 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 && 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 && 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 && !(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"); +}