Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:35:04 

0001 /** 
0002 */
0003 
0004 #include "Calibration/Tools/interface/InvMatrixUtils.h"
0005 #include "Calibration/Tools/interface/InvMatrixCommonDefs.h"
0006 #include "TStyle.h"
0007 #include "TROOT.h"
0008 #include "CLHEP/Geometry/Point3D.h"
0009 //#include "ConfigParser.h"
0010 #include "Calibration/Tools/interface/matrixSaver.h"
0011 
0012 #include <iostream>
0013 #include <string>
0014 #include <fstream>
0015 #include <sstream>
0016 #include <vector>
0017 #include <cassert>
0018 
0019 /** set the style for the printout*/
0020 void setStyle() {
0021   gROOT->SetStyle("Plain");
0022   gStyle->SetTextSize(0.5);
0023   //gStyle->SetOptStat (1111111) ;
0024   gStyle->SetOptStat(0);
0025   //gStyle->SetOptFit (1111) ;
0026   gStyle->SetOptFit(0);
0027   gStyle->SetTitleBorderSize(0);
0028   gStyle->SetTitleX(0.08);
0029   gStyle->SetTitleY(0.97);
0030   gStyle->SetPalette(1, nullptr);
0031   gStyle->SetStatColor(0);
0032   gStyle->SetFrameFillStyle(0);
0033   gStyle->SetFrameFillColor(0);
0034   return;
0035 }
0036 
0037 // -----------------------------------------------------------------
0038 
0039 TCanvas *getGlobalCanvas(std::string name) {
0040   setStyle();
0041   TCanvas *globalCanvas = static_cast<TCanvas *>(gROOT->FindObject(name.c_str()));
0042   if (globalCanvas) {
0043     globalCanvas->Clear();
0044     globalCanvas->UseCurrentStyle();
0045     globalCanvas->SetWindowSize(700, 600);
0046 
0047   } else {
0048     globalCanvas = new TCanvas(name.c_str(), name.c_str(), 700, 600);
0049   }
0050   return globalCanvas;
0051 }
0052 
0053 // -----------------------------------------------------------------
0054 
0055 TFile *getGlobalTFile(std::string name) {
0056   //    std::cout << "writing " << name << std::endl ;
0057   //    setStyle () ;
0058   TFile *globalTFile = (TFile *)gROOT->FindObject(name.c_str());
0059   if (!globalTFile) {
0060     //        std::cout << "does not exist. creating it " << std::endl;
0061     globalTFile = new TFile(name.c_str(), "RECREATE");
0062   }
0063 
0064   return globalTFile;
0065 }
0066 
0067 // -----------------------------------------------------------------
0068 
0069 int saveGlobalTFile(std::string name) {
0070   TFile *globalTFile = static_cast<TFile *>(gROOT->FindObject(name.c_str()));
0071   if (!globalTFile)
0072     return 1;
0073   globalTFile->Write();
0074   globalTFile->Close();
0075   delete globalTFile;
0076   return 0;
0077 }
0078 
0079 // -----------------------------------------------------------------
0080 
0081 CLHEP::HepMatrix *getSavedMatrix(const std::string &name) {
0082   matrixSaver reader;
0083   CLHEP::HepMatrix *savedMatrix;
0084   if (reader.touch(name)) {
0085     savedMatrix = static_cast<CLHEP::HepMatrix *>(reader.getMatrix(name));
0086   } else {
0087     savedMatrix = new CLHEP::HepMatrix(SCMaxEta, SCMaxPhi, 0);
0088   }
0089 
0090   return savedMatrix;
0091 }
0092 
0093 //========================================================================
0094 //da usare
0095 //HepGeom::Point3D<double> TBimpactPoint = TBposition (amplitude,m_beamEnergy) ;
0096 
0097 /* 
0098 dove trovare il codice di Chiara in CMSSW, per la ricostruzione
0099 della posizione:
0100 
0101 http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/RecoTBCalo/EcalTBAnalysisCoreTools/src/TBPositionCalc.cc?rev=1.1&cvsroot=CMSSW&content-type=text/vnd.viewcvs-markup
0102 
0103 */
0104 
0105 /** return the impact position of the electron over ECAL*/
0106 HepGeom::Point3D<Float_t> TBposition(const Float_t amplit[7][7],
0107                                      const Float_t beamEne,
0108                                      const Float_t w0,
0109                                      const Float_t x0,
0110                                      const Float_t a0,
0111                                      const Float_t sideX,  // crystal geometry, in mm
0112                                      const Float_t sideY) {
0113   // variables
0114   Float_t caloX = 0.;
0115   Float_t caloY = 0.;
0116   Float_t sumWeight = 0.;
0117   Float_t depth = x0 * (log(beamEne) + a0);  // shower depthh, in mm
0118   Float_t sin3 = 0.052335956;                // sin (3 degrees) , sin3 = sin(3.*3.141592654/180.)
0119 
0120   Float_t invE3x3 = 1. / get3x3(amplit);
0121 
0122   // loop over 3x3 crystals
0123   for (int eta = 2; eta <= 4; eta++) {
0124     for (int phi = 2; phi <= 4; phi++) {
0125       Float_t weight = log(amplit[eta][phi] * invE3x3) + w0;
0126       if (weight > 0) {
0127         caloX += (eta - 3) * sideX * weight;
0128         caloY -= (phi - 3) * sideY * weight;
0129         sumWeight += weight;
0130       }
0131     }
0132   }
0133 
0134   caloX /= sumWeight;
0135   caloY /= sumWeight;
0136 
0137   // correction for shower depthh
0138   caloX -= depth * sin3;
0139   caloY -= depth * sin3;
0140 
0141   // FIXME the z set to zero
0142   HepGeom::Point3D<Float_t> TBposition(caloX, caloY, 0);
0143 
0144   return TBposition;
0145 }
0146 
0147 // -----------------------------------------------------------------
0148 
0149 /** get the energy in the 5x5 
0150 from the 7x7 array around the most energetic crystal*/
0151 double get5x5(const Float_t energy[7][7]) {
0152   double total = 0.;
0153 
0154   for (int eta = 1; eta < 6; ++eta)
0155     for (int phi = 1; phi < 6; ++phi)
0156       total += energy[eta][phi];
0157 
0158   return total;
0159 }
0160 
0161 // -----------------------------------------------------------------
0162 
0163 /** get the energy in the 3x3 
0164 from the 7x7 array around the most energetic crystal*/
0165 double get3x3(const Float_t energy[7][7]) {
0166   double total = 0.;
0167 
0168   for (int eta = 2; eta < 5; ++eta)
0169     for (int phi = 2; phi < 5; ++phi)
0170       total += energy[eta][phi];
0171 
0172   return total;
0173 }
0174 
0175 // -----------------------------------------------------------------
0176 
0177 /**to get the parameters from a congiguration file*/
0178 /* int parseConfigFile (const TString& config)
0179 {
0180     if (gConfigParser) return 1 ;
0181     
0182     std::cout << "parsing "
0183         << config << " file"
0184         << std::endl ;
0185     
0186     
0187     gConfigParser = new ConfigParser () ;
0188     if ( !(gConfigParser->init(config)) )
0189         {
0190         std::cout << "Analysis::parseConfigFile: Could not open configuration file "
0191         << config << std::endl;
0192         perror ("Analysis::parseConfigFile: ") ;
0193         exit (-1) ;
0194         }
0195     gConfigParser->print () ;
0196     return 0 ;
0197 }*/
0198 
0199 // -----------------------------------------------------------------
0200 
0201 // per un certo eta, il cristallo puo' essere qualsiasi intero
0202 // Calibrationtra xmin e xmax
0203 // lo posso fissare solo sfruttando anche phi
0204 int xtalFromEtaPhi(const int &myEta, const int &myPhi) {
0205   int xMin = 20 * myEta + 1;
0206   int xMax = 20 * (myEta + 1) + 1;
0207 
0208   int myCryst = 999999;
0209 
0210   for (int x = xMin; x < xMax; x++) {
0211     if (phiFromXtal(x) == myPhi)
0212       myCryst = x;
0213   }
0214   return myCryst;
0215 }
0216 
0217 // -----------------------------------------------------------------
0218 
0219 int xtalFromiEtaiPhi(const int &iEta, const int &iPhi) {
0220   assert(iEta >= 1);
0221   assert(iEta <= 85);
0222   assert(iPhi >= 1);
0223   assert(iPhi <= 20);
0224   return 20 * (iEta - 1) + 21 - iPhi;
0225 }
0226 
0227 // -----------------------------------------------------------------
0228 
0229 int etaFromXtal(const int &xtal) {
0230   //  return floor (static_cast<double> ((xtal-1) / 20)) ;
0231   return int(floor((xtal - 1) / 20));
0232 }
0233 
0234 // -----------------------------------------------------------------
0235 
0236 int phiFromXtal(const int &xtal) {
0237   int phi = (xtal - 1) - 20 * etaFromXtal(xtal);
0238   return (20 - phi - 1);
0239 }
0240 
0241 // -----------------------------------------------------------------
0242 
0243 int ietaFromXtal(const int &xtal) { return etaFromXtal(xtal) + 1; }
0244 
0245 // -----------------------------------------------------------------
0246 
0247 int iphiFromXtal(const int &xtal) { return phiFromXtal(xtal) + 1; }
0248 
0249 // -----------------------------------------------------------------
0250 
0251 int extract(std::vector<int> *output, const std::string &dati) {
0252   std::ifstream _dati(dati.c_str());
0253   // loop over the file
0254   while (!_dati.eof()) {
0255     // get the line
0256     std::string dataline;
0257     do {
0258       getline(_dati, dataline, '\n');
0259     } while (*dataline.begin() == '#');
0260     std::stringstream linea(dataline);
0261     // loop over the line
0262     while (!linea.eof()) {
0263       int buffer = -1;
0264       linea >> buffer;
0265       if (buffer != -1)
0266         output->push_back(buffer);
0267     }  // loop over the line
0268   }  // loop over the file
0269   return output->size();
0270 }
0271 
0272 // -----------------------------------------------------------------
0273 
0274 // FIXME questi eta, phi sono quelli della matrice CLHEP,
0275 // FIXME non quelli del super-modulo, giusto?
0276 int writeCalibTxt(const CLHEP::HepMatrix &AmplitudeMatrix,
0277                   const CLHEP::HepMatrix &SigmaMatrix,
0278                   const CLHEP::HepMatrix &StatisticMatrix,
0279                   std::string fileName) {
0280   // look for the reference crystal
0281   double reference = 0.;
0282   for (int eta = 0; eta < SCMaxEta; ++eta)
0283     for (int phi = 0; phi < SCMaxPhi; ++phi) {
0284       if (AmplitudeMatrix[eta][phi] && SigmaMatrix[eta][phi] < 100 /*FIXME sigmaCut*/) {
0285         reference = AmplitudeMatrix[eta][phi];
0286         std::cout << "[InvMatrixUtils][writeCalibTxt] reference crystal: "
0287                   << "(" << eta << "," << phi << ") -> " << reference << "\n";
0288         break;
0289       }
0290     }
0291   if (!reference) {
0292     std::cerr << "ERROR: no calibration coefficients found" << std::endl;
0293     return 1;
0294   }
0295 
0296   // open the file for output
0297   std::ofstream txt_outfile;
0298   txt_outfile.open(fileName.c_str());
0299   txt_outfile << "# xtal\tcoeff\tsigma\tevt\tisGood\n";
0300 
0301   // loop over the crystals
0302   for (int eta = 0; eta < SCMaxEta; ++eta)
0303     for (int phi = 0; phi < SCMaxPhi; ++phi) {
0304       int isGood = 1;
0305       if (AmplitudeMatrix[eta][phi] == 0)
0306         isGood = 0;
0307       if (SigmaMatrix[eta][phi] > 100 /*FIXME sigmaCut*/)
0308         isGood = 0;
0309       txt_outfile << xtalFromEtaPhi(eta, phi) << "\t" << AmplitudeMatrix[eta][phi] / reference << "\t"
0310                   << SigmaMatrix[eta][phi] << "\t" << StatisticMatrix[eta][phi] << "\t" << isGood << "\n";
0311     }
0312 
0313   // save and close the file
0314   txt_outfile.close();
0315   return 0;
0316 }
0317 
0318 // -----------------------------------------------------------------
0319 
0320 int writeCMSSWCoeff(const CLHEP::HepMatrix &amplMatrix,
0321                     double calibThres,
0322                     float ERef,
0323                     const CLHEP::HepMatrix &sigmaMatrix,
0324                     const CLHEP::HepMatrix &statisticMatrix,
0325                     std::string fileName,
0326                     std::string genTag,
0327                     std::string method,
0328                     std::string version,
0329                     std::string type) {
0330   // open the file for output
0331   std::ofstream txt_outfile;
0332   txt_outfile.open(fileName.c_str());
0333   txt_outfile << "1\n";   // super-module number
0334   txt_outfile << "-1\n";  // number of events
0335   txt_outfile << genTag << "\n";
0336   txt_outfile << method << "\n";
0337   txt_outfile << version << "\n";
0338   txt_outfile << type << "\n";
0339 
0340   double reference = ERef;
0341 
0342   // loop over crystals
0343   for (int eta = 0; eta < SCMaxEta; ++eta)
0344     for (int phi = 0; phi < SCMaxPhi; ++phi) {
0345       if (amplMatrix[eta][phi] <= calibThres)
0346         txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << 1 << "\t" << -1 << "\t" << -1 << "\t" << 0 << "\n";
0347       else
0348         txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << reference / amplMatrix[eta][phi] << "\t"
0349                     << sigmaMatrix[eta][phi] << "\t" << statisticMatrix[eta][phi] << "\t" << 1 << "\n";
0350     }  // loop over crystals
0351 
0352   // save and close the file
0353   txt_outfile.close();
0354   return 0;
0355 }
0356 
0357 // -----------------------------------------------------------------
0358 
0359 int writeCMSSWCoeff(const CLHEP::HepMatrix &amplMatrix,
0360                     double calibThres,
0361                     int etaRef,
0362                     int phiRef,
0363                     const CLHEP::HepMatrix &sigmaMatrix,
0364                     const CLHEP::HepMatrix &statisticMatrix,
0365                     std::string fileName,
0366                     std::string genTag,
0367                     std::string method,
0368                     std::string version,
0369                     std::string type) {
0370   // open the file for output
0371   std::ofstream txt_outfile;
0372   txt_outfile.open(fileName.c_str());
0373   txt_outfile << "1\n";   // super-module number
0374   txt_outfile << "-1\n";  // number of events
0375   txt_outfile << genTag << "\n";
0376   txt_outfile << method << "\n";
0377   txt_outfile << version << "\n";
0378   txt_outfile << type << "\n";
0379 
0380   if (amplMatrix[etaRef - 1][phiRef - 1] == 0) {
0381     std::cerr << "The reference crystal: (" << etaRef << "," << phiRef << ") is out of range\n";
0382     return 1;
0383   }
0384   double reference = amplMatrix[etaRef - 1][phiRef - 1];
0385 
0386   // loop over crystals
0387   for (int eta = 0; eta < SCMaxEta; ++eta)
0388     for (int phi = 0; phi < SCMaxPhi; ++phi) {
0389       if (amplMatrix[eta][phi] <= calibThres)
0390         txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << 1 << "\t" << -1 << "\t" << -1 << "\t" << 0 << "\n";
0391       else
0392         txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << reference / amplMatrix[eta][phi] << "\t"
0393                     << sigmaMatrix[eta][phi] << "\t" << statisticMatrix[eta][phi] << "\t" << 1 << "\n";
0394     }  // loop over crystals
0395 
0396   // save and close the file
0397   txt_outfile.close();
0398   return 0;
0399 }
0400 
0401 // -----------------------------------------------------------------
0402 
0403 int translateCoeff(const CLHEP::HepMatrix &calibcoeff,
0404                    const CLHEP::HepMatrix &sigmaMatrix,
0405                    const CLHEP::HepMatrix &statisticMatrix,
0406                    std::string SMnumber,
0407                    double calibThres,
0408                    std::string fileName,
0409                    std::string genTag,
0410                    std::string method,
0411                    std::string version,
0412                    std::string type) {
0413   // open the file for output
0414   std::ofstream txt_outfile;
0415   txt_outfile.open(fileName.c_str());
0416   txt_outfile << SMnumber << "\n";  // super-module number
0417   txt_outfile << "-1\n";            // number of events
0418   txt_outfile << genTag << "\n";
0419   txt_outfile << method << "\n";
0420   txt_outfile << version << "\n";
0421   txt_outfile << type << "\n";
0422 
0423   // loop over crystals
0424   for (int eta = 0; eta < SCMaxEta; ++eta)
0425     for (int phi = 0; phi < SCMaxPhi; ++phi) {
0426       if (calibcoeff[eta][phi] < calibThres) {
0427         txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << 1 << "\t" << -1 << "\t" << -1 << "\t" << 0 << "\n";
0428         std::cout << "[translateCoefff][" << SMnumber << "]\t WARNING crystal " << xtalFromiEtaiPhi(eta + 1, phi + 1)
0429                   << " calib coeff below threshold: "
0430                   << "\t" << 1 << "\t" << -1 << "\t" << -1 << "\t" << 0 << "\n";
0431       } else
0432         txt_outfile << xtalFromiEtaiPhi(eta + 1, phi + 1) << "\t" << calibcoeff[eta][phi] << "\t"
0433                     << sigmaMatrix[eta][phi] << "\t" << statisticMatrix[eta][phi] << "\t" << 1 << "\n";
0434     }  // loop over crystals
0435 
0436   // save and close the file
0437   txt_outfile.close();
0438   return 0;
0439 }
0440 
0441 // -----------------------------------------------------------------
0442 
0443 int readCMSSWcoeff(CLHEP::HepMatrix &calibcoeff, const std::string &inputFileName, double defaultVal) {
0444   std::ifstream CMSSWfile;
0445   CMSSWfile.open(inputFileName.c_str());
0446   std::string buffer;
0447   CMSSWfile >> buffer;
0448   CMSSWfile >> buffer;
0449   CMSSWfile >> buffer;
0450   CMSSWfile >> buffer;
0451   CMSSWfile >> buffer;
0452   CMSSWfile >> buffer;
0453   while (!CMSSWfile.eof()) {
0454     int xtalnum;
0455     CMSSWfile >> xtalnum;
0456     double coeff;
0457     CMSSWfile >> coeff;
0458     double buffer;
0459     CMSSWfile >> buffer;
0460     int good;
0461     CMSSWfile >> good;
0462     CMSSWfile >> good;
0463     if (!good)
0464       coeff = defaultVal;  //FIXME 0 o 1?
0465     calibcoeff[etaFromXtal(xtalnum)][phiFromXtal(xtalnum)] = coeff;
0466   }
0467   return 0;
0468 }
0469 
0470 // -----------------------------------------------------------------
0471 
0472 int readCMSSWcoeffForComparison(CLHEP::HepMatrix &calibcoeff, const std::string &inputFileName) {
0473   std::ifstream CMSSWfile;
0474   CMSSWfile.open(inputFileName.c_str());
0475   std::string buffer;
0476   CMSSWfile >> buffer;
0477   CMSSWfile >> buffer;
0478   CMSSWfile >> buffer;
0479   CMSSWfile >> buffer;
0480   CMSSWfile >> buffer;
0481   CMSSWfile >> buffer;
0482   while (!CMSSWfile.eof()) {
0483     int xtalnum;
0484     CMSSWfile >> xtalnum;
0485     double coeff;
0486     CMSSWfile >> coeff;
0487     double buffer;
0488     CMSSWfile >> buffer;
0489     int good;
0490     CMSSWfile >> good;
0491     CMSSWfile >> good;
0492     if (!good)
0493       coeff = 0.;  //FIXME 0 o 1?
0494     calibcoeff[etaFromXtal(xtalnum)][phiFromXtal(xtalnum)] = coeff;
0495   }
0496   return 0;
0497 }
0498 
0499 // -----------------------------------------------------------------
0500 
0501 TH1D *smartProfile(TH2F *strip, double width) {
0502   TProfile *stripProfile = strip->ProfileX();
0503 
0504   // (from FitSlices of TH2.h)
0505 
0506   double xmin = stripProfile->GetXaxis()->GetXmin();
0507   double xmax = stripProfile->GetXaxis()->GetXmax();
0508   int profileBins = stripProfile->GetNbinsX();
0509 
0510   std::string name = strip->GetName();
0511   name += "_smart";
0512   TH1D *prof = new TH1D(name.c_str(), strip->GetTitle(), profileBins, xmin, xmax);
0513 
0514   int cut = 0;  // minimum number of entries per fitted bin
0515   int nbins = strip->GetXaxis()->GetNbins();
0516   int binmin = 1;
0517   int ngroup = 1;  // bins per step
0518   int binmax = nbins;
0519 
0520   // loop over the strip bins
0521   for (int bin = binmin; bin <= binmax; bin += ngroup) {
0522     TH1D *hpy = strip->ProjectionY("_temp", bin, bin + ngroup - 1, "e");
0523     if (hpy == nullptr)
0524       continue;
0525     int nentries = Int_t(hpy->GetEntries());
0526     if (nentries == 0 || nentries < cut) {
0527       delete hpy;
0528       continue;
0529     }
0530 
0531     Int_t biny = bin + ngroup / 2;
0532 
0533     hpy->GetXaxis()->SetRangeUser(hpy->GetMean() - width * hpy->GetRMS(), hpy->GetMean() + width * hpy->GetRMS());
0534     prof->Fill(strip->GetXaxis()->GetBinCenter(biny), hpy->GetMean());
0535     prof->SetBinError(biny, hpy->GetRMS());
0536 
0537     delete hpy;
0538   }  // loop over the bins
0539 
0540   delete stripProfile;
0541   return prof;
0542 }
0543 
0544 // -----------------------------------------------------------------
0545 
0546 TH1D *smartGausProfile(TH2F *strip, double width) {
0547   TProfile *stripProfile = strip->ProfileX();
0548 
0549   // (from FitSlices of TH2.h)
0550 
0551   double xmin = stripProfile->GetXaxis()->GetXmin();
0552   double xmax = stripProfile->GetXaxis()->GetXmax();
0553   int profileBins = stripProfile->GetNbinsX();
0554 
0555   std::string name = strip->GetName();
0556   name += "_smartGaus";
0557   TH1D *prof = new TH1D(name.c_str(), strip->GetTitle(), profileBins, xmin, xmax);
0558 
0559   int cut = 0;  // minimum number of entries per fitted bin
0560   int nbins = strip->GetXaxis()->GetNbins();
0561   int binmin = 1;
0562   int ngroup = 1;  // bins per step
0563   int binmax = nbins;
0564 
0565   // loop over the strip bins
0566   for (int bin = binmin; bin <= binmax; bin += ngroup) {
0567     TH1D *hpy = strip->ProjectionY("_temp", bin, bin + ngroup - 1, "e");
0568     if (hpy == nullptr)
0569       continue;
0570     int nentries = Int_t(hpy->GetEntries());
0571     if (nentries == 0 || nentries < cut) {
0572       delete hpy;
0573       continue;
0574     }
0575 
0576     Int_t biny = bin + ngroup / 2;
0577 
0578     TF1 *gaussian =
0579         new TF1("gaussian", "gaus", hpy->GetMean() - width * hpy->GetRMS(), hpy->GetMean() + width * hpy->GetRMS());
0580     gaussian->SetParameter(1, hpy->GetMean());
0581     gaussian->SetParameter(2, hpy->GetRMS());
0582     hpy->Fit("gaussian", "RQL");
0583 
0584     hpy->GetXaxis()->SetRangeUser(hpy->GetMean() - width * hpy->GetRMS(), hpy->GetMean() + width * hpy->GetRMS());
0585     prof->Fill(strip->GetXaxis()->GetBinCenter(biny), gaussian->GetParameter(1));
0586     prof->SetBinError(biny, gaussian->GetParameter(2));
0587 
0588     delete gaussian;
0589     delete hpy;
0590   }  // loop over the bins
0591 
0592   delete stripProfile;
0593   return prof;
0594 }
0595 
0596 // -----------------------------------------------------------------
0597 
0598 TH1D *smartError(TH1D *strip) {
0599   double xmin = strip->GetXaxis()->GetXmin();
0600   double xmax = strip->GetXaxis()->GetXmax();
0601   int stripsBins = strip->GetNbinsX();
0602 
0603   std::string name = strip->GetName();
0604   name += "_error";
0605   TH1D *error = new TH1D(name.c_str(), strip->GetTitle(), stripsBins, xmin, xmax);
0606 
0607   int binmin = 1;
0608   int ngroup = 1;  // bins per step
0609   int binmax = stripsBins;
0610   for (int bin = binmin; bin <= binmax; bin += ngroup) {
0611     double dummyError = strip->GetBinError(bin);
0612     error->SetBinContent(bin, dummyError);
0613   }
0614   return error;
0615 }
0616 
0617 // -----------------------------------------------------------------
0618 
0619 double effectiveSigma(TH1F &histogram, int vSteps) {
0620   double totInt = histogram.Integral();
0621   int maxBin = histogram.GetMaximumBin();
0622   int maxBinVal = int(histogram.GetBinContent(maxBin));
0623   int totBins = histogram.GetNbinsX();
0624   double area = totInt;
0625   double threshold = 0;
0626   double vStep = maxBinVal / vSteps;
0627   int leftBin = 1;
0628   int rightBin = totBins - 1;
0629   //loop over the vertical range
0630   while (area / totInt > 0.683) {
0631     threshold += vStep;
0632     // loop toward the left
0633     for (int back = maxBin; back > 0; --back) {
0634       if (histogram.GetBinContent(back) < threshold) {
0635         leftBin = back;
0636         break;
0637       }
0638     }  // loop toward the left
0639 
0640     // loop toward the right
0641     for (int fwd = maxBin; fwd < totBins; ++fwd) {
0642       if (histogram.GetBinContent(fwd) < threshold) {
0643         rightBin = fwd;
0644         break;
0645       }
0646     }  // loop toward the right
0647     area = histogram.Integral(leftBin, rightBin);
0648   }  //loop over the vertical range
0649 
0650   histogram.GetXaxis()->SetRange(leftBin, rightBin);
0651   // double sigmaEff = histogram.GetRMS () ;
0652   double halfWidthRange = 0.5 * (histogram.GetBinCenter(rightBin) - histogram.GetBinCenter(leftBin));
0653   return halfWidthRange;
0654 }
0655 
0656 // -----------------------------------------------------------------
0657 
0658 std::pair<int, int> findSupport(TH1F &histogram, double thres) {
0659   int totBins = histogram.GetNbinsX();
0660   if (thres >= histogram.GetMaximum())
0661     return std::pair<int, int>(0, totBins);
0662 
0663   int leftBin = totBins - 1;
0664   // search from left for the minimum
0665   for (int bin = 1; bin < totBins; ++bin) {
0666     if (histogram.GetBinContent(bin) > thres) {
0667       leftBin = bin;
0668       break;
0669     }
0670   }  // search from left for the minimum
0671   int rightBin = 1;
0672   // search from right for the maximum
0673   for (int bin = totBins - 1; bin > 0; --bin) {
0674     if (histogram.GetBinContent(bin) > thres) {
0675       rightBin = bin;
0676       break;
0677     }
0678   }  // search from right for the maximum
0679   return std::pair<int, int>(leftBin, rightBin);
0680 }
0681 
0682 // -----------------------------------------------------------------
0683 
0684 void mtrTransfer(double output[SCMaxEta][SCMaxPhi], CLHEP::HepMatrix *input, double Default) {
0685   for (int eta = 0; eta < SCMaxEta; ++eta)
0686     for (int phi = 0; phi < SCMaxPhi; ++phi) {
0687       if ((*input)[eta][phi])
0688         output[eta][phi] = (*input)[eta][phi];
0689       else
0690         output[eta][phi] = Default;
0691     }
0692   return;
0693 }
0694 
0695 // -----------------------------------------------------------------
0696 
0697 double etaCorrE1E25(int eta) {
0698   double p0 = 0.807883;
0699   double p1 = 0.000182551;
0700   double p2 = -5.76961e-06;
0701   double p3 = 7.41903e-08;
0702   double p4 = -2.25384e-10;
0703 
0704   double corr;
0705   if (eta < 6)
0706     corr = p0;
0707   else
0708     corr = p0 + p1 * eta + p2 * eta * eta + p3 * eta * eta * eta + p4 * eta * eta * eta * eta;
0709   return corr / p0;
0710 }
0711 // -----------------------------------------------------------------
0712 
0713 double etaCorrE1E49(int eta) {
0714   double p0 = 0.799895;
0715   double p1 = 0.000235487;
0716   double p2 = -8.26496e-06;
0717   double p3 = 1.21564e-07;
0718   double p4 = -4.83286e-10;
0719 
0720   double corr;
0721   if (eta < 8)
0722     corr = p0;
0723   else
0724     corr = p0 + p1 * eta + p2 * eta * eta + p3 * eta * eta * eta + p4 * eta * eta * eta * eta;
0725   return corr / p0;
0726 }
0727 // -----------------------------------------------------------------
0728 
0729 double etaCorrE1E9(int eta) {
0730   if (eta < 4)
0731     return 1.0;
0732   // grazie Paolo
0733   double p0 = 0.834629;
0734   double p1 = 0.00015254;
0735   double p2 = -4.91784e-06;
0736   double p3 = 6.54652e-08;
0737   double p4 = -2.4894e-10;
0738 
0739   double corr;
0740   if (eta < 6)
0741     corr = p0;
0742   else
0743     corr = p0 + p1 * eta + p2 * eta * eta + p3 * eta * eta * eta + p4 * eta * eta * eta * eta;
0744   return corr / p0;
0745 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!