EcalSimAlgos/test/testEcalShape.cpp

0001 #include "SimCalorimetry/EcalSimAlgos/interface/EcalSimParameterMap.h"
0002 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0003 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0004 #include "SimCalorimetry/EcalSimAlgos/interface/APDShape.h"
0005 #include "SimCalorimetry/EcalSimAlgos/interface/EBShape.h"
0006 #include "SimCalorimetry/EcalSimAlgos/interface/EEShape.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include <iostream>
0009 #include <iomanip>
0010
0011 #include "TROOT.h"
0012 #include "TStyle.h"
0013 #include "TH1F.h"
0014 #include "TCanvas.h"
0015 #include "TF1.h"
0016
0017 int main() {
0018   edm::MessageDrop::instance()->debugEnabled = false;
0019
0020   const EcalSimParameterMap parameterMap;
0021
0022   // for the purpose of testing the EcalShape class, doesn't need to (should not) fetch a shape from the DB but just to assume the testbeam one
0023   const APDShape theAPDShape;
0024   const EBShape theEBShape;
0025   const EEShape theEEShape;
0026
0027   const int nsamp = 500;  // phase I hardcoded ECAL shapes arrays
0028   const int tconv = 10;   // kNBinsPerNSec = 10
0029
0030   const unsigned int histsiz = nsamp * tconv;
0031
0032   for (unsigned int i(0); i != 3; ++i) {
0033     const DetId id(0 == i || 2 == i ? (DetId)EBDetId(1, 1) : (DetId)EEDetId(1, 50, 1));
0034
0035     const EcalShapeBase* theShape(0 == i ? (EcalShapeBase*)&theEBShape
0036                                          : (1 == i ? (EcalShapeBase*)&theEEShape : (EcalShapeBase*)&theAPDShape));
0037
0038     const double ToM(theShape->timeOfMax());
0039     const double T0(theShape->timeOfThr());
0040     const double rT(theShape->timeToRise());
0041
0042     // standard display of the implemented shape function
0043     const int csize = 500;
0044     TCanvas* showShape = new TCanvas("showShape", "showShape", 2 * csize, csize);
0045     /*
0046
0047 //  const std::vector<double>& nt = theShape.getTimeTable();
0048 //  const std::vector<double>& ntd = theShape.getDerivTable();
0049
0050   TH1F* shape1 = new TH1F("shape1","Tabulated Ecal MGPA shape",histsiz,0.,(float)(histsiz));
0051   TH1F* deriv1 = new TH1F("deriv1","Tabulated Ecal MGPA derivative",histsiz,0.,(float)(histsiz));
0052
0053
0054
0055   std::cout << "interpolated ECAL pulse shape and its derivative \n" << std::endl;
0056   for ( unsigned int i = 0; i < histsiz; ++i )
0057   {
0058      const double time ( (i-0.5)*0.1 - T0 ) ;
0059      const double myShape ( theShape( time ) ) ;
0060      const double myDeriv ( theShape.derivative( time ) ) ;
0061      shape1->Fill((float)(i+0.5),(float)myShape );
0062      deriv1->Fill((float)(i+0.5),(float)myDeriv );
0063      std::cout << " time (ns) = " << std::fixed << std::setw(6) << std::setprecision(2) << time + T0 + 0.1
0064            << " shape = " << std::setw(11) << std::setprecision(8) << myShape
0065            << " derivative = " << std::setw(11) << std::setprecision(8) << myDeriv << std::endl;
0066   }
0067
0068   showShape->Divide(2,1);
0069   showShape->cd(1);
0070   shape1->Draw();
0071   showShape->cd(2);
0072   deriv1->Draw();
0073   showShape->SaveAs("EcalShape.jpg");
0074   showShape->Clear("");
0075
0076   delete shape1;
0077   delete deriv1;
0078 */
0079
0080     const std::string name(0 == i ? "Barrel" : (1 == i ? "Endcap" : "APD"));
0081
0082     std::cout << "\n ********************* " << name << "************************";
0083
0084     std::cout << "\n Maximum time from tabulated values = " << std::fixed << std::setw(6) << std::setprecision(2) << ToM
0085               << std::endl;
0086
0087     std::cout << "\n Tzero from tabulated values        = " << std::fixed << std::setw(6) << std::setprecision(2) << T0
0088               << std::endl;
0089
0090     std::cout << "\n Rising time from tabulated values  = " << std::fixed << std::setw(6) << std::setprecision(2) << rT
0091               << std::endl;
0092
0093     // signal used with the nominal parameters and no jitter
0094
0095     std::cout << "\n computed ECAL " << name << " pulse shape and its derivative (LHC timePhaseShift = 1) \n"
0096               << std::endl;
0097     const double tzero = rT - (parameterMap.simParameters(id).binOfMaximum() - 1.) * 25.;
0098     double x = tzero;
0099
0100     const std::string title1("Computed Ecal " + name + " MGPA shape");
0101     const std::string title2("Computed Ecal " + name + " MGPA derivative");
0102
0103     TH1F* shape2 = new TH1F("shape2", title1.c_str(), nsamp, 0., (float)nsamp);
0104     TH1F* deriv2 = new TH1F("deriv2", title2.c_str(), nsamp, 0., (float)nsamp);
0105     double y = 0.;
0106     double dy = 0.;
0107
0108     for (unsigned int i(0); i != histsiz; ++i) {
0109       y = (*theShape)(x);
0110       dy = theShape->derivative(x);
0111       shape2->Fill((float)(x - tzero), (float)y);
0112       deriv2->Fill((float)(x - tzero), (float)dy);
0113       std::cout << " time (ns) = " << std::fixed << std::setw(6) << std::setprecision(2) << x - tzero
0114                 << " shape = " << std::setw(11) << std::setprecision(5) << y << " derivative = " << std::setw(11)
0115                 << std::setprecision(5) << dy << std::endl;
0116       x = x + 1. / (double)tconv;
0117     }
0118
0119     for (unsigned int iSample(0); iSample != 10; ++iSample) {
0120       std::cout << (*theShape)(tzero + iSample * 25.0) << std::endl;
0121     }
0122
0123     showShape->Divide(2, 1);
0124     showShape->cd(1);
0125     gPad->SetGrid();
0126     shape2->GetXaxis()->SetNdivisions(10, kFALSE);
0127     shape2->Draw();
0128
0129     showShape->cd(2);
0130     gPad->SetGrid();
0131     deriv2->GetXaxis()->SetNdivisions(10, kFALSE);
0132     deriv2->Draw();
0133
0134     const std::string fname(name + "EcalShapeUsed.jpg");
0135     showShape->SaveAs(fname.c_str());
0136
0137     delete shape2;
0138     delete deriv2;
0139     delete showShape;
0140   }
0141
0142   return 0;
0143 }