X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=STEER%2FAliLego.cxx;h=4a4bff9bf6a1e355e47461a29484e09c3e63d9e6;hb=ef278eae2ab8c4df7a67b5488fb05ae29f9e3ccb;hp=003856342ca41552e78ce1db66df2383ab72b240;hpb=00719c1b1e1da0adeb2052d41395e01cd224e80d;p=u%2Fmrichter%2FAliRoot.git diff --git a/STEER/AliLego.cxx b/STEER/AliLego.cxx index 003856342ca..4a4bff9bf6a 100644 --- a/STEER/AliLego.cxx +++ b/STEER/AliLego.cxx @@ -13,15 +13,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -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$ */ ////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////// @@ -53,235 +45,400 @@ Introduction of the Copyright and cvs Log // ////////////////////////////////////////////////////////////// +#include "TClonesArray.h" +#include "TH2.h" #include "TMath.h" +#include "TString.h" +#include "TVirtualMC.h" + +#include "AliConst.h" +#include "AliDebugVolume.h" #include "AliLego.h" +#include "AliLegoGenerator.h" #include "AliRun.h" -#include "AliConst.h" +#include "AliMC.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) { - fHistRadl = 0; - fHistAbso = 0; - fHistGcm2 = 0; - fHistReta = 0; + // + // Default constructor + // } -//___________________________________________ -AliLego::AliLego(const char *name, const char *title) - : TNamed(name,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) { - fHistRadl = 0; - fHistAbso = 0; - fHistGcm2 = 0; - fHistReta = 0; + // + // Copy constructor + // + lego.Copy(*this); } -//___________________________________________ -AliLego::~AliLego() + +//_______________________________________________________________________ +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) { - delete fHistRadl; - delete fHistAbso; - delete fHistGcm2; - delete fHistReta; + // + // 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 = gAlice->GetDebug(); } -//___________________________________________ -void AliLego::GenerateKinematics() +//_______________________________________________________________________ +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) { -// 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; - -// --- Set to 0 radiation length, absorption length and g/cm2 --- - fTotRadl = 0; - fTotAbso = 0; - fTotGcm2 = 0; - - fCurTheta = (fThetaMin+(fThetaBin-0.5)*(fThetaMax-fThetaMin)/fNtheta)*kDegrad; - fCurPhi = (fPhiMin+(fPhiBin-0.5)*(fPhiMax-fPhiMin)/fNphi)*kDegrad; - 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; - } + // 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); + // + // -// --- We do start here - Float_t polar[3]={0.,0.,0.}; - Int_t ntr; - gAlice->SetTrack(1, 0, mpart, pmom, orig, polar, 0, "LEGO ray", ntr); + fVolumesFwd = new TClonesArray("AliDebugVolume",1000); + fVolumesBwd = new TClonesArray("AliDebugVolume",1000); + fDebug = gAlice->GetDebug(); } -//___________________________________________ -void AliLego::Init(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) +//_______________________________________________________________________ +AliLego::~AliLego() { -// specify the angular limits and the size of the rectangular box - fNtheta = ntheta; - fThetaMin = themin; - fThetaMax = themax; - fNphi = nphi; - fPhiMin = phimin; - fPhiMax = phimax; - fRadMin = rmin; - fRadMax = rmax; - fZMax = zmax; - Float_t etamin = -TMath::Log(TMath::Tan(TMath::Min((Double_t)fThetaMax*kDegrad/2,TMath::Pi()/2-1.e-10))); - Float_t etamax = -TMath::Log(TMath::Tan(TMath::Max((Double_t)fThetaMin*kDegrad/2, 1.e-10))); + // + // Destructor + // + delete fHistRadl; + delete fHistAbso; + delete fHistGcm2; + delete fGener; + delete fVolumesFwd; + delete fVolumesBwd; +} - fHistRadl = new TH2F("hradl","Radiation length map", nphi,phimin,phimax,ntheta,themin,themax); - fHistAbso = new TH2F("habso","Interaction length map", nphi,phimin,phimax,ntheta,themin,themax); - 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); +//_______________________________________________________________________ +void AliLego::BeginEvent() +{ + // + // --- Set to 0 radiation length, absorption length and g/cm2 --- + // + 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; + } } -//___________________________________________ -Float_t AliLego::PropagateCylinder(Float_t *x, Float_t *v, Float_t r, Float_t z) +//_______________________________________________________________________ +void AliLego::FinishEvent() { -// 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; + // + // 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); +} -// ---> 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 +//_______________________________________________________________________ +void AliLego::FinishRun() +{ + // + // Store histograms in current Root file + // + fHistRadl->Write(); + fHistAbso->Write(); + fHistGcm2->Write(); -// ---> 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; + // Delete histograms from memory + fHistRadl->Delete(); fHistRadl=0; + fHistAbso->Delete(); fHistAbso=0; + fHistGcm2->Delete(); fHistGcm2=0; + // + if (fErrorCondition) DumpVolumes(); } -//___________________________________________ -void AliLego::Run() +//_______________________________________________________________________ +void AliLego::Copy(TObject&) const { - // loop on phi,theta bins - gMC->InitLego(); - Float_t thed, phid, eta; - for (fPhiBin=1; fPhiBin<=fNphi; fPhiBin++) { - printf("AliLego Generating rays in phi bin:%d\n",fPhiBin); - for (fThetaBin=1; fThetaBin<=fNtheta; fThetaBin++) { - gMC->Gtrigi(); - gMC->Gtrigc(); - GenerateKinematics(); - gMC->Gtreve_root(); - - thed = fCurTheta*kRaddeg; - phid = fCurPhi*kRaddeg; - eta = -TMath::Log(TMath::Tan(TMath::Max( - TMath::Min((Double_t)fCurTheta/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(); + // + // Copy *this onto lego -- not implemented + // + 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]; - 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]) > fZMax || - pos[0]*pos[0] +pos[1]*pos[1] > fRadMax*fRadMax) { - 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 = PropagateCylinder(vect,dir,fRadMax,fZMax); - - fTotAbso += step/absl; - fTotRadl += step/radl; - fTotGcm2 += step*dens; + 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 + 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()); + + 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"); +}