Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages
CalDetTracker

CDSimpleMC.cxx

Go to the documentation of this file.
00001 
00002 
00003 // Program name: CDSimpleMC.cxx
00004 //                                                                     
00005 // Package: CalDetTracker   
00006 //                                                                    
00007 // Coded by Jeff Hartnell Jul-Dec/2003                             
00008 //                                                                    
00009 // Purpose: To make plots and analyse the Tracker*.root files
00010 //                                                                    
00011 // Contact: jeffrey.hartnell@physics.ox.ac.uk                         
00013 
00014 #include <math.h>
00015 
00016 #include "MessageService/MsgService.h"
00017 
00018 #include "CalDetTracker/CDSimpleMC.h"
00019 
00020 using namespace std;
00021 
00022 ClassImp(CDSimpleMC)
00023 
00024 CVSID("$Id: CDSimpleMC.cxx,v 1.8 2006/05/18 21:52:56 admarino Exp $");
00025 
00026 //......................................................................
00027 
00028 CDSimpleMC::CDSimpleMC()
00029 {
00030   MSG("CDSimpleMC",Msg::kDebug)
00031     <<"Running CDSimpleMC Constructor..."<<endl;
00032 
00033   fS="";//general purpose string
00034   fParticleMomentum=1800;//set a default momentum
00035   fMuMass=105.658369;
00036   fSlices=1000;
00037 
00038   //set the default thicknesses
00039   fThkFe=2.50;//metric iron at CalDet
00040   fThkAl=0.0508;//module skin thickness=0.508 mm (0.020 inches)
00041   fThkSc=0.95;//(1.0-2*fThkScTiO2) in TDR and GMINOS
00042   fThkScTiO2=0.025;//0.025 in TDR and GMINOS
00043   fThkAir=2.4;//Not completely sure of this one, but it doesn't matter
00044   this->CalcPlaneThk();
00045 
00046   //set the default results
00047   fTotalScEnLoss=0;
00048   fTotalNumScPlanesHit=0;
00049 
00050   //set the default process to use for dE_dx
00051   fELossProcess=Process::eIonization;
00052   //case eIonization   
00053   //case ePairProduction 
00054   //case eBremsstrahlung 
00055   //case eNuclearInteraction
00056   //case eAll     
00057 
00058   //set the standard materials
00059   fStdMatFe=MinosMaterial::eIronCD;//CalDet iron
00060   fStdMatAl=MinosMaterial::eAluminium;
00061   fStdMatSc=MinosMaterial::ePolystyreneMinos;
00062   fStdMatScTiO2=MinosMaterial::ePolystyreneMinos;
00063 
00064   //set the densities
00065   fRoFe=MinosMaterial::Density(fStdMatFe);//Was 7.874;
00066   fRoAl=MinosMaterial::Density(fStdMatAl);//Was 2.70;
00067   fRoScTiO2=MinosMaterial::Density(fStdMatScTiO2);//Was 1.06;
00068   fRoSc=MinosMaterial::Density(fStdMatSc);//Was 1.06;
00069 
00070   //get new materials
00071   fFe=new Material(fStdMatFe);
00072   fAl=new Material(fStdMatAl);
00073   fSc=new Material(fStdMatSc);
00074   
00075   //get new energyLoss objects
00076   fIonFe=new MuEnergyLoss(*fFe,fModelSelector);
00077   fIonAl=new MuEnergyLoss(*fAl,fModelSelector);
00078   fIonSc=new MuEnergyLoss(*fSc,fModelSelector);
00079   fIonScTiO2=new MuEnergyLoss(*fSc,fModelSelector);//assume Scint
00080 
00081   MSG("CDSimpleMC",Msg::kDebug)
00082     <<"Finished CDSimpleMC Constructor"<<endl;
00083 }
00084 
00085 //......................................................................
00086 
00087 CDSimpleMC::~CDSimpleMC()
00088 {
00089   MSG("CDSimpleMC",Msg::kDebug)
00090       <<"Running CDSimpleMC Destructor..."<<endl;
00091 
00092   //clean up...
00093 
00094   //standard materials
00095   if (fFe) delete fFe;
00096   if (fAl) delete fAl;
00097   if (fSc) delete fSc;
00098   
00099   //energyLoss objects
00100   if (fIonFe) delete fIonFe;
00101   if (fIonAl) delete fIonAl;
00102   if (fIonSc) delete fIonSc;
00103   if (fIonScTiO2) delete fIonScTiO2;
00104 }
00105 
00106 //......................................................................
00107 
00108 void CDSimpleMC::SetThkFe(Double_t thkFe)
00109 {
00110   fThkFe=thkFe;
00111   this->CalcPlaneThk();
00112 }
00113 
00114 //......................................................................
00115 
00116 void CDSimpleMC::CalcPlaneThk()
00117 {
00118   fThkPlane=fThkAl+fThkScTiO2+fThkSc+fThkScTiO2+fThkAl+fThkFe+fThkAir;
00119 }
00120 
00121 //......................................................................
00122 
00123 Double_t CDSimpleMC::GetParticleEnergy() const
00124 {
00125   //convert momentum to energy
00126   return sqrt(fParticleMomentum*fParticleMomentum+fMuMass*fMuMass);
00127 }
00128 
00129 //......................................................................
00130 
00131 void CDSimpleMC::PrintDensities() const
00132 {
00133   MSG("CDSimpleMC",Msg::kInfo)
00134     <<"Densities: Fe="<<fRoFe
00135     <<", Al="<<fRoAl
00136     <<", Sc="<<fRoSc
00137     <<", ScTiO2="<<fRoScTiO2<<endl;
00138 }
00139 
00140 //......................................................................
00141 
00142 void CDSimpleMC::PrintThicknesses() const
00143 {
00144   MSG("CDSimpleMC",Msg::kInfo)
00145     <<"Thicknesses: Fe="<<fThkFe
00146     <<", Al="<<fThkAl
00147     <<", Sc="<<fThkSc
00148     <<", ScTiO2="<<fThkScTiO2
00149     <<", Air="<<fThkAir
00150     <<", Plane="<<fThkPlane
00151     <<endl;
00152 }
00153 
00154 //......................................................................
00155 
00156 Double_t CDSimpleMC::dE_dxIterative(MuEnergyLoss& mat,
00157                                     Double_t energyRemaining,
00158                                     Double_t totalThk,
00159                                     Double_t rho,
00160                                     Double_t energyCutOff,
00161                                     Int_t slices,
00162                                     Double_t& fractionThrough) const
00163 {
00166 
00167   Double_t totalEnergyDeposited=0;
00168   Double_t sliceThk=totalThk/slices;
00169 
00170   //set the fractionThrough to 0 as default
00171   //then you can always add the fraction through to the plane
00172   fractionThrough=0;
00173 
00174   for (Int_t s=0;s<slices;s++){
00175 
00176     if (energyRemaining<energyCutOff){
00177       //calc fractionThrough material
00178       fractionThrough=s*sliceThk/totalThk;
00179 
00180       MSG("CDSimpleMC",Msg::kDebug) 
00181         <<"All energy lost at slice "<<s+1<<"/"<<slices
00182         <<" ("<<100.*fractionThrough<<"% through material)"<<endl;
00183       return -1;
00184     }
00185 
00186     //calculate energy deposited in the slice
00187     Double_t energyDeposited=mat.dE_dx
00188       (energyRemaining,fELossProcess)*rho*sliceThk;
00189     totalEnergyDeposited+=energyDeposited;
00190 
00191     //increment the energy remaining
00192     energyRemaining-=energyDeposited;
00193   }
00194   
00195   return totalEnergyDeposited;
00196 }
00197 
00198 //......................................................................
00199 
00200 void CDSimpleMC::CalcRatioScToFe()
00201 {
00202   fEnRatio.clear();
00203   fXRatio.clear();
00204   fXPlanesRatio.clear();
00205 
00206   //get beam energy
00207   Float_t energyRemaining=this->GetParticleEnergy();
00208   Float_t distTravelled=0;
00209   Double_t fractionThrough=0;
00210 
00211   while (energyRemaining>fMuMass){
00212     
00213     MSG("CDAnalysis",Msg::kDebug) 
00214       <<"CalcRatioScToFe: energyRemaining="<<energyRemaining<<endl;
00215     
00217     //calculate the ratio of energy deposited in Fe to Sc
00219     Float_t eFe=dE_dxIterative(*fIonFe,energyRemaining,
00220                                1,fRoFe,fMuMass,fSlices,fractionThrough);
00221     Float_t eSc=dE_dxIterative(*fIonSc,energyRemaining,
00222                                1,fRoSc,fMuMass,fSlices,fractionThrough);
00223 
00224     //jump out of loop if particle has run out of energy
00225     if (energyRemaining-eFe-eSc<fMuMass || 
00226         eFe==-1 || eSc==-1) break;
00227 
00228     fEnRatio.push_back(eFe/eSc);
00229     fXRatio.push_back(distTravelled);
00230     fXPlanesRatio.push_back(distTravelled/fThkPlane);
00231 
00232     energyRemaining-=eFe;
00233     energyRemaining-=eSc;
00234     distTravelled+=fThkPlane;
00235   }
00236 }
00237 
00238 //......................................................................
00239 
00240 void CDSimpleMC::RunMC()
00241 {
00242   //clear all the vectors so this can be used in a loop
00243   fEn.clear();
00244   fMom.clear();
00245   fdE_dxFe.clear();
00246   fdE_dxAl.clear();
00247   fdE_dxSc.clear();
00248 
00249   fEnDeposited.clear();
00250   fX.clear();
00251   fXPlanes.clear();
00252 
00253   fEnFeDeposited.clear();
00254   fXFe.clear();
00255   fXFePlanes.clear();
00256 
00257   fEnAlDeposited.clear();
00258   fXAl.clear();
00259   fXAlPlanes.clear();
00260 
00261   fEnScTiO2Deposited.clear();
00262   fXScTiO2.clear();
00263   fXScTiO2Planes.clear();
00264 
00265   fEnNotScDeposited.clear();
00266   fXNotSc.clear();
00267   fXNotScPlanes.clear();
00268 
00269   fEnScDeposited.clear();
00270   fEnSc.clear();
00271   fMomSc.clear();
00272   fXSc.clear();
00273   fXScPlanes.clear();
00274 
00275   fEnDepositedPl.clear();
00276 
00277   //print out useful info
00278   MSG("CDSimpleMC",Msg::kInfo)
00279     <<"This run of the simple MC is configured with the following:"
00280     <<endl
00281     <<"Initial particle momentum="<<fParticleMomentum<<endl;
00282   this->PrintDensities();
00283   this->PrintThicknesses();
00284 
00285   //get beam energy
00286   Float_t energyRemaining=this->GetParticleEnergy();
00287   Float_t distTravelled=0.001;//cm
00288   Float_t sumScEnLoss=0;
00289   Int_t plCounter=0;
00290   Double_t fractionThrough=0;
00291   Double_t fractionThroughFe=0;
00292 
00293   //calcuate the dE/dx
00294   for (Float_t energy=107;energy<energyRemaining;energy+=30){
00295     fEn.push_back(energy);
00296     Double_t mom=sqrt(pow(energy,2)-pow(fMuMass,2));
00297     MAXMSG("CDSimpleMC",Msg::kDebug,100) 
00298       <<"energy="<<energy<<", mom="<<mom<<endl;
00299     fMom.push_back(mom);
00300     fMomGeV.push_back(mom/1000);
00301     fdE_dxFe.push_back(fIonFe->dE_dx(energy,fELossProcess));
00302     fdE_dxAl.push_back(fIonAl->dE_dx(energy,fELossProcess));
00303     fdE_dxSc.push_back(fIonSc->dE_dx(energy,fELossProcess));
00304   }
00305 
00306   //loop over the particle energy until it has no kinetic energy left
00307   while (energyRemaining>fMuMass){
00308 
00309     MSG("CDSimpleMC",Msg::kDebug) 
00310       <<"energyRemaining="<<energyRemaining<<endl;
00311 
00313     //calculate the energy deposited
00314     //iterate energy remaining at each stage
00316 
00317     //variables to store the sum of the energy dep in aluminium
00318     Double_t alx2En=0;
00319     Double_t alx2X=0;
00320     Double_t tiO2x2En=0;
00321     Double_t tiO2x2X=0;
00322 
00323     //***************************
00324     //calc energy deposited in 1st Al module skin
00325     //***************************
00326     Double_t energyDeposited=dE_dxIterative(*fIonAl,energyRemaining,
00327                                             fThkAl,fRoAl,fMuMass,fSlices,
00328                                             fractionThrough);
00329     MSG("CDSimpleMC",Msg::kDebug) 
00330       <<"1st Al ELoss: "<<energyDeposited
00331       <<" ("<<energyRemaining<<")"<<endl;
00332 
00333     if (energyRemaining-energyDeposited<fMuMass || 
00334         energyDeposited==-1) break;
00335 
00336     //get the first al energy deposit
00337     alx2En+=energyDeposited;
00338     alx2X+=distTravelled;
00339 
00340     fEnDeposited.push_back(energyDeposited);
00341     fX.push_back(distTravelled);
00342     fXPlanes.push_back(distTravelled/fThkPlane);
00343     fEnAlDeposited.push_back(energyDeposited);
00344     fXAl.push_back(distTravelled);
00345     fXAlPlanes.push_back(distTravelled/fThkPlane);
00346     fEnDepositedPl.push_back(energyDeposited);
00347     vector<Double_t>::iterator enDepositedPlIter=fEnDepositedPl.end()-1;
00348 
00349     //increment the distance travelled
00350     distTravelled+=fThkAl;
00351     Double_t plEnergyDeposited=energyDeposited;
00352     energyRemaining-=energyDeposited;
00353 
00354     //***************************
00355     //calc energy deposited in 1st TiO2 reflective layer
00356     //***************************
00357     energyDeposited=dE_dxIterative(*fIonScTiO2,energyRemaining,
00358                                    fThkScTiO2,fRoScTiO2,fMuMass,fSlices,
00359                                    fractionThrough);
00360     MSG("CDSimpleMC",Msg::kDebug) 
00361       <<"1st ScTiO2 ELoss: "<<energyDeposited
00362       <<" ("<<energyRemaining<<")"<<endl;
00363 
00364     if (energyRemaining-energyDeposited<fMuMass || 
00365         energyDeposited==-1) break;
00366 
00367     //get the first TiO2 energy deposit
00368     tiO2x2En+=energyDeposited;
00369     tiO2x2X+=distTravelled;
00370 
00371     fEnDeposited.push_back(energyDeposited);
00372     fX.push_back(distTravelled);
00373     fXPlanes.push_back(distTravelled/fThkPlane);
00374     fEnScTiO2Deposited.push_back(energyDeposited);
00375     fXScTiO2.push_back(distTravelled);
00376     fXScTiO2Planes.push_back(distTravelled/fThkPlane);
00377     *enDepositedPlIter+=energyDeposited;
00378 
00379     //increment the distance travelled
00380     distTravelled+=fThkScTiO2;
00381     plEnergyDeposited+=energyDeposited;
00382     energyRemaining-=energyDeposited;
00383 
00384     //***************************
00385     //calc energy deposited in Sc
00386     //***************************
00387     energyDeposited=dE_dxIterative(*fIonSc,energyRemaining,
00388                                    fThkSc,fRoSc,fMuMass,fSlices,
00389                                    fractionThrough);
00390     MSG("CDSimpleMC",Msg::kDebug) 
00391       <<"Sc ELoss: "<<energyDeposited<<" ("<<energyRemaining<<")"<<endl;
00392 
00393     //count the planes here since it is active planes that matter
00394     //the particle could have only just penetrated the Scint but
00395     //that will probably be rare
00396     plCounter++;
00397 
00398     if (energyRemaining-energyDeposited<fMuMass || 
00399         energyDeposited==-1) break;
00400 
00401     fEnDeposited.push_back(energyDeposited);
00402     fX.push_back(distTravelled);
00403     fXPlanes.push_back(distTravelled/fThkPlane);
00404     fEnScDeposited.push_back(energyDeposited);
00405     fEnSc.push_back(energyRemaining);
00406     Double_t momemtum=sqrt(pow(energyRemaining,2)-pow(fMuMass,2));
00407     fMomSc.push_back(momemtum);
00408     fMomScGeV.push_back(momemtum/1000);
00409     fXSc.push_back(distTravelled);
00410     fXScPlanes.push_back(distTravelled/fThkPlane);
00411     *enDepositedPlIter+=energyDeposited;
00412 
00413     //increment the distance travelled
00414     distTravelled+=fThkSc;
00415     //plEnergyDeposited+=energyDeposited;//not Sc
00416     energyRemaining-=energyDeposited;
00417     //keep running total of energy lost in Scint
00418     sumScEnLoss+=energyDeposited;
00419 
00420     //***************************
00421     //calc energy deposited in 2nd TiO2 reflective layer
00422     //***************************
00423     energyDeposited=dE_dxIterative(*fIonScTiO2,energyRemaining,
00424                                    fThkScTiO2,fRoScTiO2,fMuMass,fSlices,
00425                                    fractionThrough);
00426     MSG("CDSimpleMC",Msg::kDebug) 
00427       <<"2nd ScTiO2 ELoss: "<<energyDeposited
00428       <<" ("<<energyRemaining<<")"<<endl;
00429 
00430     if (energyRemaining-energyDeposited<fMuMass || 
00431         energyDeposited==-1) break;
00432 
00433     //get the second TiO2 energy deposit
00434     tiO2x2En+=energyDeposited;
00435     tiO2x2X+=distTravelled;
00436     tiO2x2X/=2;//divide the distance by 2 to get the average
00437     fEnTiO2x2Deposited.push_back(tiO2x2En);
00438     fXTiO2x2.push_back(tiO2x2X);
00439 
00440     fEnDeposited.push_back(energyDeposited);
00441     fX.push_back(distTravelled);
00442     fXPlanes.push_back(distTravelled/fThkPlane);
00443     fEnScTiO2Deposited.push_back(energyDeposited);
00444     fXScTiO2.push_back(distTravelled);
00445     fXScTiO2Planes.push_back(distTravelled/fThkPlane);
00446     *enDepositedPlIter+=energyDeposited;
00447 
00448     //increment the distance travelled
00449     distTravelled+=fThkScTiO2;
00450     plEnergyDeposited+=energyDeposited;
00451     energyRemaining-=energyDeposited;
00452 
00453     //***************************
00454     //calc energy deposited in 2nd Al module skin
00455     //***************************
00456     energyDeposited=dE_dxIterative(*fIonAl,energyRemaining,
00457                                    fThkAl,fRoAl,fMuMass,fSlices,
00458                                    fractionThrough);
00459     MSG("CDSimpleMC",Msg::kDebug) 
00460       <<"2nd Al ELoss: "<<energyDeposited
00461       <<" ("<<energyRemaining<<")"<<endl;
00462 
00463     if (energyRemaining-energyDeposited<fMuMass || 
00464         energyDeposited==-1) break;
00465 
00466     //get the second al energy deposit
00467     alx2En+=energyDeposited;
00468     alx2X+=distTravelled;
00469     alx2X/=2;//divide the distance by 2 to get the average
00470     fEnAlx2Deposited.push_back(alx2En);
00471     fXAlx2.push_back(alx2X);
00472 
00473     fEnDeposited.push_back(energyDeposited);
00474     fX.push_back(distTravelled);
00475     fXPlanes.push_back(distTravelled/fThkPlane);
00476     fEnAlDeposited.push_back(energyDeposited);
00477     fXAl.push_back(distTravelled);
00478     fXAlPlanes.push_back(distTravelled/fThkPlane);
00479     *enDepositedPlIter+=energyDeposited;
00480 
00481     //increment the distance travelled
00482     distTravelled+=fThkAl;
00483     plEnergyDeposited+=energyDeposited;
00484     energyRemaining-=energyDeposited;
00485 
00486     //***************************
00487     //calc energy deposited in Fe
00488     //***************************
00489     energyDeposited=dE_dxIterative(*fIonFe,energyRemaining,
00490                                    fThkFe,fRoFe,fMuMass,fSlices,
00491                                    fractionThrough);
00492     MSG("CDSimpleMC",Msg::kDebug) 
00493       <<"Fe ELoss: "<<energyDeposited<<" ("<<energyRemaining<<")"<<endl;
00494 
00495     //set the fractionThrough Fe variable
00496     fractionThroughFe=fractionThrough;
00497 
00498     if (energyRemaining-energyDeposited<fMuMass || 
00499         energyDeposited==-1) break;
00500 
00501     fEnDeposited.push_back(energyDeposited);
00502     fX.push_back(distTravelled);
00503     fXPlanes.push_back(distTravelled/fThkPlane);
00504     fEnFeDeposited.push_back(energyDeposited);
00505     fXFe.push_back(distTravelled);
00506     fXFePlanes.push_back(distTravelled/fThkPlane);
00507     *enDepositedPlIter+=energyDeposited;//sum energy in plane
00508 
00509     //increment the distance travelled and energy remaining
00510     distTravelled+=fThkFe;
00511     plEnergyDeposited+=energyDeposited;
00512     energyRemaining-=energyDeposited;
00513 
00514     //calc energy deposited in Air
00515     energyDeposited=0;
00516     //eDeposited.push_back(energyDeposited);
00517     //x.push_back(distTravelled);
00518     //xPlanes.push_back(distTravelled/fThkPlane);
00519     //*eDepositedPlIter+=energyDeposited;//sum energy in plane
00520 
00521     //increment the distance travelled
00522     distTravelled+=fThkAir;
00523     energyRemaining-=energyDeposited;
00524 
00525     //add the not scint energy to its vector
00526     fEnNotScDeposited.push_back(plEnergyDeposited);
00527     fXNotSc.push_back(distTravelled);
00528     fXNotScPlanes.push_back(distTravelled/fThkPlane-1);//1st time=pl 0
00529   }
00530 
00531   //set the data members
00532   fTotalScEnLoss=sumScEnLoss;
00533   //add on the fractionThrough Fe to the last plane
00534   //this is not 100% right but close enough since Sc/Fe=2/30
00535   fTotalNumScPlanesHit=plCounter+fractionThroughFe;
00536 
00537 }
00538 
00539 //......................................................................
Generated on Mon Feb 15 11:06:30 2010 for loon by  doxygen 1.3.9.1