CMSSW Reference Manual

typedef std::pair<XYZPoint,double> HDShower::Spot

Definition at line 29 of file HDShower.h.

typedef std::pair<unsigned int, double> HDShower::Step

Definition at line 30 of file HDShower.h.

typedef Steps::const_iterator HDShower::step_iterator

Definition at line 32 of file HDShower.h.

typedef std::vector<Step> HDShower::Steps

Definition at line 31 of file HDShower.h.

typedef math::XYZVector HDShower::XYZPoint

Definition at line 27 of file HDShower.h.

HDShower::HDShower	(	const RandomEngine *	engine,
		HDShowerParametrization *	myParam,
		EcalHitMaker *	myGrid,
		HcalHitMaker *	myHcalHitMaker,
		int	onECAL,
		double	epart
	)

Definition at line 30 of file HDShower.cc.

00036   : theParam(myParam), 
00037     theGrid(myGrid),
00038     theHcalHitMaker(myHcalHitMaker),
00039     onEcal(onECAL),
00040     e(epart),
00041     random(engine)
00042 { 
00043   // To get an access to constants read in FASTCalorimeter
00044   //  FASTCalorimeter * myCalorimeter= FASTCalorimeter::instance();
00045 
00046   // Values taken from FamosGeneric/FamosCalorimeter/src/FASTCalorimeter.cc
00047   lossesOpt       = myParam->hsParameters()->getHDlossesOpt();
00048   nDepthSteps     = myParam->hsParameters()->getHDnDepthSteps();
00049   nTRsteps        = myParam->hsParameters()->getHDnTRsteps();
00050   transParam      = myParam->hsParameters()->getHDtransParam();
00051   eSpotSize       = myParam->hsParameters()->getHDeSpotSize();
00052   depthStep       = myParam->hsParameters()->getHDdepthStep();
00053   criticalEnergy  = myParam->hsParameters()->getHDcriticalEnergy();
00054   maxTRfactor     = myParam->hsParameters()->getHDmaxTRfactor();
00055   balanceEH       = myParam->hsParameters()->getHDbalanceEH();
00056   hcalDepthFactor = myParam->hsParameters()->getHDhcalDepthFactor();
00057 
00058   // Special tr.size fluctuations 
00059   transParam *= (1. + random->flatShoot()); 
00060 
00061   // Special long. fluctuations
00062   hcalDepthFactor +=  0.05 * (2.* random->flatShoot() - 1.);
00063 
00064   transFactor = 1.;   // normally 1, in HF - might be smaller 
00065                       // to take into account
00066                       // a narrowness of the HF shower (Cherenkov light) 
00067 
00068   // simple protection ...
00069   if(e < 0) e = 0.;
00070 
00071   // Get the Famos Histos pointer
00072   //  myHistos = FamosHistos::instance();
00073   //  std::cout << " Hello FamosShower " << std::endl;
00074   
00075   theECALproperties = theParam->ecalProperties();
00076   theHCALproperties = theParam->hcalProperties();
00077 
00078   double emax = theParam->emax(); 
00079   double emid = theParam->emid(); 
00080   double emin = theParam->emin(); 
00081   double effective = e;
00082     
00083 
00084   if( e < emid ) {
00085     theParam->setCase(1);
00086     // avoid "underflow" below Emin (for parameters calculation only)
00087     if(e < emin) effective = emin; 
00088   } 
00089   else 
00090     theParam->setCase(2);
00091  
00092   // Avoid "overflow" beyond Emax (for parameters)
00093   if(effective > 0.5 * emax) {
00094     eSpotSize *= 2.5;
00095     if(effective > emax) {
00096       effective = emax; 
00097       eSpotSize *= 2.; 
00098       depthStep *= 2.;
00099     }
00100   }
00101 
00102 
00103   if(debug == 2 )
00104     LogDebug("FastCalorimetry") <<  " HDShower : " << std::endl 
00105          << "       Energy   "  <<               e << std::endl     
00106          << "      lossesOpt "  <<       lossesOpt << std::endl  
00107          << "    nDepthSteps "  <<     nDepthSteps << std::endl
00108          << "       nTRsteps "  <<        nTRsteps << std::endl
00109          << "     transParam "  <<      transParam << std::endl
00110          << "      eSpotSize "  <<       eSpotSize << std::endl
00111          << " criticalEnergy "  <<  criticalEnergy << std::endl
00112          << "    maxTRfactor "  <<     maxTRfactor << std::endl
00113          << "      balanceEH "  <<       balanceEH << std::endl
00114          << "hcalDepthFactor "  << hcalDepthFactor << std::endl;
00115 
00116 
00117   double alpEM1 = theParam->alpe1();
00118   double alpEM2 = theParam->alpe2();
00119 
00120   double betEM1 = theParam->bete1();
00121   double betEM2 = theParam->bete2();
00122 
00123   double alpHD1 = theParam->alph1();
00124   double alpHD2 = theParam->alph2();
00125 
00126   double betHD1 = theParam->beth1();
00127   double betHD2 = theParam->beth2();
00128 
00129   double part1 = theParam->part1();
00130   double part2 = theParam->part2();
00131 
00132   aloge = std::log(effective);
00133  
00134   double edpar = (theParam->e1() + aloge * theParam->e2()) * effective;
00135   double aedep = std::log(edpar);
00136 
00137   if(debug == 2)
00138     LogDebug("FastCalorimetry") << " HDShower : " << std::endl
00139          << "     edpar " <<   edpar << "   aedep " << aedep << std::endl 
00140          << "    alpEM1 " <<  alpEM1 << std::endl  
00141          << "    alpEM2 " <<  alpEM2 << std::endl  
00142          << "    betEM1 " <<  betEM1 << std::endl  
00143          << "    betEM2 " <<  betEM2 << std::endl  
00144          << "    alpHD1 " <<  alpHD1 << std::endl  
00145          << "    alpHD2 " <<  alpHD2 << std::endl  
00146          << "    betHD1 " <<  betHD1 << std::endl  
00147          << "    betHD2 " <<  betHD2 << std::endl  
00148          << "     part1 " <<   part1 << std::endl  
00149          << "     part2 " <<   part2 << std::endl; 
00150 
00151   // private members to set
00152   theR1  = theParam->r1();
00153   theR2  = theParam->r2();
00154   theR3  = theParam->r3();
00155 
00156   alpEM  = alpEM1 + alpEM2 * aedep;
00157   tgamEM = tgamma(alpEM);
00158   betEM  = betEM1 - betEM2 * aedep;
00159   alpHD  = alpHD1 + alpHD2 * aedep;
00160   tgamHD = tgamma(alpHD);
00161   betHD  = betHD1 - betHD2 * aedep;
00162   part   = part1  -  part2 * aedep;
00163   if(part > 1.) part = 1.;          // protection - just in case of 
00164 
00165   if(debug  == 2 )
00166     LogDebug("FastCalorimetry") << " HDShower : " << std::endl 
00167          << "    alpEM " <<  alpEM << std::endl  
00168          << "   tgamEM " << tgamEM << std::endl
00169          << "    betEM " <<  betEM << std::endl
00170          << "    alpHD " <<  alpHD << std::endl
00171          << "   tgamHD " << tgamHD << std::endl
00172          << "    betHD " <<  betHD << std::endl
00173          << "     part " <<   part << std::endl;
00174        
00175 
00176   if(onECAL){
00177     lambdaEM = theParam->ecalProperties()->interactionLength();
00178     x0EM     = theParam->ecalProperties()->radLenIncm();
00179   } 
00180   else {
00181     lambdaEM = 0.;
00182     x0EM     = 0.;
00183   }
00184   lambdaHD = theParam->hcalProperties()->interactionLength();
00185   x0HD     = theParam->hcalProperties()->radLenIncm();
00186 
00187   if(debug == 2)
00188     LogDebug("FastCalorimetry") << " HDShower e " << e        << std::endl
00189          << "          x0EM = " << x0EM     << std::endl 
00190          << "          x0HD = " << x0HD     << std::endl 
00191          << "         lamEM = " << lambdaEM << std::endl
00192          << "         lamHD = " << lambdaHD << std::endl;
00193 
00194 
00195   // Starting point of the shower
00196   // try first with ECAL lambda  
00197 
00198   double sum1   = 0.;  // lambda path from the ECAL/HF entrance;  
00199   double sum2   = 0.;  // lambda path from the interaction point;
00200   double sum3   = 0.;  // x0     path from the interaction point;
00201   int  nsteps   = 0;   // full number of longitudinal steps (counter);
00202 
00203   int nmoresteps;      // how many longitudinal steps in addition to 
00204                        // one (if interaction happens there) in ECAL
00205 
00206   mip = 1;             // just to initiate particle as MIP in ECAL    
00207 
00208   if(e < criticalEnergy ) nmoresteps = 1;  
00209   else                    nmoresteps = nDepthSteps;
00210   
00211   double depthECAL  = 0.;
00212   double depthGAP   = 0.;
00213   double depthGAPx0 = 0.; 
00214   if(onECAL ) {
00215     depthECAL  = theGrid->ecalTotalL0();         // ECAL depth segment
00216     depthGAP   = theGrid->ecalHcalGapTotalL0();  // GAP  depth segment
00217     depthGAPx0 = theGrid->ecalHcalGapTotalX0();  // GAP  depth x0
00218   }
00219   
00220   double depthHCAL   = theGrid->hcalTotalL0();    // HCAL depth segment
00221   double depthToHCAL = depthECAL + depthGAP;    
00222 
00223   //---------------------------------------------------------------------------
00224   // Depth simulation & various protections, among them
00225   // if too deep - get flat random in the allowed region
00226   // if no HCAL material behind - force to deposit in ECAL
00227   double maxDepth    = depthToHCAL + depthHCAL - 1.1 * depthStep;
00228   double depthStart  = std::log(1./random->flatShoot()); // starting point lambda unts
00229 
00230   if(e < emin) {
00231     if(debug)
00232       LogDebug("FastCalorimetry") << " FamosHDShower : e <emin ->  depthStart = 0" << std::endl; 
00233     depthStart = 0.;
00234   }
00235  
00236   if(depthStart > maxDepth) {
00237     if(debug) LogDebug("FastCalorimetry") << " FamosHDShower : depthStart too big ...   = " 
00238                    << depthStart << std::endl; 
00239     
00240     depthStart = maxDepth *  random->flatShoot();
00241     if( depthStart < 0.) depthStart = 0.;
00242     if(debug) LogDebug("FastCalorimetry") << " FamosHDShower : depthStart re-calculated = " 
00243                    << depthStart << std::endl; 
00244   }
00245   
00246   if( onECAL && e < emid ) {
00247     if((depthECAL - depthStart)/depthECAL > 0.2 && depthECAL > depthStep ) {
00248       
00249       depthStart = 0.5 * depthECAL * random->flatShoot();
00250       if(debug) 
00251         LogDebug("FastCalorimetry") << " FamosHDShower : small energy, "
00252              << " depthStart reduced to = " << depthStart << std::endl; 
00253       
00254     }
00255   }
00256 
00257   if( depthHCAL < depthStep) {
00258     if(debug) LogDebug("FastCalorimetry") << " FamosHDShower : depthHCAL  too small ... = " 
00259                    << depthHCAL << " depthStart -> forced to 0 !!!" 
00260                    << std::endl;
00261     depthStart = 0.;    
00262     nmoresteps = 0;
00263     
00264     if(depthECAL < depthStep) {
00265       nsteps = -1;
00266       LogInfo("FastCalorimetry") << " FamosHDShower : too small ECAL and HCAL depths - " 
00267            << " particle is lost !!! " << std::endl; 
00268     }
00269   }
00270 
00271 
00272 
00273   if(debug)
00274     LogDebug("FastCalorimetry") << " FamosHDShower  depths(lam) - "  << std::endl 
00275          << "          ECAL = " << depthECAL  << std::endl
00276          << "           GAP = " << depthGAP   << std::endl
00277          << "          HCAL = " << depthHCAL  << std::endl
00278          << "  starting point = " << depthStart << std::endl; 
00279 
00280   if( onEcal ) {
00281     if(debug) LogDebug("FastCalorimetry") << " FamosHDShower : onECAL" << std::endl;
00282     if(depthStart < depthECAL) {
00283       if(debug) LogDebug("FastCalorimetry") << " FamosHDShower : depthStart < depthECAL" << std::endl;
00284       if((depthECAL - depthStart)/depthECAL > 0.1 && depthECAL > depthStep) {
00285         if(debug) LogDebug("FastCalorimetry") << " FamosHDShower : enough space to make ECAL step"
00286                        << std::endl;
00287         //  ECAL - one step
00288         nsteps++; 
00289         sum1   += depthECAL;             // at the end of step
00290         sum2   += depthECAL-depthStart; 
00291         sum3   += sum2 * lambdaEM / x0EM;
00292         lamtotal.push_back(sum1);
00293         lamdepth.push_back(sum2);
00294         lamcurr.push_back(lambdaEM);
00295         lamstep.push_back(depthECAL-depthStart);
00296         x0depth.push_back(sum3);
00297         x0curr.push_back(x0EM);
00298         detector.push_back(1);
00299         mip = 0;        
00300 
00301         if(debug) LogDebug("FastCalorimetry") << " FamosHDShower : " << " in ECAL sum1, sum2 "
00302                        << sum1 << " " << sum2 << std::endl;
00303         
00304         //                           // Gap - no additional step after ECAL
00305         //                           // just move further to HCAL over the gap
00306         sum1   += depthGAP;          
00307         sum2   += depthGAP; 
00308         sum3   += depthGAPx0;
00309       }
00310       // Just shift starting point to HCAL
00311       else { 
00312         //      cout << " FamosHDShower : not enough space to make ECAL step" << std::endl;
00313         if(debug)  LogDebug("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;
00314 
00315         depthStart = depthToHCAL;
00316         sum1 += depthStart;     
00317       }
00318     }
00319     else { // GAP or HCAL   
00320       
00321       if(depthStart >= depthECAL &&  depthStart < depthToHCAL ) {
00322       depthStart = depthToHCAL; // just a shift to HCAL for simplicity
00323       }
00324       sum1 += depthStart;
00325 
00326       if(debug) LogDebug("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;
00327     }
00328   }
00329   else {   // Forward 
00330     if(debug)  LogDebug("FastCalorimetry") << " FamosHDShower : forward" << std::endl;
00331     sum1 += depthStart;
00332     transFactor = 0.5;   // makes narower tresverse size of shower     
00333   }
00334  
00335   for (int i = 0; i < nmoresteps ; i++) {
00336     sum1 += depthStep;
00337     if (sum1 > (depthECAL + depthGAP + depthHCAL)) break; 
00338     sum2 += depthStep;
00339     sum3 += sum2 * lambdaHD / x0HD;
00340     lamtotal.push_back(sum1);
00341     lamdepth.push_back(sum2);
00342     lamcurr.push_back(lambdaHD);
00343     lamstep.push_back(depthStep);
00344     x0depth.push_back(sum3);
00345     x0curr.push_back(x0HD);
00346     detector.push_back(3);
00347     nsteps++;
00348   }
00349   
00350 
00351   // Make fractions of energy and transverse radii at each step 
00352 
00353   if(nsteps > 0) { makeSteps(nsteps); }
00354 
00355 }

virtual HDShower::~HDShower ( ) [inline, virtual]

Definition at line 43 of file HDShower.h.

00043 {;}

bool HDShower::compute ( )

Compute the shower longitudinal and lateral development.

Definition at line 464 of file HDShower.cc.

References EcalHitMaker::addHit(), HcalHitMaker::addHit(), count, debug, detector, ReconstructionGR_cff::ecal, lat::endl(), eStep, RandomEngine::flatShoot(), EcalHitMaker::getPads(), hcalDepthFactor, i, index, indexFinder(), infinity, j, lamstep, lamtotal, LogDebug, lossesOpt, maxTRfactor, nspots, nTRsteps, phi, radius(), random, HLT_VtxMuL3::result, rlamStep, HcalHitMaker::setDepth(), HcalHitMaker::setSpotEnergy(), EcalHitMaker::setSpotEnergy(), StDecayID::status, theGrid, theHcalHitMaker, and transProb().

Referenced by CalorimetryManager::HDShowerSimulation().

00464                        {
00465   
00466   //  TimeMe theT("FamosHDShower::compute");
00467 
00468   bool status = false;
00469   int numLongit = eStep.size();
00470   if(debug)
00471     LogDebug("FastCalorimetry") << " FamosHDShower::compute - " 
00472             << " N_long.steps required : " << numLongit << std::endl;
00473 
00474   if(numLongit > 0) {
00475 
00476     status = true;    
00477     // Prepare the trsanverse probability function
00478     std::vector<double> Fhist;
00479     std::vector<double> rhist; 
00480     for (int j = 0; j < nTRsteps + 1; j++) {
00481       rhist.push_back(maxTRfactor * j / nTRsteps );  
00482       Fhist.push_back(transProb(maxTRfactor,1.,rhist[j]));
00483       if(debug == 3) 
00484         LogDebug("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
00485     }
00486     
00487     // Longitudinal steps
00488     for (int i = 0; i < numLongit ; i++) {
00489       
00490       double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
00491       // vary the longitudinal profile if needed
00492       if(detector[i] != 1) currentDepthL0 *= hcalDepthFactor;                     
00493       if(debug)
00494         LogDebug("FastCalorimetry") << " FamosHDShower::compute - detector = "
00495                                     << detector[i]
00496                                     << "  currentDepthL0 = " 
00497                                     << currentDepthL0 << std::endl;
00498       
00499       double maxTRsize   = maxTRfactor * rlamStep[i];     // in lambda units
00500       double rbinsize    = maxTRsize / nTRsteps; 
00501       double espot       = eStep[i] / (double)nspots[i];  // re-adjust espot
00502 
00503       if(espot > 2. || espot < 0. ) 
00504         LogDebug("FastCalorimetry") << " FamosHDShower::compute - unphysical espot = " 
00505              << espot << std::endl;
00506 
00507       int ecal = 0;
00508       if(detector[i] != 1) { 
00509         bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
00510         
00511         if(debug)
00512           LogDebug("FastCalorimetry") << " FamosHDShower::compute - status of " 
00513                << " theHcalHitMaker->setDepth(currentDepthL0) is " 
00514                << setHDdepth << std::endl;
00515         
00516         if(!setHDdepth) 
00517           {
00518             currentDepthL0 -= lamstep[i];
00519             setHDdepth =  theHcalHitMaker->setDepth(currentDepthL0);
00520           }
00521         if(!setHDdepth) continue;
00522         
00523         theHcalHitMaker->setSpotEnergy(espot);        
00524       }
00525       else {
00526         ecal = 1;
00527         bool status = theGrid->getPads(currentDepthL0);   
00528         
00529         if(debug)
00530           LogDebug("FastCalorimetry") << " FamosHDShower::compute - status of Grid = " 
00531                << status << std::endl;
00532         
00533         if(!status) continue; 
00534 
00535         theGrid->setSpotEnergy(espot);
00536       }  
00537 
00538       
00539       // Transverse distribition
00540       int nok   = 0;                          // counter of OK  
00541       int count = 0;
00542       int inf   = infinity;
00543       if(lossesOpt) inf = nspots[i];          // losses are enabled, otherwise
00544       // only OK points are counted ...
00545       for (int j = 0; j < inf; j++) {
00546         if(nok == nspots[i]) break;
00547         count ++;
00548         
00549         double prob   = random->flatShoot();
00550         int index     = indexFinder(prob,Fhist);
00551         double radius = rlamStep[i] * rhist[index] +
00552           random->flatShoot() * rbinsize; // in-bin  
00553         double phi = 2.*M_PI*random->flatShoot();
00554         
00555         if(debug == 2)
00556           LogDebug("FastCalorimetry") << std::endl << " FamosHDShower::compute " << " r = " << radius 
00557                << "    phi = " << phi << std::endl;
00558         
00559         bool result;
00560         if(ecal) {
00561           result = theGrid->addHit(radius,phi,0);
00562           
00563           if(debug == 2)
00564             LogDebug("FastCalorimetry") << " FamosHDShower::compute - " 
00565                  << " theGrid->addHit result = " 
00566                  << result << std::endl;
00567         }
00568         else {
00569           result = theHcalHitMaker->addHit(radius,phi,0); 
00570           
00571           if(debug == 2)
00572             LogDebug("FastCalorimetry") << " FamosHDShower::compute - " 
00573                  << " theHcalHitMaker->addHit result = " 
00574                  << result << std::endl;
00575         }    
00576         if(result) nok ++; 
00577         
00578       } // end of tranverse simulation
00579       if(count == infinity) { 
00580         status = false; 
00581         if(debug)
00582           LogDebug("FastCalorimetry") << "*** FamosHDShower::compute " << " maximum number of" 
00583                << " transverse points " << count << " is used !!!" << std::endl; 
00584         break;
00585       }
00586 
00587       if(debug)
00588         LogDebug("FastCalorimetry") << " FamosHDShower::compute " << " long.step No." << i 
00589              << "   Ntry, Nok = " << count
00590              << " " << nok << std::endl; 
00591       
00592     } // end of longitudinal steps
00593   } // end of no steps
00594   return status;
00595 
00596 }

double HDShower::gam	(	double	x,
		double	a
	)			const `[inline, private]`

Definition at line 51 of file HDShower.h.

References funct::exp(), and funct::pow().

Referenced by makeSteps().

00051 { return pow(x,a-1.)*exp(-x); }

int HDShower::getmip ( ) [inline]

Definition at line 41 of file HDShower.h.

References mip.

Referenced by CalorimetryManager::HDShowerSimulation().

00041 {return mip;}

int HDShower::indexFinder	(	double	x,
		const std::vector< double > &	Fhist
	)			`[private]`

Definition at line 598 of file HDShower.cc.

References debug, lat::endl(), iter, LogDebug, size, and cmsRelvalreportInput::step.

Referenced by compute().

00598                                                                    {
00599   // binary search in the vector of doubles
00600   int size = Fhist.size();
00601 
00602   int curr = size / 2;
00603   int step = size / 4;
00604   int iter;
00605   int prevdir = 0; 
00606   int actudir = 0; 
00607 
00608   for (iter = 0; iter < size ; iter++) {
00609 
00610     if( curr >= size || curr < 1 )
00611       LogError("FastCalorimetry") << " FamosHDShower::indexFinder - wrong current index = " 
00612            << curr << " !!!" << std::endl;
00613 
00614     if ((x <= Fhist[curr]) && (x > Fhist[curr-1])) break;
00615     prevdir = actudir;
00616     if(x > Fhist[curr]) {actudir =  1;}
00617     else                {actudir = -1;}
00618     if(prevdir * actudir < 0) { if(step > 1) step /= 2;}
00619     curr += actudir * step;
00620     if(curr > size) curr = size;
00621     else { if(curr < 1) {curr = 1;}}
00622  
00623     if(debug == 3)
00624       LogDebug("FastCalorimetry") << " indexFinder - end of iter." << iter 
00625            << " curr, F[curr-1], F[curr] = "
00626            << curr << " " << Fhist[curr-1] << " " << Fhist[curr] << std::endl;
00627     
00628   }
00629 
00630   if(debug == 3)
00631     LogDebug("FastCalorimetry") << " indexFinder x = " << x << "  found index = " << curr-1
00632          << std::endl;
00633 
00634 
00635   return curr-1;
00636 }

void HDShower::makeSteps ( int nsteps ) [private]

Definition at line 357 of file HDShower.cc.

References aloge, alpEM, alpHD, and, balanceEH, betEM, betHD, count, debug, detector, e, lat::endl(), eSpotSize, eStep, RandomEngine::flatShoot(), gam(), i, infinity, int, lamcurr, lamdepth, lamstep, LogDebug, nspots, random, rlamStep, sum(), pyDBSRunClass::temp, tgamEM, tgamHD, theR1, theR2, theR3, transFactor, transParam, x, x0curr, x0depth, and y.

Referenced by HDShower().

00357                                    {
00358 
00359   double sumes = 0.;
00360   double sum   = 0.;
00361   std::vector<double> temp;
00362 
00363 
00364   if(debug)
00365     LogDebug("FastCalorimetry") << " FamosHDShower::makeSteps - " 
00366        << " nsteps required : " << nsteps << std::endl;
00367 
00368   int count = 0;
00369   for (int i = 0; i < nsteps; i++) {    
00370     
00371     double deplam = lamdepth[i] - 0.5 * lamstep[i];
00372     double depx0  = x0depth[i]  - 0.5 * lamstep[i] / x0curr[i]; 
00373     double     x = betEM * depx0;
00374     double     y = betHD * deplam;
00375     
00376     if(debug == 2)
00377       LogDebug("FastCalorimetry") << " FamosHDShower::makeSteps " 
00378                                   << " - step " << i
00379                                   << "   depx0, x = " << depx0 << ", " << x 
00380                                   << "   deplam, y = " << deplam << ", "
00381                                   << y << std::endl;
00382     
00383     double est = (part * betEM * gam(x,alpEM) * lamcurr[i] /
00384                   (x0curr[i] * tgamEM) + 
00385                   (1.-part) * betHD * gam(y,alpHD) / tgamHD) * lamstep[i];
00386     
00387     // protection ...
00388     if(est < 0.) {
00389       LogDebug("FastCalorimetry") << "*** FamosHDShower::makeSteps " << " - negative step energy !!!" 
00390            << std::endl;
00391       est = 0.;
00392       break ; 
00393     }
00394 
00395     // for estimates only
00396     sum += est;
00397     int nPest = (int) (est * e / sum / eSpotSize) ;
00398 
00399     if(debug == 2)
00400       LogDebug("FastCalorimetry") << " FamosHDShower::makeSteps - nPoints estimate = " 
00401            <<  nPest << std::endl;
00402 
00403     if(nPest <= 1 && count !=0 ) break;
00404 
00405     // good step - to proceed
00406 
00407     temp.push_back(est);
00408     sumes += est;
00409     
00410     rlamStep.push_back(transParam * (theR1 + (theR2 - theR3 * aloge))
00411                        * deplam * transFactor); 
00412     count ++;
00413   }
00414 
00415   // fluctuations in ECAL and re-distribution of remaining energy in HCAL
00416   if(detector[0] == 1 && count > 1) {
00417     double oldECALenergy = temp[0];
00418     double oldHCALenergy = sumes - oldECALenergy ;
00419     double newECALenergy = 2. * sumes;
00420     for (int i = 0; newECALenergy > sumes && i < infinity; i++)
00421       newECALenergy = 2.* balanceEH * random->flatShoot() * oldECALenergy; 
00422      
00423     if(debug == 2)
00424       LogDebug("FastCalorimetry") << "*** FamosHDShower::makeSteps " << " ECAL fraction : old/new - "
00425            << oldECALenergy/sumes << "/" << newECALenergy/sumes << std::endl;
00426 
00427     temp[0] = newECALenergy;
00428     double newHCALenergy = sumes - newECALenergy;
00429     double newHCALreweight =  newHCALenergy / oldHCALenergy;
00430 
00431     for (int i = 1; i < count; i++) {
00432       temp[i] *= newHCALreweight;
00433     }
00434     
00435   }
00436 
00437   
00438   // final re-normalization of the energy fractions  
00439   for (int i = 0; i < count ; i++) {
00440     eStep.push_back(temp[i] * e / sumes );
00441     nspots.push_back((int)(eStep[i]/eSpotSize)+1);
00442 
00443    if(debug)
00444      LogDebug("FastCalorimetry") << i << "  xO and lamdepth at the end of step = " 
00445           << x0depth[i] << " " 
00446           << lamdepth[i] << "   Estep func = " <<  eStep[i] 
00447           << "   Rstep = " << rlamStep[i] << "  Nspots = " <<  nspots[i]
00448           << std::endl; 
00449 
00450   }
00451 
00452   // The only step is in ECAL - let's make the size bigger ...  
00453   if(count == 1 and detector[0] == 1) rlamStep[0] *= 2.;
00454 
00455 
00456   if(debug) {
00457     if(eStep[0] > 0.95 * e && detector[0] == 1) 
00458       LogDebug("FastCalorimetry") << " FamosHDShower::makeSteps - " << "ECAL energy = " << eStep[0]
00459            << " out of total = " << e << std::endl;  
00460   }
00461 
00462 }

double HDShower::transProb	(	double	factor,
		double	R,
		double	r
	)			`[inline, private]`

Definition at line 56 of file HDShower.h.

Referenced by compute().

00056                                                       {
00057     double fsq = factor * factor; 
00058     return ((fsq + 1.)/fsq) * r * r / (r*r + R*R) ; 
00059   }

double HDShower::aloge [private]

Definition at line 79 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::alpEM [private]

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::alpHD [private]

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::balanceEH [private]

Definition at line 122 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::betEM [private]

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::betHD [private]

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::criticalEnergy [private]

Definition at line 118 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthStart [private]

Definition at line 78 of file HDShower.h.

Referenced by HDShower().

double HDShower::depthStep [private]

Definition at line 116 of file HDShower.h.

Referenced by HDShower().

std::vector<int> HDShower::detector [private]

Definition at line 81 of file HDShower.h.

Referenced by compute(), HDShower(), and makeSteps().

double HDShower::e [private]

Definition at line 101 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::eSpotSize [private]

Definition at line 114 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::eStep [private]

Definition at line 82 of file HDShower.h.

Referenced by compute(), and makeSteps().

double HDShower::hcalDepthFactor [private]

Definition at line 124 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::infinity [private]

Definition at line 86 of file HDShower.h.

Referenced by compute(), and makeSteps().

double HDShower::lambdaEM [private]

Definition at line 77 of file HDShower.h.

Referenced by HDShower().

double HDShower::lambdaHD [private]

Definition at line 77 of file HDShower.h.

Referenced by HDShower().

std::vector<double> HDShower::lamcurr [private]

Definition at line 84 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::lamdepth [private]

Definition at line 84 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::lamstep [private]

Definition at line 84 of file HDShower.h.

Referenced by compute(), HDShower(), and makeSteps().

std::vector<double> HDShower::lamtotal [private]

Definition at line 84 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::lossesOpt [private]

Definition at line 104 of file HDShower.h.

Referenced by compute(), and HDShower().

double HDShower::maxTRfactor [private]

Definition at line 120 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::mip [private]

Definition at line 98 of file HDShower.h.

Referenced by getmip(), and HDShower().

int HDShower::nDepthSteps [private]

Definition at line 106 of file HDShower.h.

Referenced by HDShower().

std::vector<int> HDShower::nspots [private]

Definition at line 81 of file HDShower.h.

Referenced by compute(), and makeSteps().

int HDShower::nTRsteps [private]

Definition at line 108 of file HDShower.h.

Referenced by compute(), and HDShower().

int HDShower::onEcal [private]

Definition at line 95 of file HDShower.h.

Referenced by HDShower().

double HDShower::part [private]

Definition at line 74 of file HDShower.h.

const RandomEngine* HDShower::random [private]

Definition at line 127 of file HDShower.h.

Referenced by compute(), HDShower(), and makeSteps().

std::vector<double> HDShower::rlamStep [private]

Definition at line 82 of file HDShower.h.

Referenced by compute(), and makeSteps().

double HDShower::tgamEM [private]

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::tgamHD [private]

Definition at line 74 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

const ECALProperties* HDShower::theECALproperties [private]

Definition at line 69 of file HDShower.h.

Referenced by HDShower().

EcalHitMaker* HDShower::theGrid [private]

Definition at line 89 of file HDShower.h.

Referenced by compute(), and HDShower().

HcalHitMaker* HDShower::theHcalHitMaker [private]

Definition at line 92 of file HDShower.h.

Referenced by compute().

const HCALProperties* HDShower::theHCALproperties [private]

Definition at line 70 of file HDShower.h.

Referenced by HDShower().

HDShowerParametrization* HDShower::theParam [private]

Definition at line 66 of file HDShower.h.

Referenced by HDShower().

double HDShower::theR1 [private]

Definition at line 73 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::theR2 [private]

Definition at line 73 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::theR3 [private]

Definition at line 73 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::transFactor [private]

Definition at line 112 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::transParam [private]

Definition at line 110 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::x0curr [private]

Definition at line 83 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

std::vector<double> HDShower::x0depth [private]

Definition at line 83 of file HDShower.h.

Referenced by HDShower(), and makeSteps().

double HDShower::x0EM [private]

Definition at line 77 of file HDShower.h.

Referenced by HDShower().

double HDShower::x0HD [private]

Definition at line 77 of file HDShower.h.

Referenced by HDShower().

HDShower Class Reference

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Types
typedef std::pair< XYZPoint, double >	Spot
typedef std::pair< unsigned int, double >	Step
typedef Steps::const_iterator	step_iterator
typedef std::vector< Step >	Steps
typedef math::XYZVector	XYZPoint
Public Member Functions
bool	compute ()
	Compute the shower longitudinal and lateral development.
int	getmip ()
	HDShower (const RandomEngine engine, HDShowerParametrization myParam, EcalHitMaker myGrid, HcalHitMaker myHcalHitMaker, int onECAL, double epart)
virtual	~HDShower ()
Private Member Functions
double	gam (double x, double a) const
int	indexFinder (double x, const std::vector< double > &Fhist)
void	makeSteps (int nsteps)
double	transProb (double factor, double R, double r)
Private Attributes
double	aloge
double	alpEM
double	alpHD
double	balanceEH
double	betEM
double	betHD
double	criticalEnergy
double	depthStart
double	depthStep
std::vector< int >	detector
double	e
double	eSpotSize
std::vector< double >	eStep
double	hcalDepthFactor
int	infinity
double	lambdaEM
double	lambdaHD
std::vector< double >	lamcurr
std::vector< double >	lamdepth
std::vector< double >	lamstep
std::vector< double >	lamtotal
int	lossesOpt
double	maxTRfactor
int	mip
int	nDepthSteps
std::vector< int >	nspots
int	nTRsteps
int	onEcal
double	part
const RandomEngine *	random
std::vector< double >	rlamStep
double	tgamEM
double	tgamHD
const ECALProperties *	theECALproperties
EcalHitMaker *	theGrid
HcalHitMaker *	theHcalHitMaker
const HCALProperties *	theHCALproperties
HDShowerParametrization *	theParam
double	theR1
double	theR2
double	theR3
double	transFactor
double	transParam
std::vector< double >	x0curr
std::vector< double >	x0depth
double	x0EM
double	x0HD