CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC4_patch1/src/FastSimulation/ShowerDevelopment/src/HDShower.cc

Go to the documentation of this file.
00001 //FastSimulation Headers
00002 #include "FastSimulation/ShowerDevelopment/interface/HDShower.h"
00003 //#include "FastSimulation/Utilities/interface/Histos.h"
00004 #include "FastSimulation/Utilities/interface/RandomEngine.h"
00005 
00006 
00008 // What's this?
00009 //#include "FastSimulation/FamosCalorimeters/interface/FASTCalorimeter.h"
00010 
00011 #include "FastSimulation/CaloHitMakers/interface/EcalHitMaker.h"
00012 #include "FastSimulation/CaloHitMakers/interface/HcalHitMaker.h"
00013 
00014 // CMSSW headers
00015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00016 
00018 // And This???? Doesn't seem to be needed
00019 // #include "Calorimetry/CaloDetector/interface/CellGeometry.h"
00020 
00021 #include "DQMServices/Core/interface/DQMStore.h"
00022 #include "DQMServices/Core/interface/MonitorElement.h"
00023 
00024 #include <cmath>
00025 
00026 // number attempts for transverse distribution if exit on a spec. condition
00027 #define infinity 10000
00028 // debugging flag ( 0, 1, 2, 3)
00029 #define debug 0
00030 
00031 using namespace edm;
00032 
00033 HDShower::HDShower(const RandomEngine* engine,
00034                    HDShowerParametrization* myParam, 
00035                    EcalHitMaker* myGrid,
00036                    HcalHitMaker* myHcalHitMaker,
00037                    int onECAL,
00038                    double epart,
00039            DQMStore * const dbeIn)
00040   : theParam(myParam), 
00041     theGrid(myGrid),
00042     theHcalHitMaker(myHcalHitMaker),
00043     onEcal(onECAL),
00044     e(epart),
00045     random(engine),
00046         dbe(dbeIn)
00047 { 
00048   // To get an access to constants read in FASTCalorimeter
00049   //  FASTCalorimeter * myCalorimeter= FASTCalorimeter::instance();
00050 
00051   // Values taken from FamosGeneric/FamosCalorimeter/src/FASTCalorimeter.cc
00052   lossesOpt       = myParam->hsParameters()->getHDlossesOpt();
00053   nDepthSteps     = myParam->hsParameters()->getHDnDepthSteps();
00054   nTRsteps        = myParam->hsParameters()->getHDnTRsteps();
00055   transParam      = myParam->hsParameters()->getHDtransParam();
00056   eSpotSize       = myParam->hsParameters()->getHDeSpotSize();
00057   depthStep       = myParam->hsParameters()->getHDdepthStep();
00058   criticalEnergy  = myParam->hsParameters()->getHDcriticalEnergy();
00059   maxTRfactor     = myParam->hsParameters()->getHDmaxTRfactor();
00060   balanceEH       = myParam->hsParameters()->getHDbalanceEH();
00061   hcalDepthFactor = myParam->hsParameters()->getHDhcalDepthFactor();
00062 
00063   // Special tr.size fluctuations 
00064   transParam *= (1. + random->flatShoot()); 
00065 
00066   // Special long. fluctuations
00067   hcalDepthFactor +=  0.05 * (2.* random->flatShoot() - 1.);
00068 
00069   transFactor = 1.;   // normally 1, in HF - might be smaller 
00070                       // to take into account
00071                       // a narrowness of the HF shower (Cherenkov light) 
00072 
00073   // simple protection ...
00074   if(e < 0) e = 0.;
00075 
00076   // Get the Famos Histos pointer
00077   //  myHistos = FamosHistos::instance();
00078   //  std::cout << " Hello FamosShower " << std::endl;
00079   
00080   theECALproperties = theParam->ecalProperties();
00081   theHCALproperties = theParam->hcalProperties();
00082 
00083   if (dbe) {
00084     dbe->cd();             
00085     if (!dbe->get("HDShower/ParticlesEnergy")) {}//std::cout << "NOT FOUND IN Shower.cc" << std::endl;}
00086     else {
00087       dbe->get("HDShower/ParticlesEnergy")->Fill(log10(e));
00088     }
00089   }
00090 
00091   double emax = theParam->emax(); 
00092   double emid = theParam->emid(); 
00093   double emin = theParam->emin(); 
00094   double effective = e;
00095 
00096   if( e < emid ) {
00097     theParam->setCase(1);
00098     // avoid "underflow" below Emin (for parameters calculation only)
00099     if(e < emin) effective = emin; 
00100   } 
00101   else 
00102     theParam->setCase(2);
00103  
00104   // Avoid "overflow" beyond Emax (for parameters)
00105   if(effective > 0.5 * emax) {
00106     eSpotSize *= 2.5;
00107     if(effective > emax) {
00108       effective = emax; 
00109       eSpotSize *= 4.; 
00110       depthStep *= 2.;
00111       if(effective > 1000.)
00112       eSpotSize *= 2.; 
00113     }
00114   }
00115 
00116   if(debug == 2 )
00117     LogInfo("FastCalorimetry") <<  " HDShower : " << std::endl 
00118          << "       Energy   "  <<               e << std::endl     
00119          << "      lossesOpt "  <<       lossesOpt << std::endl  
00120          << "    nDepthSteps "  <<     nDepthSteps << std::endl
00121          << "       nTRsteps "  <<        nTRsteps << std::endl
00122          << "     transParam "  <<      transParam << std::endl
00123          << "      eSpotSize "  <<       eSpotSize << std::endl
00124          << " criticalEnergy "  <<  criticalEnergy << std::endl
00125          << "    maxTRfactor "  <<     maxTRfactor << std::endl
00126          << "      balanceEH "  <<       balanceEH << std::endl
00127          << "hcalDepthFactor "  << hcalDepthFactor << std::endl;
00128 
00129 
00130   double alpEM1 = theParam->alpe1();
00131   double alpEM2 = theParam->alpe2();
00132 
00133   double betEM1 = theParam->bete1();
00134   double betEM2 = theParam->bete2();
00135 
00136   double alpHD1 = theParam->alph1();
00137   double alpHD2 = theParam->alph2();
00138 
00139   double betHD1 = theParam->beth1();
00140   double betHD2 = theParam->beth2();
00141 
00142   double part1 = theParam->part1();
00143   double part2 = theParam->part2();
00144 
00145   aloge = std::log(effective);
00146  
00147   double edpar = (theParam->e1() + aloge * theParam->e2()) * effective;
00148   double aedep = std::log(edpar);
00149 
00150   if(debug == 2)
00151     LogInfo("FastCalorimetry") << " HDShower : " << std::endl
00152          << "     edpar " <<   edpar << "   aedep " << aedep << std::endl 
00153          << "    alpEM1 " <<  alpEM1 << std::endl  
00154          << "    alpEM2 " <<  alpEM2 << std::endl  
00155          << "    betEM1 " <<  betEM1 << std::endl  
00156          << "    betEM2 " <<  betEM2 << std::endl  
00157          << "    alpHD1 " <<  alpHD1 << std::endl  
00158          << "    alpHD2 " <<  alpHD2 << std::endl  
00159          << "    betHD1 " <<  betHD1 << std::endl  
00160          << "    betHD2 " <<  betHD2 << std::endl  
00161          << "     part1 " <<   part1 << std::endl  
00162          << "     part2 " <<   part2 << std::endl; 
00163 
00164   // private members to set
00165   theR1  = theParam->r1();
00166   theR2  = theParam->r2();
00167   theR3  = theParam->r3();
00168 
00169   alpEM  = alpEM1 + alpEM2 * aedep;
00170   tgamEM = tgamma(alpEM);
00171   betEM  = betEM1 - betEM2 * aedep;
00172   alpHD  = alpHD1 + alpHD2 * aedep;
00173   tgamHD = tgamma(alpHD);
00174   betHD  = betHD1 - betHD2 * aedep;
00175   part   = part1  -  part2 * aedep;
00176   if(part > 1.) part = 1.;          // protection - just in case of 
00177 
00178   if(debug  == 2 )
00179     LogInfo("FastCalorimetry") << " HDShower : " << std::endl 
00180          << "    alpEM " <<  alpEM << std::endl  
00181          << "   tgamEM " << tgamEM << std::endl
00182          << "    betEM " <<  betEM << std::endl
00183          << "    alpHD " <<  alpHD << std::endl
00184          << "   tgamHD " << tgamHD << std::endl
00185          << "    betHD " <<  betHD << std::endl
00186          << "     part " <<   part << std::endl;
00187        
00188 
00189   if(onECAL){
00190     lambdaEM = theParam->ecalProperties()->interactionLength();
00191     x0EM     = theParam->ecalProperties()->radLenIncm();
00192   } 
00193   else {
00194     lambdaEM = 0.;
00195     x0EM     = 0.;
00196   }
00197   lambdaHD = theParam->hcalProperties()->interactionLength();
00198   x0HD     = theParam->hcalProperties()->radLenIncm();
00199 
00200   if(debug == 2)
00201     LogInfo("FastCalorimetry") << " HDShower e " << e        << std::endl
00202          << "          x0EM = " << x0EM     << std::endl 
00203          << "          x0HD = " << x0HD     << std::endl 
00204          << "         lamEM = " << lambdaEM << std::endl
00205          << "         lamHD = " << lambdaHD << std::endl;
00206 
00207 
00208   // Starting point of the shower
00209   // try first with ECAL lambda  
00210 
00211   double sum1   = 0.;  // lambda path from the ECAL/HF entrance;  
00212   double sum2   = 0.;  // lambda path from the interaction point;
00213   double sum3   = 0.;  // x0     path from the interaction point;
00214   int  nsteps   = 0;   // full number of longitudinal steps (counter);
00215 
00216   int nmoresteps;      // how many longitudinal steps in addition to 
00217                        // one (if interaction happens there) in ECAL
00218 
00219   mip = 1;             // just to initiate particle as MIP in ECAL    
00220 
00221   if(e < criticalEnergy ) nmoresteps = 1;  
00222   else                    nmoresteps = nDepthSteps;
00223   
00224   depthECAL  = 0.;
00225   depthGAP   = 0.;
00226   depthGAPx0 = 0.; 
00227   if(onECAL ) {
00228     depthECAL  = theGrid->ecalTotalL0();         // ECAL depth segment
00229         //depthECAL  = theGrid->ecalTotalL0() + theGrid->ps1TotalL0() + theGrid->ps2TotalL0() + theGrid->ps2eeTotalL0(); //TEST: include preshower
00230     depthGAP   = theGrid->ecalHcalGapTotalL0();  // GAP  depth segment
00231     depthGAPx0 = theGrid->ecalHcalGapTotalX0();  // GAP  depth x0
00232   }
00233   
00234   depthHCAL   = theGrid->hcalTotalL0();    // HCAL depth segment
00235   depthToHCAL = depthECAL + depthGAP;    
00236 
00237   //---------------------------------------------------------------------------
00238   // Depth simulation & various protections, among them
00239   // if too deep - get flat random in the allowed region
00240   // if no HCAL material behind - force to deposit in ECAL
00241   double maxDepth    = depthToHCAL + depthHCAL - 1.1 * depthStep;
00242   double depthStart  = std::log(1./random->flatShoot()); // starting point lambda unts
00243 
00244   if(e < emin) {
00245     if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : e <emin ->  depthStart = 0" << std::endl; 
00246     depthStart = 0.;
00247   }
00248  
00249   if(depthStart > maxDepth) {
00250     if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : depthStart too big ...   = " << depthStart << std::endl;     
00251     depthStart = maxDepth *  random->flatShoot();
00252     if(depthStart < 0.) depthStart = 0.;
00253     if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : depthStart re-calculated = " << depthStart << std::endl; 
00254   }
00255   
00256   if(onECAL && e < emid) {
00257     if(depthECAL > depthStep && (depthECAL - depthStart)/depthECAL > 0.2) {
00258       depthStart = 0.5 * depthECAL * random->flatShoot();
00259       if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : small energy, " << " depthStart reduced to = " << depthStart << std::endl; 
00260     }
00261   }
00262 
00263   if(depthHCAL < depthStep) {
00264     if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : depthHCAL  too small ... = " << depthHCAL << " depthStart -> forced to 0 !!!" << std::endl;
00265     depthStart = 0.;    
00266     nmoresteps = 0;
00267     
00268     if(depthECAL < depthStep) {
00269       nsteps = -1;
00270       LogInfo("FastCalorimetry") << " FamosHDShower : too small ECAL and HCAL depths - " << " particle is lost !!! " << std::endl; 
00271     }
00272   }
00273 
00274   if(debug)
00275     LogInfo("FastCalorimetry") << " FamosHDShower  depths(lam) - "  << std::endl 
00276          << "          ECAL = " << depthECAL  << std::endl
00277          << "           GAP = " << depthGAP   << std::endl
00278          << "          HCAL = " << depthHCAL  << std::endl
00279          << "  starting point = " << depthStart << std::endl; 
00280 
00281   if( onEcal ) {
00282     if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : onECAL" << std::endl;
00283     if(depthStart < depthECAL) {
00284       if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : depthStart < depthECAL" << std::endl;
00285       if(depthECAL > depthStep && (depthECAL - depthStart)/depthECAL > 0.1) {
00286             if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : enough space to make ECAL step" << std::endl;
00287             
00288             //  ECAL - one step
00289             nsteps++; 
00290             sum1   += depthECAL;             // at the end of step
00291             sum2   += depthECAL-depthStart; 
00292             sum3   += sum2 * lambdaEM / x0EM;
00293             lamtotal.push_back(sum1);
00294             lamdepth.push_back(sum2);
00295             lamcurr.push_back(lambdaEM);
00296             lamstep.push_back(depthECAL-depthStart);
00297             x0depth.push_back(sum3);
00298             x0curr.push_back(x0EM);
00299             detector.push_back(1);
00300             mip = 0;    
00301         
00302             if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : " << " in ECAL sum1, sum2 " << 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       else { // Just shift starting point to HCAL
00311             //  cout << " FamosHDShower : not enough space to make ECAL step" << std::endl;
00312             if(debug)  LogInfo("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;
00313 
00314             depthStart = depthToHCAL;
00315             sum1 += depthStart;     
00316       }
00317     }
00318     else { // GAP or HCAL         
00319       if(depthStart >= depthECAL &&  depthStart < depthToHCAL ) {
00320         depthStart = depthToHCAL; // just a shift to HCAL for simplicity
00321       }
00322       sum1 += depthStart;
00323       if(debug) LogInfo("FastCalorimetry") << " FamosHDShower : goto HCAL" << std::endl;
00324     }
00325   }
00326   else {   // Forward 
00327     if(debug)  LogInfo("FastCalorimetry") << " FamosHDShower : forward" << std::endl;
00328     sum1 += depthStart;
00329     transFactor = 0.5;   // makes narower tresverse size of shower     
00330   }
00331  
00332   for (int i = 0; i < nmoresteps ; i++) {
00333     sum1 += depthStep;
00334     if (sum1 > (depthECAL + depthGAP + depthHCAL)) break; 
00335     sum2 += depthStep;
00336     sum3 += sum2 * lambdaHD / x0HD;
00337     lamtotal.push_back(sum1);
00338     lamdepth.push_back(sum2);
00339     lamcurr.push_back(lambdaHD);
00340     lamstep.push_back(depthStep);
00341     x0depth.push_back(sum3);
00342     x0curr.push_back(x0HD);
00343     detector.push_back(3);
00344     nsteps++;
00345   }
00346 
00347   // Make fractions of energy and transverse radii at each step 
00348 
00349   if(nsteps > 0) { makeSteps(nsteps); }
00350 
00351 }
00352 
00353 void HDShower::makeSteps(int nsteps) {
00354 
00355   double sumes = 0.;
00356   double sum   = 0.;
00357   std::vector<double> temp;
00358 
00359   if(debug)
00360     LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - " 
00361        << " nsteps required : " << nsteps << std::endl;
00362 
00363   int count = 0;
00364   for (int i = 0; i < nsteps; i++) {    
00365     
00366     double deplam = lamdepth[i] - 0.5 * lamstep[i];
00367     double depx0  = x0depth[i]  - 0.5 * lamstep[i] / x0curr[i]; 
00368     double     x = betEM * depx0;
00369     double     y = betHD * deplam;
00370     
00371     if(debug == 2)
00372       LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps " 
00373                                   << " - step " << i
00374                                   << "   depx0, x = " << depx0 << ", " << x 
00375                                   << "   deplam, y = " << deplam << ", "
00376                                   << y << std::endl;
00377     
00378     double est = (part * betEM * gam(x,alpEM) * lamcurr[i] /
00379                   (x0curr[i] * tgamEM) + 
00380                   (1.-part) * betHD * gam(y,alpHD) / tgamHD) * lamstep[i];
00381     
00382     // protection ...
00383     if(est < 0.) {
00384       LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps " << " - negative step energy !!!" 
00385            << std::endl;
00386       est = 0.;
00387       break ; 
00388     }
00389 
00390     // for estimates only
00391     sum += est;
00392     int nPest = (int) (est * e / sum / eSpotSize) ;
00393 
00394     if(debug == 2)
00395       LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - nPoints estimate = " 
00396            <<  nPest << std::endl;
00397 
00398     if(nPest <= 1 && count !=0 ) break;
00399 
00400     // good step - to proceed
00401 
00402     temp.push_back(est);
00403     sumes += est;
00404     
00405     rlamStep.push_back(transParam * (theR1 + (theR2 - theR3 * aloge))
00406                        * deplam * transFactor); 
00407     count ++;
00408   }
00409 
00410   // fluctuations in ECAL and re-distribution of remaining energy in HCAL
00411   if(detector[0] == 1 && count > 1) {
00412     double oldECALenergy = temp[0];
00413     double oldHCALenergy = sumes - oldECALenergy ;
00414     double newECALenergy = 2. * sumes;
00415     for (int i = 0; newECALenergy > sumes && i < infinity; i++)
00416       newECALenergy = 2.* balanceEH * random->flatShoot() * oldECALenergy; 
00417      
00418     if(debug == 2)
00419       LogInfo("FastCalorimetry") << "*** FamosHDShower::makeSteps " << " ECAL fraction : old/new - "
00420            << oldECALenergy/sumes << "/" << newECALenergy/sumes << std::endl;
00421 
00422     temp[0] = newECALenergy;
00423     double newHCALenergy = sumes - newECALenergy;
00424     double newHCALreweight =  newHCALenergy / oldHCALenergy;
00425 
00426     for (int i = 1; i < count; i++) {
00427       temp[i] *= newHCALreweight;
00428     }
00429     
00430   }
00431   
00432   // final re-normalization of the energy fractions  
00433   for (int i = 0; i < count ; i++) {
00434     eStep.push_back(temp[i] * e / sumes );
00435     nspots.push_back((int)(eStep[i]/eSpotSize)+1);
00436 
00437    if(debug)
00438      LogInfo("FastCalorimetry") 
00439        << " step " << i
00440        << "  det: " << detector[i]   
00441        << "  xO and lamdepth at the end of step = " 
00442        << x0depth[i] << " "  
00443        << lamdepth[i] << "   Estep func = " <<  eStep[i]
00444        << "   Rstep = " << rlamStep[i] << "  Nspots = " <<  nspots[i]
00445        << "  espot = " <<  eStep[i] / (double)nspots[i]
00446        << std::endl; 
00447 
00448   }
00449 
00450   // The only step is in ECAL - let's make the size bigger ...  
00451   if(count == 1 and detector[0] == 1) rlamStep[0] *= 2.;
00452 
00453   if(debug) {
00454     if(eStep[0] > 0.95 * e && detector[0] == 1) 
00455       LogInfo("FastCalorimetry") << " FamosHDShower::makeSteps - " << "ECAL energy = " << eStep[0]
00456            << " out of total = " << e << std::endl;  
00457   }
00458 
00459 }
00460 
00461 bool HDShower::compute() {
00462   
00463   //  TimeMe theT("FamosHDShower::compute");
00464 
00465   bool status = false;
00466   int numLongit = eStep.size();
00467   if(debug)
00468     LogInfo("FastCalorimetry") << " FamosHDShower::compute - " 
00469             << " N_long.steps required : " << numLongit << std::endl;
00470 
00471   if(numLongit > 0) {
00472 
00473     status = true;    
00474     // Prepare the trsanverse probability function
00475     std::vector<double> Fhist;
00476     std::vector<double> rhist; 
00477     for (int j = 0; j < nTRsteps + 1; j++) {
00478       rhist.push_back(maxTRfactor * j / nTRsteps );  
00479       Fhist.push_back(transProb(maxTRfactor,1.,rhist[j]));
00480       if(debug == 3) LogInfo("FastCalorimetry") << "indexFinder - i, Fhist[i] = " << j << " " << Fhist[j] << std::endl;
00481     }
00482     
00483     // Longitudinal steps
00484     for (int i = 0; i < numLongit ; i++) {
00485       
00486       double currentDepthL0 = lamtotal[i] - 0.5 * lamstep[i];
00487       // vary the longitudinal profile if needed
00488       if(detector[i] != 1) currentDepthL0 *= hcalDepthFactor;                     
00489       if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute - detector = " << detector[i] << "  currentDepthL0 = " << currentDepthL0 << std::endl;
00490       
00491       double maxTRsize   = maxTRfactor * rlamStep[i];     // in lambda units
00492       double rbinsize    = maxTRsize / nTRsteps; 
00493       double espot       = eStep[i] / (double)nspots[i];  // re-adjust espot
00494 
00495       if(espot > 2. || espot < 0. ) LogInfo("FastCalorimetry") << " FamosHDShower::compute - unphysical espot = " << espot << std::endl;
00496 
00497       int ecal = 0;
00498       if(detector[i] != 1) { 
00499             bool setHDdepth = theHcalHitMaker->setDepth(currentDepthL0);
00500         
00501             if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of " << " theHcalHitMaker->setDepth(currentDepthL0) is " 
00502                                         << setHDdepth << std::endl;
00503         
00504             if(!setHDdepth) {
00505               currentDepthL0 -= lamstep[i];
00506               setHDdepth =  theHcalHitMaker->setDepth(currentDepthL0);
00507             }
00508             
00509                 if(!setHDdepth) continue;
00510 
00511             theHcalHitMaker->setSpotEnergy(espot);
00512                 
00513             //fill hcal longitudinal distribution histogram
00514         if (dbe) {
00515               dbe->cd();             
00516               if (!dbe->get("HDShower/LongitudinalShapeHCAL")) {}//std::cout << "NOT FOUND IN Shower.cc" << std::endl;}
00517               else {
00518                 //bins of 0.1 L0
00519                 double dt = 0.1;
00520                 // eStep is a real energy - scale by particle energy e
00521                         // subtract distance to hcal from current depth
00522                 dbe->get("HDShower/LongitudinalShapeHCAL")->Fill(currentDepthL0 - depthToHCAL, 1/e*eStep[i]/dt);
00523               }
00524         }
00525                 
00526       }
00527       else {
00528             ecal = 1;
00529             bool status = theGrid->getPads(currentDepthL0);   
00530         
00531             if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute - status of Grid = " << status << std::endl;
00532         
00533         if(!status) continue; 
00534 
00535         int ntry = nspots[i] * 10;  
00536         if( ntry >= infinity ) {  // use max allowed in case of too many spots
00537               nspots[i] = 0.5 * infinity;  
00538               espot *= 0.1 * (double)ntry / double(nspots[i]);
00539             }     
00540             else {
00541               espot *= 0.1;  // fine-grain energy spots in ECAL
00542                          // to avoid false ECAL clustering 
00543           nspots[i] = ntry;
00544             }
00545 
00546         theGrid->setSpotEnergy(espot);
00547                 
00548             //fill ecal longitudinal distribution histogram
00549         if (dbe) {
00550               dbe->cd();             
00551               if (!dbe->get("HDShower/LongitudinalShapeECAL")) {}//std::cout << "NOT FOUND IN Shower.cc" << std::endl;}
00552               else {
00553                 //bins of 0.1 L0
00554                 double dt = 0.1;
00555                 // eStep is a real energy - scale by particle energy e
00556                 dbe->get("HDShower/LongitudinalShapeECAL")->Fill(currentDepthL0, 1/e*eStep[i]/dt);
00557               }
00558         }
00559       }  
00560       
00561       // Transverse distribition
00562       int nok   = 0;                          // counter of OK  
00563       int count = 0;
00564       int inf   = infinity;
00565       if(lossesOpt) inf = nspots[i];          // if losses are enabled, otherwise
00566       // only OK points are counted ...
00567       if(nspots[i] > inf ) std::cout << " FamosHDShower::compute - at long.step " << i << "  too many spots required : "  <<  nspots[i] << " !!! " << std::endl;
00568 
00569       for (int j = 0; j < inf; j++) {
00570             if(nok == nspots[i]) break;
00571             count ++;
00572             
00573             double prob   = random->flatShoot();
00574             int index     = indexFinder(prob,Fhist);
00575             double radius = rlamStep[i] * rhist[index] + random->flatShoot() * rbinsize; // in-bin  
00576             double phi = 2.*M_PI*random->flatShoot();
00577         
00578             if(debug == 2)  LogInfo("FastCalorimetry") << std::endl << " FamosHDShower::compute " << " r = " << radius 
00579                                                 << "    phi = " << phi << std::endl;
00580         
00581             bool result;
00582             if(ecal) {
00583               result = theGrid->addHit(radius,phi,0);
00584           
00585               if(debug == 2) LogInfo("FastCalorimetry") << " FamosHDShower::compute - " << " theGrid->addHit result = " << result << std::endl;
00586                   
00587                   //fill ecal transverse distribution histogram
00588           if (dbe) {
00589             dbe->cd();             
00590                 if (!dbe->get("HDShower/TransverseShapeECAL")) {}//std::cout << "NOT FOUND IN Shower.cc" << std::endl;}
00591                 else {
00592                       double drho = 0.1;
00593                       // espot is a real energy - scale by particle energy
00594                       dbe->get("HDShower/TransverseShapeECAL")->Fill(radius,1/e*espot/drho);
00595                 }
00596               }
00597             }
00598             else {
00599               result = theHcalHitMaker->addHit(radius,phi,0); 
00600  
00601               if(debug == 2) LogInfo("FastCalorimetry") << " FamosHDShower::compute - " << " theHcalHitMaker->addHit result = " << result << std::endl;
00602                   
00603                   //fill hcal transverse distribution histogram
00604           if (dbe) {
00605             dbe->cd();             
00606                 if (!dbe->get("HDShower/TransverseShapeHCAL")) {}//std::cout << "NOT FOUND IN Shower.cc" << std::endl;}
00607                 else {
00608                       double drho = 0.1;
00609                       // espot is a real energy - scale by particle energy
00610                       dbe->get("HDShower/TransverseShapeHCAL")->Fill(radius,1/e*espot/drho);
00611                 }
00612               }
00613             }    
00614             
00615                 if(result) nok ++; 
00616 
00617 
00618 
00619       } // end of tranverse simulation
00620           
00621           
00622       if(count == infinity) { 
00623             if(debug) LogInfo("FastCalorimetry") << " FamosHDShower::compute " << " maximum number of" 
00624                                         << " transverse points " << count << " is used !!!" << std::endl; 
00625       }
00626 
00627       if(debug) LogInfo("FastCalorimetry")  << " FamosHDShower::compute " << " long.step No." 
00628                                 << i << "   Ntry, Nok = " << count << " " << nok << std::endl; 
00629     } // end of longitudinal steps
00630   } // end of no steps
00631 
00632   return status;
00633 
00634 }
00635 
00636 int HDShower::indexFinder(double x, const std::vector<double> & Fhist) {
00637   // binary search in the vector of doubles
00638   int size = Fhist.size();
00639 
00640   int curr = size / 2;
00641   int step = size / 4;
00642   int iter;
00643   int prevdir = 0; 
00644   int actudir = 0; 
00645 
00646   for (iter = 0; iter < size ; iter++) {
00647 
00648     if( curr >= size || curr < 1 )
00649       LogWarning("FastCalorimetry") << " FamosHDShower::indexFinder - wrong current index = " 
00650            << curr << " !!!" << std::endl;
00651 
00652     if ((x <= Fhist[curr]) && (x > Fhist[curr-1])) break;
00653     prevdir = actudir;
00654     if(x > Fhist[curr]) {actudir =  1;}
00655     else                {actudir = -1;}
00656     if(prevdir * actudir < 0) { if(step > 1) step /= 2;}
00657     curr += actudir * step;
00658     if(curr > size) curr = size;
00659     else { if(curr < 1) {curr = 1;}}
00660  
00661     if(debug == 3)
00662       LogInfo("FastCalorimetry") << " indexFinder - end of iter." << iter 
00663            << " curr, F[curr-1], F[curr] = "
00664            << curr << " " << Fhist[curr-1] << " " << Fhist[curr] << std::endl;
00665     
00666   }
00667 
00668   if(debug == 3)
00669     LogInfo("FastCalorimetry") << " indexFinder x = " << x << "  found index = " << curr-1
00670          << std::endl;
00671 
00672 
00673   return curr-1;
00674 }