Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:30:07 

0001 #include <TCanvas.h>
0002 #include <TChain.h>
0003 #include <TFile.h>
0004 #include <TFitResult.h>
0005 #include <TFitResultPtr.h>
0006 #include <TF1.h>
0007 #include <TH1D.h>
0008 #include <TH2D.h>
0009 #include <TLegend.h>
0010 #include <TPaveStats.h>
0011 #include <TPaveText.h>
0012 #include <TProfile.h>
0013 #include <TProfile2D.h>
0014 #include <TROOT.h>
0015 #include <TStyle.h>
0016 
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <iomanip>
0020 #include <iostream>
0021 #include <fstream>
0022 #include <string>
0023 #include <vector>
0024 
0025 // Compare plots with and without beamline counters
0026 void PlotHistTBCompare(std::string infile1, std::string infile2, 
0027                std::string text, std::string prefix, 
0028                int maxMod=16, bool save=false);
0029 // Plot energy distribution in a given layer
0030 void PlotHistTBHitEn(std::string infile, std::string text, std::string prefix, 
0031              int type, bool separate, int rebin, int maxMod=16, 
0032              bool save=false);
0033 // Plot general TB plots at SIM/Digi/Reco labels
0034 void PlotHistSimDigRec(std::string fname="TBGenSimDigiReco.root", 
0035                std::string dirnm="HGCalTBAnalyzer", int type=0,
0036                std::string prefix="EL32", int maxMod=16,
0037                bool save=false);
0038 // Class to manipulate the Tree produced by the analysis code
0039 // HGCTB l1(std::string infile, std::string outfile);
0040 // l1.Loop(std::string calFile, std::string wtFile, unsigned int nLayers,
0041 //         float beamE, float cutETotByEBeam, int ntotEvent, int clusterEn, 
0042 //         int nPrintEvent)
0043 
0044 // There are 3 utility classes:
0045 // HexTopology : Useful for clusterize within a wafer around a local maximum
0046 //               HexTopology topo(fine);
0047 //               unsigned int id0 = topo.localMax(ids, energies, layer);
0048 //               double enCluster = topo.cluster(ids, energies, id, extent);
0049 //               std::vector<int> idx = topo.neighbours(id, extent);
0050 //               std::pair<int,int> rowcol = topo.findRowCol(id);
0051 //               Here "fine" is a boolean set to true if the wafer has 240 cells
0052 //                    "ids"  is a list of cell indices and "energies" the
0053 //                     corresponding energy deposits
0054 //                    "id" is a specific cell index and "extent" is number of
0055 //                    rows/columns for a nxn matrix around the cell "id"
0056 // HexCellGeometry: Useful to get position in an unrotated wafer 
0057 //               HexCellGeometry hg(fine);
0058 //               std::pair<double,double> xy = hg.position(id);
0059 //               Here "fine" is a boolean set to true if the wafer has 240 cells
0060 //                    "id" is the cell index and "xy" position (x,y) in mm in
0061 //                    local coordinate system
0062 // HexTBGeometry: Useful to get position in 2017 TB setup
0063 //               HexTBGeometry hg;
0064 //               std::pair<double,double> xy = hg.position(iw,ic,ifEE);
0065 //               Here "iw" is the wafer index, "ic" is the cell index;
0066 //                    "ifEE" is a boolean saying - true/false for EE/FH
0067 
0068 class HexTopology {
0069 
0070 public :
0071   HexTopology(bool fine);
0072   virtual ~HexTopology() {}
0073 
0074   double             cluster(const std::vector<unsigned int>& ids,
0075                  const std::vector<float>& energy, 
0076                  unsigned int id, int extent);
0077   unsigned int       localMax(const std::vector<unsigned int>& ids,
0078                   const std::vector<float>& energy, 
0079                   unsigned int layer);
0080   std::vector<int>   neighbours(int cell, int extent);
0081   std::pair<int,int> findRowCol(int cell);
0082 private :
0083   std::vector<int> leftCell, cells, addCell, subCell;
0084   unsigned int     nrows;
0085 };
0086 
0087 HexTopology::HexTopology(bool fine) {
0088   int cellCoarse[15] = {2,5,8,11,12,11,12,11,12,11,12,11,8,5,2};
0089   int addCoarse[15] =  {1,0,1, 0, 1, 0, 1, 0, 1, 0, 1, 0,1,0,1};
0090   int cellFine[20] = {3,6,9,12,15,16,15,16,15,16,15,16,15,16,15,14,11,8,5,2};
0091   int addFine[20]  = {0,1,0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0,1,0,1};
0092   int fullCoarse(11), fullFine(15);
0093   
0094   int full(0);
0095   if (fine) {
0096     for (int k=0; k<20; ++k) {
0097       cells.push_back(cellFine[k]); addCell.push_back(addFine[k]);
0098     }
0099     full = fullFine;
0100   } else {
0101     for (int k=0; k<15; ++k) {
0102       cells.push_back(cellCoarse[k]); addCell.push_back(addCoarse[k]);
0103     }
0104     full = fullCoarse;
0105   }
0106   leftCell.push_back(0);
0107   nrows = cells.size();
0108   for (unsigned int k=0; k<nrows; ++k) {
0109     leftCell.push_back((leftCell[k]+cells[k]));
0110     int ii = (cells[k]-full -1)/2;
0111     subCell.push_back(ii);
0112   }
0113   /*
0114   std::cout << "Initialize Wafer of type " << fine << " with " << nrows
0115         << " rows and " << leftCell[nrows] << " cells" << std::endl;
0116   for (unsigned int k=0; k<nrows; ++k) 
0117     std::cout << "Row[" << k << "] with " << cells[k] << ":" << addCell[k]
0118           << " cells: leftmost " << leftCell[k] << " offset "
0119           << subCell[k] << std::endl;
0120   */
0121 }
0122 
0123 std::pair<int,int> HexTopology::findRowCol(int cell) {
0124   int row(-1), col(-1);
0125   if (cell >= leftCell[0] && cell < leftCell[nrows]) {
0126     for (unsigned int i=0; i<=nrows; ++i) {
0127       if (cell < leftCell[i]) {
0128     row = i-1; col = cell-leftCell[i-1]-subCell[i-1]-addCell[i-1]; break;
0129       }
0130     }
0131   }
0132 // std::cout << "Cell " << cell << " row " << row << " column " << col << std::endl;
0133   return std::pair<int,int>(row,col);
0134 }
0135 
0136 std::vector<int> HexTopology::neighbours(int cell, int extent) {
0137   int cols[11] = {0,-1,1,-2,2,-3,3,-4,4,-5,5};
0138   std::vector<int>   listCell;
0139   std::pair<int,int> rowCol = findRowCol(cell);
0140   if (rowCol.first >= 0) {
0141     int nadd = 2*extent + 1;
0142     int addc = addCell[rowCol.first];
0143     for (int i = 0; i <= extent; ++i) {
0144       for (int j = 0; j < nadd; ++j) {
0145     int row = rowCol.first + i;
0146     if (row >= 0 && row <= (int)(nrows)) {
0147       int adc = (addc > addCell[row]) ? addc : addCell[row];
0148       int col = rowCol.second+cols[j]+leftCell[row]+subCell[row]+adc;
0149       if (col >= leftCell[row] && col < leftCell[row+1]) 
0150         listCell.push_back(col);
0151     }
0152     if (i > 0) {
0153       int row = rowCol.first - i;
0154       if (row >= 0 && row <= (int)(nrows)) {
0155         int adc = (addc > addCell[row]) ? addc : addCell[row];
0156         int col = rowCol.second+cols[j]+leftCell[row]+subCell[row]+adc;
0157         if (col >= leftCell[row] && col < leftCell[row+1]) 
0158           listCell.push_back(col);
0159       }
0160     }
0161       }
0162       nadd--;
0163     }
0164   }
0165   return listCell;
0166 }
0167 
0168 double HexTopology::cluster(const std::vector<unsigned int>& ids,
0169                 const std::vector<float>& energies, 
0170                 unsigned int id, int extent) {
0171   double energy(0);
0172   int cell = id&0xFF;
0173   std::vector<int> listCell = neighbours(cell, extent);
0174   for (unsigned int k=0; k<ids.size(); ++k) {
0175     int idcell = ids[k]&0xFF;
0176     if ((std::find(listCell.begin(),listCell.end(),idcell) != listCell.end()) &&
0177     ((id&0xFFFFFF00) == (ids[k]&0xFFFFFF00))) energy += energies[k];
0178   }
0179   return energy;
0180 }
0181 
0182 unsigned int HexTopology::localMax(const std::vector<unsigned int>& ids,
0183                    const std::vector<float>& energy, 
0184                    unsigned int layer) {
0185   float        emax(0);
0186   unsigned int idmax(0);
0187   for (unsigned int k=0; k<ids.size(); ++k) {
0188     if (layer == ((ids[k]>>19)&0x7F)) {
0189       if (energy[k] > emax) {
0190     emax  = energy[k];
0191     idmax = ids[k];
0192       }
0193     }
0194   }
0195   return idmax;
0196 }
0197   
0198 class HexCellGeometry {
0199 
0200 public :
0201   HexCellGeometry(bool fine);
0202   virtual ~HexCellGeometry() {}
0203 
0204   std::pair<double,double> position(const int cell);
0205 
0206 private :
0207   std::vector<std::pair<double,double> > xypos;
0208 
0209 };
0210 
0211 HexCellGeometry::HexCellGeometry(bool fine) {
0212   const int nC(15), nF(20);
0213   int nCoarse(11), nyCoarse(-42), nFine(15), nyFine(-56);
0214   int cellCoarse[nC] = {2,5,8,11,12,11,12,11,12,11,12,11,8,5,2};
0215   int cellFine[nF] = {3,6,9,12,15,16,15,16,15,16,15,16,15,16,15,14,11,8,5,2};
0216   double wafer(123.7);
0217 
0218   int    rows = (fine) ? nF : nC;
0219   double cell = (fine) ? wafer/nFine : wafer/nCoarse;
0220   double dx   = 0.5*cell;
0221   double dy   = 0.5*dx*tan(30.0*M_PI/180.0);
0222   int    ny   = (fine) ? nyFine : nyCoarse;
0223   for (int ir = 0; ir < rows; ++ir) {
0224     int    column = (fine) ? cellFine[ir] : cellCoarse[ir];
0225     int    nx     = 1 - column;
0226     double ypos   = dy*ny;
0227     for (int ic = 0; ic<column; ++ic) {
0228       double xpos = dx*nx;
0229       nx += 2;
0230       xypos.push_back(std::pair<double,double>(xpos,ypos));
0231     }
0232     ny += 6;
0233   }
0234   std::cout << "Initialize HexCellGeometry for " << xypos.size() << " cells"
0235         << std::endl;
0236 }
0237 
0238 std::pair<double,double> HexCellGeometry::position(const int cell) {
0239   std::pair<double,double> xy(0,0);
0240   if (cell >= 0 && cell < (int)(xypos.size())) {
0241     xy = xypos[cell];
0242   } else {
0243     xy = std::pair<double,double>(0,0);
0244   }
0245   return xy;
0246 }
0247 
0248 
0249 class HexTBGeometry {
0250 
0251 public :
0252   HexTBGeometry();
0253   virtual ~HexTBGeometry() {}
0254   
0255   std::pair<double,double> position(const int wafer, const int cell, 
0256                     const bool ifEE) const;
0257   std::pair<double,double> position_cell(const int cell) const;
0258   std::pair<double,double> position_wafer(const int wafer) const;
0259 
0260 private :
0261   std::vector<std::pair<double,double> > xypos_cell;
0262   std::vector<std::pair<double,double> > xypos_wafer;
0263 
0264 };
0265 
0266 HexTBGeometry::HexTBGeometry() {
0267 
0268   // Cell Geometry
0269   const int nC(15), nxCoarse(11), nyCoarse(21);
0270   int cellCoarse[nC] = {2,5,8,11,12,11,12,11,12,11,12,11,8,5,2};
0271   double wafer(123.7);
0272 
0273   double cell = wafer/nxCoarse;
0274   double dx   = 0.5*cell;
0275   double dy   = 0.5*cell*tan(30.0*M_PI/180.0);
0276   int    ny(nyCoarse);
0277   for (int ir = 0; ir < nC; ++ir) {
0278     int    nx     = 1 - cellCoarse[ir];
0279     double ypos   = dy*ny;
0280     for (int ic = 0; ic<cellCoarse[ir]; ++ic) {
0281       double xpos = dx*nx;
0282       nx += 2;
0283       // The wafers are rotated by 90 degrees anticlockwise
0284       // so x becomes y and y becomes x.
0285       xypos_cell.push_back(std::pair<double,double>(ypos,xpos)); 
0286     }
0287     ny -= 3;
0288   }
0289   std::cout << "Initialize HexTBGeometry for cells with " << xypos_cell.size() 
0290         << " cells" << std::endl;
0291 
0292   // Wafer geometry
0293   const int nW(3), nyWafer(3);
0294   int cellWafer[nW] = {2,3,2};
0295   dx   = 0.5*wafer;
0296   dy   = 0.5*wafer*tan(30.0*M_PI/180.0);
0297   ny   = nyWafer;
0298   for (int ir = 0; ir < nW; ++ir) {
0299     int    nx     = 1 - cellWafer[ir];
0300     double ypos   = dy*ny;
0301     for (int ic = 0; ic<cellWafer[ir]; ++ic) {
0302       double xpos = dx*nx;
0303       nx += 2;
0304       // The wafers are rotated by 90 degrees anticlckwise
0305       // so x becomes y and y becomes x.
0306       xypos_wafer.push_back(std::pair<double,double>(ypos,xpos));
0307     }
0308     ny -= 3;
0309   }
0310   std::cout << "Initialize HexGeometry for wafers with " << xypos_wafer.size() 
0311         << " wafers" << std::endl;
0312 }
0313 
0314 std::pair<double,double> HexTBGeometry::position_cell(const int cell) const {
0315   std::pair<double,double> xy(0,0);
0316   if (cell >= 0 && cell < (int)(xypos_cell.size())) {
0317     xy = xypos_cell[cell];
0318   } else {
0319     xy = std::pair<double,double>(0,0);
0320   }
0321   return xy;
0322 }
0323 
0324 std::pair<double,double> HexTBGeometry::position_wafer(const int wafer) const {
0325   std::pair<double,double> xy(0,0);
0326   if (wafer >= 0 && wafer < (int)(xypos_wafer.size())) {
0327     xy = xypos_wafer[wafer];
0328   } else {
0329     xy = std::pair<double,double>(0,0);
0330   }
0331   return xy;
0332 }
0333 
0334 std::pair<double,double> HexTBGeometry::position(const int wafer,const int cell,
0335                          const bool ifEE) const {
0336   int wf = ifEE ? (wafer+3) : wafer;
0337   std::pair<double,double> xyw = position_wafer(wf);
0338   std::pair<double,double> xyc = position_cell(cell);
0339   double xx = xyw.first + xyc.first;
0340   double yy = xyw.second + xyc.second;
0341   return std::pair<double,double>(xx,yy);
0342 }
0343 
0344 class HGCTB {
0345 public :
0346   TTree                *fChain;   //!pointer to the analyzed TTree or TChain
0347   Int_t                 fCurrent; //!current Tree number in a TChain
0348 
0349   // Declaration of leaf types
0350   std::vector<float>        *simHitLayEn1EE;
0351   std::vector<float>        *simHitLayEn2EE;
0352   std::vector<float>        *simHitLayEn1FH;
0353   std::vector<float>        *simHitLayEn2FH;
0354   std::vector<float>        *simHitLayEn1BH;
0355   std::vector<float>        *simHitLayEn2BH;
0356   Double_t                   xBeam, yBeam, zBeam, pBeam;
0357   std::vector<unsigned int> *simHitCellIdEE;
0358   std::vector<float>        *simHitCellEnEE;
0359   std::vector<unsigned int> *simHitCellIdFH;
0360   std::vector<float>        *simHitCellEnFH;
0361   std::vector<unsigned int> *simHitCellIdBH;
0362   std::vector<float>        *simHitCellEnBH;
0363 
0364   // List of branches
0365   TBranch                   *b_simHitLayEn1EE;   //!
0366   TBranch                   *b_simHitLayEn2EE;   //!
0367   TBranch                   *b_simHitLayEn1FH;   //!
0368   TBranch                   *b_simHitLayEn2FH;   //!
0369   TBranch                   *b_simHitLayEn1BH;   //!
0370   TBranch                   *b_simHitLayEn2BH;   //!
0371   TBranch                   *b_xBeam;            //!
0372   TBranch                   *b_yBeam;            //!
0373   TBranch                   *b_zBeam;            //!
0374   TBranch                   *b_pBeam;            //!
0375   TBranch                   *b_simHitCellIdEE;   //!
0376   TBranch                   *b_simHitCellEnEE;   //!
0377   TBranch                   *b_simHitCellIdFH;   //!
0378   TBranch                   *b_simHitCellEnFH;   //!
0379   TBranch                   *b_simHitCellIdBH;   //!
0380   TBranch                   *b_simHitCellEnBH;   //!
0381 
0382   HGCTB(std::string inName, std::string outName);
0383   virtual ~HGCTB();
0384   virtual Int_t    Cut(Long64_t entry);
0385   virtual Int_t    GetEntry(Long64_t entry);
0386   virtual Long64_t LoadTree(Long64_t entry);
0387   virtual void     Init(TTree *tree);
0388   virtual void     Loop(std::string calFile, std::string wtFile, 
0389             unsigned int nLayers, float beamE, float cutETotByEBeam,
0390             int ntotEvent, int clusterEn=1, int nPrintEvent=100);
0391   virtual Bool_t   Notify();
0392   virtual void     Show(Long64_t entry = -1);
0393   std::string      outName_;
0394 };
0395 
0396 HGCTB::HGCTB(std::string inName, std::string outName) : fChain(0),
0397                             outName_(outName) {
0398   TFile      *file = new TFile(inName.c_str());
0399   TDirectory *dir  = (TDirectory*)file->FindObjectAny("HGCalTBAnalyzer");
0400   TTree      *tree = (TTree*)dir->Get("HGCTB");
0401   Init(tree);
0402 }
0403 
0404 HGCTB::~HGCTB() {
0405   if (!fChain) return;
0406   delete fChain->GetCurrentFile();
0407 }
0408 
0409 Int_t HGCTB::GetEntry(Long64_t entry) {
0410   // Read contents of entry.
0411   if (!fChain) return 0;
0412   return fChain->GetEntry(entry);
0413 }
0414 
0415 Long64_t HGCTB::LoadTree(Long64_t entry) {
0416   // Set the environment to read one entry
0417   if (!fChain) return -5;
0418   Long64_t centry = fChain->LoadTree(entry);
0419   if (centry < 0) return centry;
0420   if (!fChain->InheritsFrom(TChain::Class()))  return centry;
0421   TChain *chain = (TChain*)fChain;
0422   if (chain->GetTreeNumber() != fCurrent) {
0423     fCurrent = chain->GetTreeNumber();
0424     Notify();
0425   }
0426   return centry;
0427 }
0428 
0429 void HGCTB::Init(TTree *tree) {
0430   // The Init() function is called when the selector needs to initialize
0431   // a new tree or chain. Typically here the branch addresses and branch
0432   // pointers of the tree will be set.
0433   // It is normally not necessary to make changes to the generated
0434   // code, but the routine can be extended by the user if needed.
0435   // Init() will be called many times when running on PROOF
0436   // (once per file to be processed).
0437   
0438   // Set object pointer
0439   simHitLayEn1EE = 0;
0440   simHitLayEn2EE = 0;
0441   simHitLayEn1FH = 0;
0442   simHitLayEn2FH = 0;
0443   simHitLayEn1BH = 0;
0444   simHitLayEn2BH = 0;
0445   simHitCellIdEE = 0;
0446   simHitCellEnEE = 0;
0447   simHitCellIdFH = 0;
0448   simHitCellEnFH = 0;
0449   simHitCellIdBH = 0;
0450   simHitCellEnBH = 0;
0451   // Set branch addresses and branch pointers
0452   if (!tree) return;
0453   fChain = tree;
0454   fCurrent = -1;
0455   fChain->SetMakeClass(1);
0456   
0457   fChain->SetBranchAddress("simHitLayEn1EE", &simHitLayEn1EE,&b_simHitLayEn1EE);
0458   fChain->SetBranchAddress("simHitLayEn2EE", &simHitLayEn2EE,&b_simHitLayEn2EE);
0459   fChain->SetBranchAddress("simHitLayEn1FH", &simHitLayEn1FH,&b_simHitLayEn1FH);
0460   fChain->SetBranchAddress("simHitLayEn2FH", &simHitLayEn2FH,&b_simHitLayEn2FH);
0461   fChain->SetBranchAddress("simHitLayEn1BH", &simHitLayEn1BH,&b_simHitLayEn1BH);
0462   fChain->SetBranchAddress("simHitLayEn2BH", &simHitLayEn2BH,&b_simHitLayEn2BH);
0463   fChain->SetBranchAddress("xBeam",          &xBeam,         &b_xBeam);
0464   fChain->SetBranchAddress("yBeam",          &yBeam,         &b_yBeam);
0465   fChain->SetBranchAddress("zBeam",          &zBeam,         &b_zBeam);
0466   fChain->SetBranchAddress("pBeam",          &pBeam,         &b_pBeam);
0467   fChain->SetBranchAddress("simHitCellIdEE", &simHitCellIdEE,&b_simHitCellIdEE);
0468   fChain->SetBranchAddress("simHitCellEnEE", &simHitCellEnEE,&b_simHitCellEnEE);
0469   fChain->SetBranchAddress("simHitCellIdFH", &simHitCellIdFH,&b_simHitCellIdFH);
0470   fChain->SetBranchAddress("simHitCellEnFH", &simHitCellEnFH,&b_simHitCellEnFH);
0471   fChain->SetBranchAddress("simHitCellIdBH", &simHitCellIdBH,&b_simHitCellIdBH);
0472   fChain->SetBranchAddress("simHitCellEnBH", &simHitCellEnBH,&b_simHitCellEnBH);
0473   Notify();
0474 }
0475 
0476 Bool_t HGCTB::Notify() {
0477   // The Notify() function is called when a new file is opened. This
0478   // can be either for a new TTree in a TChain or when when a new TTree
0479   // is started when using PROOF. It is normally not necessary to make changes
0480   // to the generated code, but the routine can be extended by the
0481   // user if needed. The return value is currently not used.
0482   
0483   return kTRUE;
0484 }
0485 
0486 void HGCTB::Show(Long64_t entry) {
0487   // Print contents of entry.
0488   // If entry is not specified, print current entry
0489   if (!fChain) return;
0490   fChain->Show(entry);
0491 }
0492 
0493 Int_t HGCTB::Cut(Long64_t ) {
0494   // This function may be called from Loop.
0495   // returns  1 if entry is accepted.
0496   // returns -1 otherwise.
0497   return 1;
0498 }
0499 
0500 void HGCTB::Loop(std::string calFile, std::string wtFile, unsigned int nLayers,
0501          float beamE, float cutETotByEBeam, int ntotEvent, 
0502          int clusterEn, int nPrintEvent) {
0503   //   In a ROOT session, you can do:
0504   //      Root > .L HGCTB.C
0505   //      Root > HGCTB t
0506   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0507   //      Root > t.Show();       // Show values of entry 12
0508   //      Root > t.Show(16);     // Read and show values of entry 16
0509   //      Root > t.Loop();       // Loop on all entries
0510   //
0511 
0512   //     This is the loop skeleton where:
0513   //    jentry is the global entry number in the chain
0514   //    ientry is the entry number in the current Tree
0515   //  Note that the argument to GetEntry must be:
0516   //    jentry for TChain::GetEntry
0517   //    ientry for TTree::GetEntry and TBranch::GetEntry
0518   //
0519   //       To read only selected branches, Insert statements like:
0520   // METHOD1:
0521   //    fChain->SetBranchStatus("*",0);  // disable all branches
0522   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0523   // METHOD2: replace line
0524   //    fChain->GetEntry(jentry);       //read all branches
0525   //by  b_branchname->GetEntry(ientry); //read only this branch
0526 
0527   if (fChain == 0) return;
0528 
0529   // Read in calibration and layerweights
0530   std::vector<double> calFactor(nLayers,1.0), layerWeight(nLayers,1.0);
0531   ifstream infil1, infil2;
0532   unsigned int line1(0), line2(0);
0533   infil1.open(calFile.c_str());
0534   if (infil1.is_open()) {
0535     while (!infil1.eof()) {
0536       float MPV;      
0537       infil1 >> MPV;
0538       calFactor.push_back(MPV/1.3);
0539       std::cout << line1 << " MPV: calFactor(MIP) = " << MPV << " : " 
0540         << calFactor[line1] << std::endl;
0541       line1++;
0542       if (line1 >= nLayers) break;
0543     }
0544     infil1.close();
0545   }
0546   infil2.open(wtFile.c_str());
0547   if (infil2.is_open()) {
0548     while (!infil2.eof()) {
0549       float weight;      
0550       infil2 >> weight;
0551       layerWeight.push_back(weight);
0552       std::cout << "layerWeight[" << line2 << "] = " << layerWeight[line2] 
0553         << std::endl;
0554       line2++;
0555       if (line2 >= nLayers) break;
0556     }
0557     infil2.close();
0558   }
0559 
0560   //Create histograms
0561   TFile *fout = new TFile(outName_.c_str(), "RECREATE");
0562   int   nbins(40000);
0563   float rangeGeV(0.05);
0564   float eLayerMax= (nLayers > 5) ? rangeGeV/calFactor[5] : rangeGeV/calFactor[nLayers-1];
0565   TH1F* h_SimHitTotEGeV      = new TH1F("h_SimHitTotEdepGeV","",nbins,0.,rangeGeV*nLayers/2);
0566   TH1F* h_SimHitTotEMIP      = new TH1F("h_SimHitTotEdepMIP","",nbins,0.,eLayerMax*nLayers/2);
0567   TH1F* h_SimHitTotEtoBE     = new TH1F("h_SimHitTotEdeptoBeamE","",nbins,0.,rangeGeV*nLayers/(4*beamE));
0568   TH1F* h_SimHitTotEtoBEWt   = new TH1F("h_SimHitTotEdeptoBeamEWt","",nbins,0.,rangeGeV*nLayers/(4*beamE));
0569   TH1F* h_SimHitTotEMIPCut   = new TH1F("h_SimHitCutTEtoBETotEdepMIP","",nbins,0.,eLayerMax*nLayers/2);
0570   TH1F* h_SimHitTotEMIPCutWt = new TH1F("h_SimHitCutWtTEtoBETotEdepMIP","",nbins,0.,eLayerMax*nLayers/2);
0571   std::vector<TH1F*> h_SimHitLayerE, h_SimHitLayerECut, h_SimHitLayerECutWt;
0572   for (unsigned int i=0; i<nLayers; i++) {
0573     char str1[20];
0574     sprintf(str1, "%d", i+1);
0575     TString hstr1 = str1;
0576     h_SimHitLayerE.push_back(new TH1F("h_SimHitLayer"+hstr1+"E","",nbins,0.,eLayerMax));
0577     h_SimHitLayerECut.push_back(new TH1F("h_SimHitCutTEtoBELayer"+hstr1+"E","",nbins,0.,eLayerMax));
0578     h_SimHitLayerECutWt.push_back(new TH1F("h_SimHitCutWtTEtoBELayer"+hstr1+"E","",nbins,0.,eLayerMax));
0579   }
0580 
0581   HexTopology ht1(false);
0582 
0583   Long64_t nentries = fChain->GetEntriesFast();
0584   if (ntotEvent >= 0 && ntotEvent < nentries) nentries = ntotEvent;
0585 
0586   Long64_t nbytes = 0, nb = 0;
0587   for (Long64_t jentry=0; jentry<nentries;jentry++) {
0588     Long64_t ientry = LoadTree(jentry);
0589     if (ientry < 0) break;
0590     nb = fChain->GetEntry(jentry);   nbytes += nb;
0591     // if (Cut(ientry) < 0) continue;
0592     // Fill histograms
0593     if (jentry % nPrintEvent == 0) std::cout << "  " << jentry  << "  Events Processed... " << std::endl;
0594 
0595     float TotEdepGeV(0.), TotEdepMIP(0.), TotEdepGeVWt(0.), TotEdepMIPWt(0.);
0596     std::vector<double> clusterEn2E;
0597     unsigned int nLayMax = (simHitLayEn2EE->size() > nLayers) ? nLayers : simHitLayEn2EE->size();
0598     for (unsigned int ilayer=0; ilayer<nLayMax; ilayer++) {
0599       unsigned int locMaxId = ht1.localMax((*simHitCellIdEE),(*simHitCellEnEE),ilayer+1);
0600       double clusterE2 = ht1.cluster((*simHitCellIdEE),(*simHitCellEnEE),locMaxId,2);
0601       clusterEn2E.push_back(clusterE2);
0602       float edep  = (clusterEn == 0) ? simHitLayEn2EE->at(ilayer) : clusterE2;
0603       float CaliE = edep/calFactor[ilayer];
0604       TotEdepGeV   += edep;
0605       TotEdepGeVWt += (edep*layerWeight[ilayer]);
0606       TotEdepMIP   += CaliE;
0607       TotEdepMIPWt += CaliE*layerWeight[ilayer];
0608       h_SimHitLayerE[ilayer]->Fill(CaliE);
0609       if (jentry % nPrintEvent == 0) 
0610     std::cout << ilayer+1 << " local Max Id : cluster Energy 2 " 
0611           << locMaxId << " : " << clusterE2 << " : " 
0612           << (*simHitLayEn2EE)[ilayer] << std::endl;
0613     }
0614     h_SimHitTotEGeV->Fill(TotEdepGeV);
0615     h_SimHitTotEMIP->Fill(TotEdepMIP);
0616     h_SimHitTotEtoBE->Fill(TotEdepGeV/beamE);
0617     h_SimHitTotEtoBEWt->Fill(TotEdepGeVWt/beamE);
0618     if ((TotEdepGeV/beamE) > cutETotByEBeam) {
0619       h_SimHitTotEMIPCut->Fill(TotEdepMIP);
0620       h_SimHitTotEMIPCutWt->Fill(TotEdepMIPWt);
0621     }
0622     for (unsigned int ilayer=0; ilayer<nLayMax; ++ilayer) {
0623       float CaliE = (clusterEn == 0) ? (simHitLayEn2EE->at(ilayer)/calFactor[ilayer]) : (clusterEn2E[ilayer]/calFactor[ilayer]);
0624       if ((TotEdepGeV/beamE) > cutETotByEBeam) {
0625     h_SimHitLayerECut[ilayer]->Fill(CaliE);
0626     h_SimHitLayerECutWt[ilayer]->Fill(CaliE*layerWeight[ilayer]);
0627       }
0628     }
0629   }
0630   fout->cd(); 
0631   h_SimHitTotEGeV->Write();     h_SimHitTotEMIP->Write();
0632   h_SimHitTotEtoBE->Write();    h_SimHitTotEtoBEWt->Write();
0633   h_SimHitTotEMIPCut->Write();  h_SimHitTotEMIPCutWt->Write();
0634   for (unsigned int i=0; i<nLayers; i++) {
0635     h_SimHitLayerE[i]->Write();
0636     h_SimHitLayerECut[i]->Write();
0637     h_SimHitLayerECutWt[i]->Write();
0638   }
0639   fout->Write(); fout->Close();
0640 }
0641 
0642 void PlotHistTBCompare(std::string infile1, std::string infile2, 
0643                std::string text, std::string prefix, int maxMod,
0644                bool save) {
0645 
0646   std::string name1[6] = {"BeamP","SimHitEn","SimHitTm","DigiADC","DigiLng",
0647               "RecHitEn"};
0648   std::string titl1[6] = {"","SIM","SIM","DIGI","DIGI","RECO"};
0649   std::string xttl1[6] = {"Beam Momentum (GeV/c)", "Energy (GeV)", "Time (ns)",
0650               "ADC", "Layer #", "Energy (GeV)"};
0651   std::string yttl1[6] = {"Event","Hit","Hit","Hit","Hit"};
0652   int         type1[6] = {0,1,1,0,0,1};
0653   int         rebin[6] = {1,1,5,1,1,1};
0654   double      xmin1[6] = {0,0,0,0,0,0};
0655   double      xmax1[6] = {150,0.01,100.,100,15,1.0};
0656   std::string name2[2] = {"DigiOcc","RecHitOcc"};
0657   std::string xttl2[2] = {"x (cm)", "x (cm)"};
0658   std::string yttl2[2] = {"y (cm)", "y (cm)"};
0659   std::string name3[5] = {"SimHitLng","SimHitLng1", "SimHitLng2","RecHitLng",
0660               "RecHitLng1"};
0661   std::string titl3[5] = {"SIM","SIM","SIM","RECO","RECO"};
0662   std::string xttl3[5] = {"z (cm)","Layer #","Layer #","z (cm)","Layer #"};
0663   std::string yttl3[5] = {"Mean Energy (GeV)","Mean Energy (Gev)",
0664               "Mean Energy (GeV)","Mean Energy (GeV)",
0665               "Mean Energy (GeV)"};
0666   double      xmin3[5] = {10,0,0,0,0};
0667   double      xmax30[5]= {25,15,15,15,15};
0668   double      xmax31[5]= {25,50,50,50,50};
0669   std::string name4[2] = {"SimHitLat","RecHitLat"};
0670   std::string titl4[2] = {"SIM","RECO"};
0671   std::string xttl4[2] = {"x (mm)", "x (cm)"};
0672   std::string yttl4[2] = {"y (mm)", "y (cm)"};
0673   std::string detect[3]= {"HGCalEESensitive","HGCalHEFSensitive","AHCal"};
0674   std::string label[2] = {"Without","With"};
0675 
0676   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0677   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
0678   gStyle->SetOptTitle(0);
0679   gStyle->SetOptStat(110);  gStyle->SetOptFit(0);
0680   TFile      *file1 = new TFile(infile1.c_str());
0681   TFile      *file2 = new TFile(infile2.c_str());
0682   bool isopen = file2->IsOpen();
0683   char name[100], namep[100];
0684   int  color[2] = {2,4};
0685   int  marker[2]= {20,21};
0686   for (int k=0; k<6; ++k) {
0687     if (k == 0) {
0688       sprintf (name, "%s", name1[k].c_str());
0689     } else {
0690       sprintf (name, "%s%s", name1[k].c_str(), detect[0].c_str());
0691     }
0692     TH1D* hist1 = (TH1D*)file1->FindObjectAny(name);
0693     if (hist1 != 0) {
0694       TLegend *legend = new TLegend(0.50, 0.70, 0.90, 0.79);
0695       legend->SetFillColor(kWhite);
0696       sprintf (namep, "c_%s%s", name1[k].c_str(), prefix.c_str());
0697       TCanvas *pad = new TCanvas(namep, namep, 700, 500);
0698       pad->SetRightMargin(0.10);
0699       pad->SetTopMargin(0.10);
0700       if (type1[k] != 0) pad->SetLogy();
0701       if (rebin[k] > 1) hist1->Rebin(rebin[k]);
0702       hist1->GetXaxis()->SetTitle(xttl1[k].c_str());
0703       hist1->GetYaxis()->SetTitle(yttl1[k].c_str());
0704       hist1->GetXaxis()->SetLabelOffset(0.005);
0705       hist1->GetXaxis()->SetTitleOffset(1.40);
0706       hist1->GetYaxis()->SetLabelOffset(0.005);
0707       hist1->GetYaxis()->SetTitleOffset(1.40);
0708       hist1->GetXaxis()->SetRangeUser(xmin1[k],xmax1[k]);
0709       hist1->SetLineColor(color[0]);
0710       hist1->SetMarkerStyle(marker[0]);
0711       hist1->SetMarkerColor(color[0]);
0712       sprintf (namep, "%s (without beam line)", text.c_str());
0713       legend->AddEntry(hist1,namep,"lp");
0714       hist1->Draw();
0715       pad->Update();
0716       TPaveStats* st1 = (TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats");
0717       if (st1 != NULL) {
0718     st1->SetLineColor(color[0]);
0719     st1->SetTextColor(color[0]);
0720     st1->SetY1NDC(0.85); st1->SetY2NDC(0.90);
0721     st1->SetX1NDC(0.65); st1->SetX2NDC(0.90);
0722       }
0723       TPaveText *txt1 = new TPaveText(0.80,0.64,0.90,0.69,"blNDC");
0724       txt1->SetFillColor(0);
0725       char txt[100];
0726       sprintf (txt, "%s", titl1[k].c_str());
0727       txt1->AddText(txt);
0728       txt1->Draw("same");
0729       pad->Modified();
0730       pad->Update();
0731       if (isopen) {
0732     hist1 = (TH1D*)file2->FindObjectAny(name);
0733     if (hist1 != 0) {
0734       if (rebin[k] > 1) hist1->Rebin(rebin[k]);
0735       hist1->GetXaxis()->SetTitle(xttl1[k].c_str());
0736       hist1->GetYaxis()->SetTitle(yttl1[k].c_str());
0737       hist1->GetXaxis()->SetLabelOffset(0.005);
0738       hist1->GetXaxis()->SetTitleOffset(1.40);
0739       hist1->GetYaxis()->SetLabelOffset(0.005);
0740       hist1->GetYaxis()->SetTitleOffset(1.40);
0741       hist1->GetXaxis()->SetRangeUser(xmin1[k],xmax1[k]);
0742       hist1->SetLineColor(color[1]);
0743       hist1->SetMarkerStyle(marker[1]);
0744       hist1->SetMarkerColor(color[1]);
0745       hist1->Draw("sames");
0746       sprintf (namep, "%s (with beam line)", text.c_str());
0747       legend->AddEntry(hist1,namep,"lp");
0748       pad->Update();
0749       st1 = (TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats");
0750       if (st1 != NULL) {
0751         st1->SetLineColor(color[1]);
0752         st1->SetTextColor(color[1]);
0753         st1->SetY1NDC(0.80); st1->SetY2NDC(0.85);
0754         st1->SetX1NDC(0.65); st1->SetX2NDC(0.90);
0755       }
0756       pad->Modified();
0757       pad->Update();
0758     }
0759     legend->Draw("same");
0760     pad->Update();
0761       }
0762       if (save) {
0763     sprintf (name, "%s.jpg", pad->GetName());
0764     pad->Print(name);
0765       } 
0766     }
0767   }
0768 
0769   for (int k=0; k<5; ++k) {
0770     sprintf (name, "%s%s", name3[k].c_str(), detect[0].c_str());
0771     TProfile* hist1 = (TProfile*)file1->FindObjectAny(name);
0772     if (hist1 != 0) {
0773       TLegend *legend = new TLegend(0.50, 0.65, 0.90, 0.73);
0774       legend->SetFillColor(kWhite);
0775       sprintf (namep, "c_%s%s", name3[k].c_str(), prefix.c_str());
0776       TCanvas *pad = new TCanvas(namep, namep, 700, 500);
0777       pad->SetRightMargin(0.10);
0778       pad->SetTopMargin(0.10);
0779       hist1->GetXaxis()->SetTitle(xttl3[k].c_str());
0780       hist1->GetYaxis()->SetTitle(yttl3[k].c_str());
0781       hist1->GetXaxis()->SetLabelOffset(0.005);
0782       hist1->GetXaxis()->SetTitleOffset(1.40);
0783       hist1->GetYaxis()->SetLabelOffset(0.005);
0784       hist1->GetYaxis()->SetTitleOffset(1.40);
0785       double xmax3 = (maxMod == 4) ? xmax30[k] : xmax31[k];
0786       hist1->GetXaxis()->SetRangeUser(xmin3[k],xmax3);
0787       hist1->SetLineColor(color[0]);
0788       hist1->SetMarkerStyle(marker[0]);
0789       hist1->SetMarkerColor(color[0]);
0790       sprintf (namep, "%s (without beam line)", text.c_str());
0791       legend->AddEntry(hist1,namep,"lp");
0792       hist1->Draw();
0793       pad->Update();
0794       TPaveStats* st1 = (TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats");
0795       if (st1 != NULL) {
0796     st1->SetLineColor(color[0]);
0797     st1->SetTextColor(color[0]);
0798     st1->SetY1NDC(0.82); st1->SetY2NDC(0.90);
0799     st1->SetX1NDC(0.65); st1->SetX2NDC(0.90);
0800       }
0801       TPaveText *txt1 = new TPaveText(0.80,0.59,0.90,0.64,"blNDC");
0802       txt1->SetFillColor(0);
0803       char txt[100];
0804       sprintf (txt, "%s", titl3[k].c_str());
0805       txt1->AddText(txt);
0806       txt1->Draw("same");
0807       pad->Modified();
0808       pad->Update();
0809       if (isopen) {
0810     hist1 = (TProfile*)file2->FindObjectAny(name);
0811     if (hist1 != 0) {
0812       hist1->GetXaxis()->SetTitle(xttl3[k].c_str());
0813       hist1->GetYaxis()->SetTitle(yttl3[k].c_str());
0814       hist1->GetXaxis()->SetLabelOffset(0.005);
0815       hist1->GetXaxis()->SetTitleOffset(1.40);
0816       hist1->GetYaxis()->SetLabelOffset(0.005);
0817       hist1->GetYaxis()->SetTitleOffset(1.40);
0818       double xmax3 = (maxMod == 4) ? xmax30[k] : xmax31[k];
0819       hist1->GetXaxis()->SetRangeUser(xmin3[k],xmax3);
0820       hist1->SetLineColor(color[1]);
0821       hist1->SetMarkerStyle(marker[1]);
0822       hist1->SetMarkerColor(color[1]);
0823       hist1->Draw("sames");
0824       sprintf (namep, "%s (with beam line)", text.c_str());
0825       legend->AddEntry(hist1,namep,"lp");
0826       pad->Update();
0827       st1 = (TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats");
0828       if (st1 != NULL) {
0829         st1->SetLineColor(color[1]);
0830         st1->SetTextColor(color[1]);
0831         st1->SetY1NDC(0.74); st1->SetY2NDC(0.82);
0832         st1->SetX1NDC(0.65); st1->SetX2NDC(0.90);
0833       }
0834       pad->Modified();
0835       pad->Update();
0836     }
0837     legend->Draw("same");
0838     pad->Update();
0839       }
0840       if (save) {
0841     sprintf (name, "%s.jpg", pad->GetName());
0842     pad->Print(name);
0843       } 
0844     }
0845   }
0846 
0847   for (int j=0; j<2; ++j) {
0848     for (int k=0; k<2; ++k) {
0849       sprintf (name, "%s%s", name2[k].c_str(), detect[0].c_str());
0850       TH2D* hist1(0);
0851       if (j==0) hist1 = (TH2D*)file1->FindObjectAny(name);
0852       else      hist1 = (TH2D*)file2->FindObjectAny(name);
0853       if (hist1 != 0) {
0854     TLegend *legend = new TLegend(0.10, 0.86, 0.50, 0.90);
0855     legend->SetFillColor(kWhite);
0856     sprintf (namep, "c_%s%s%s", name2[k].c_str(), prefix.c_str(),
0857          label[j].c_str());
0858     TCanvas *pad = new TCanvas(namep, namep, 500, 500);
0859     pad->SetRightMargin(0.10);
0860     pad->SetTopMargin(0.10);
0861     hist1->GetXaxis()->SetTitle(xttl2[k].c_str());
0862     hist1->GetYaxis()->SetTitle(yttl2[k].c_str());
0863     hist1->GetXaxis()->SetLabelOffset(0.005);
0864     hist1->GetXaxis()->SetTitleOffset(1.40);
0865     hist1->GetYaxis()->SetLabelOffset(0.005);
0866     hist1->GetYaxis()->SetTitleOffset(1.40);
0867     hist1->SetLineColor(color[j]);
0868     hist1->SetMarkerStyle(marker[j]);
0869     hist1->SetMarkerColor(color[j]);
0870     sprintf (namep, "%s (%s beam line)", text.c_str(), label[j].c_str());
0871     legend->AddEntry(hist1,namep,"lp");
0872     hist1->Draw("lego");
0873     pad->Update();
0874     TPaveStats* st1 = (TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats");
0875     if (st1 != NULL) {
0876       st1->SetLineColor(color[j]);
0877       st1->SetTextColor(color[j]);
0878       st1->SetY1NDC(0.82); st1->SetY2NDC(0.90);
0879       st1->SetX1NDC(0.65); st1->SetX2NDC(0.90);
0880     }
0881     TPaveText *txt1 = new TPaveText(0.80,0.76,0.90,0.81,"blNDC");
0882     txt1->SetFillColor(0);
0883     char txt[100];
0884     sprintf (txt, "%s", titl3[k].c_str());
0885     txt1->AddText(txt);
0886     txt1->Draw("same");
0887     pad->Modified();
0888     pad->Update();
0889     legend->Draw("same");
0890     pad->Update();
0891     if (save) {
0892       sprintf (name, "%s.jpg", pad->GetName());
0893       pad->Print(name);
0894     }
0895       }
0896     }
0897   }
0898 
0899   for (int j=0; j<2; ++j) {
0900     for (int k=0; k<2; ++k) {
0901       sprintf (name, "%s%s", name4[k].c_str(), detect[0].c_str());
0902       TProfile2D* hist1(0);
0903       if (j==0) hist1 = (TProfile2D*)file1->FindObjectAny(name);
0904       else      hist1 = (TProfile2D*)file2->FindObjectAny(name);
0905       if (hist1 != 0) {
0906     TLegend *legend = new TLegend(0.10, 0.86, 0.50, 0.90);
0907     legend->SetFillColor(kWhite);
0908     sprintf (namep, "c_%s%s%s", name4[k].c_str(), prefix.c_str(),
0909          label[j].c_str());
0910     TCanvas *pad = new TCanvas(namep, namep, 500, 500);
0911     pad->SetRightMargin(0.10);
0912     pad->SetTopMargin(0.10);
0913     hist1->GetXaxis()->SetTitle(xttl4[k].c_str());
0914     hist1->GetYaxis()->SetTitle(yttl4[k].c_str());
0915     hist1->GetXaxis()->SetLabelOffset(0.001);
0916     hist1->GetXaxis()->SetLabelSize(0.028);
0917     hist1->GetXaxis()->SetTitleOffset(1.40);
0918     hist1->GetYaxis()->SetLabelOffset(0.001);
0919     hist1->GetYaxis()->SetLabelSize(0.028);
0920     hist1->GetYaxis()->SetTitleOffset(1.40);
0921     hist1->SetLineColor(color[j]);
0922     hist1->SetMarkerStyle(marker[j]);
0923     hist1->SetMarkerColor(color[j]);
0924     sprintf (namep, "%s (%s beam line)", text.c_str(), label[j].c_str());
0925     legend->AddEntry(hist1,namep,"lp");
0926     hist1->Draw("lego");
0927     pad->Update();
0928     TPaveStats* st1 = (TPaveStats*)hist1->GetListOfFunctions()->FindObject("stats");
0929     if (st1 != NULL) {
0930       st1->SetLineColor(color[j]);
0931       st1->SetTextColor(color[j]);
0932       st1->SetY1NDC(0.82); st1->SetY2NDC(0.90);
0933       st1->SetX1NDC(0.65); st1->SetX2NDC(0.90);
0934     }
0935     TPaveText *txt1 = new TPaveText(0.80,0.76,0.90,0.81,"blNDC");
0936     txt1->SetFillColor(0);
0937     char txt[100];
0938     sprintf (txt, "%s", titl4[k].c_str());
0939     txt1->AddText(txt);
0940     txt1->Draw("same");
0941     pad->Modified();
0942     pad->Update();
0943     legend->Draw("same");
0944     pad->Update();
0945     if (save) {
0946       sprintf (name, "%s.jpg", pad->GetName());
0947       pad->Print(name);
0948     }
0949       }
0950     }
0951   }
0952 }
0953 
0954 void PlotHistTBHitEn(std::string infile, std::string text, std::string prefix, 
0955              int type, bool separate, int rebin, int maxMod, bool save){
0956 
0957   std::string name1[2] = {"SimHitEnA","SimHitEnB"};
0958   std::string title[2] = {"SIM Layer","Layer"};
0959   std::string xtitl[2] = {"Energy (GeV)", "Energy (GeV)"};
0960   std::string ytitl[2] = {"Tracks","Tracks"};
0961   std::string detect[3]= {"HGCalEESensitive","HGCalHESiliconSensitive","AHCal"};
0962   std::string label[2] = {"Without","With"};
0963   int         start[2] = {0, 1};
0964   int         nhmax0[2]= {12,4};
0965   int         nhmax1[2]= {48,16};
0966   int         color[6] = {1,2,4,7,6,9};
0967   int         ltype[6] = {1,1,1,1,1,1};
0968   int         mtype[6] = {20,21,22,23,24,33};
0969  
0970   if (type != 1) type = 0;
0971   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0972   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
0973   gStyle->SetOptTitle(0);
0974   gStyle->SetOptStat(1110);       gStyle->SetOptFit(0);
0975   TFile      *file = new TFile(infile.c_str());
0976 
0977   for (int k=0; k<2; ++k) {
0978     int nhmax = (maxMod == 4) ? nhmax0[k] : nhmax1[k];
0979     double dxs = 0.80/nhmax;
0980     if ((dxs > 0.15) || separate) dxs = 0.15;
0981     double dyl = (separate) ? 0.025 : 0.025*nhmax;
0982     double xhi(0.90), yhi(0.81);
0983     TCanvas* pad(0);
0984     TLegend* legend(0);
0985     char     name[100], namep[100];
0986     if (!separate) {
0987       sprintf (namep,"%s%s%s",name1[k].c_str(),prefix.c_str(),label[type].c_str());
0988       pad    = new TCanvas(namep, namep, 700, 500);
0989       pad->SetRightMargin(0.10); pad->SetTopMargin(0.10); pad->SetLogy();
0990       legend = new TLegend(0.70, yhi-dyl, 0.90, yhi);
0991       legend->SetFillColor(kWhite);
0992     }
0993     TPaveText *text1 = new TPaveText(0.60,yhi-dyl-0.04,0.90,yhi-dyl-0.01,"blNDC");
0994     text1->SetFillColor(0);
0995     sprintf (namep, "%s (%s beam line)",text.c_str(),label[type].c_str());
0996     text1->AddText(namep);
0997     bool first(true);
0998     for (int j=0; j<nhmax; ++j) {
0999       sprintf(name,"%s%d%s",name1[k].c_str(),(start[k]+j),detect[0].c_str());
1000       TH1D* hist = (TH1D*)file->FindObjectAny(name);
1001       if (hist != 0) {
1002     if (separate) {
1003       sprintf (namep,"%s%s%s%d",name1[k].c_str(),prefix.c_str(),label[type].c_str(),j);
1004       pad    = new TCanvas(namep, namep, 700, 500);
1005       pad->SetRightMargin(0.10); pad->SetTopMargin(0.10); pad->SetLogy();
1006       legend = new TLegend(0.70, 0.81-dyl, 0.90, 0.81);
1007       legend->SetFillColor(kWhite);
1008     }
1009     int j1 = j%6;
1010     int j2 = (j-j1)/6;
1011     hist->Rebin(rebin);
1012     hist->GetXaxis()->SetTitle(xtitl[k].c_str());
1013     hist->GetYaxis()->SetTitle(ytitl[k].c_str());
1014     hist->GetXaxis()->SetLabelOffset(0.005);
1015     hist->GetXaxis()->SetTitleOffset(1.40);
1016     hist->GetYaxis()->SetLabelOffset(0.005);
1017     hist->GetYaxis()->SetTitleOffset(1.40);
1018     hist->SetLineColor(color[j1]);
1019     hist->SetLineStyle(ltype[j2]);
1020     hist->SetMarkerStyle(mtype[j1]);
1021     hist->SetMarkerColor(color[j1]);
1022     double xmax = (maxMod == 4) ? 0.025 : 0.040;
1023     hist->GetXaxis()->SetRangeUser(0.0,xmax);
1024     sprintf (namep, "%s %d",title[k].c_str(),(j+1));
1025     legend->AddEntry(hist,namep,"lp");
1026     if (separate || first) hist->Draw();
1027     else                   hist->Draw("sames");
1028     first = false;
1029     pad->Update();
1030     TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1031     if (st1 != NULL) {
1032       st1->SetLineColor(color[j1]);
1033       st1->SetTextColor(color[j1]);
1034       st1->SetY1NDC(0.82); st1->SetY2NDC(0.90);
1035       st1->SetX1NDC(xhi-dxs); st1->SetX2NDC(xhi);
1036       if (!separate) xhi -= dxs;
1037     }
1038     pad->Update();
1039     if (separate) {
1040       legend->Draw("same");
1041       text1->Draw("same");
1042       pad->Update();
1043       if (save) {
1044         sprintf (name, "%s.jpg", pad->GetName());
1045         pad->Print(name);
1046       }
1047     }
1048       }
1049     }
1050     if (!separate) {
1051       legend->Draw("same");
1052       text1->Draw("same");
1053       pad->Update();
1054       if (save) {
1055     sprintf (name, "%s.jpg", pad->GetName());
1056     pad->Print(name);
1057       }
1058     }
1059   }
1060 }
1061 
1062 void PlotHistSimDigRec(std::string fname, std::string text, int type,
1063                std::string prefix, int maxMod, bool save) {
1064 
1065   std::string name1[6] = {"BeamP", "SimHitEnHGCalEESensitive", 
1066               "SimHitTmHGCalEESensitive",
1067               "DigiADCHGCalEESensitive",
1068               "DigiLngHGCalEESensitive",
1069               "RecHitEnHGCalEESensitive"};
1070   std::string name2[1] = {"DigiOccHGCalEESensitive"};
1071   std::string name3[5] = {"SimHitLngHGCalEESensitive",
1072               "SimHitLng1HGCalEESensitive",
1073               "SimHitLng2HGCalEESensitive",
1074               "RecHitLngHGCalEESensitive",
1075               "RecHitLng1HGCalEESensitive"};
1076   std::string name4[2] = {"SimHitLatHGCalEESensitive",
1077               "RecHitLatHGCalEESensitive"};
1078   double xrnglo1[6] = {  0.0,0.0,  0.0, 0.0, 0.0,0.0};
1079   double xrnghi1[6] = {150.0,0.1,200.0,20.0,50.0,0.1};
1080   double xrnglo2[1] = {-10.0};
1081   double xrnghi2[1] = { 10.0};
1082   double xrnglo3[5] = {10.0, 0.0, 0.0,10.0, 0.0};
1083   double xrnghi30[5]= {50.0,20.0,20.0,50.0,20.0};
1084   double xrnghi31[5]= {50.0,50.0,50.0,50.0,20.0};
1085   double xrnglo4[2] = {-40.0,-10.0};
1086   double xrnghi4[2] = { 40.0, 10.0};
1087   int    type1[6]   = {0,1, 1,1,0,1};
1088   int    type2[6]   = {1,1,10,1,1,1};
1089   std::string xtitl1[6] = {"Beam momentum", "SimHit energy (GeV)", 
1090                "Hit time (ns)", "ADC (Digi)", 
1091                "Longitudinal profile", "RecHit energy (GeV)"};
1092   std::string ytitl1[6] = {"Events", "Hits", "Hits", "Hits", " ", "Hits"};
1093   std::string title1[6] = {"", "SIM", "SIM", "DIGI", "DIGI", "RECO"};
1094   std::string xtitl3[5] = {"z (mm)", "Layer #", "Layer #", "z (cm)", "Layer #"};
1095   std::string ytitl3[5] = {"Mean Energy (GeV)", "Mean Energy (GeV)",
1096                "Mean Energy (GeV)", "Mean Energy (GeV)",
1097                "Mean Energy (GeV)"};
1098   std::string title3[5] = {"SIM", "SIM", "SIM", "RECO", "RECO"};
1099   std::string xtitl2[1] = {"x (cm)"};
1100   std::string ytitl2[1] = {"y (cm)"};
1101   std::string title2[1] = {"Digis"};
1102   std::string xtitl4[2] = {"x (mm)", "x (cm)"};
1103   std::string ytitl4[2] = {"y (mm)", "y (cm)"};
1104   std::string title4[2] = {"SimHits","RecHits"};
1105   std::string label[2] = {"Without","With"};
1106 
1107   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
1108   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
1109   gStyle->SetOptStat(1110);
1110   TFile      *file = new TFile(fname.c_str());
1111   if (file) {
1112     char name[100];
1113     TPaveText *text1 = new TPaveText(0.15,0.91,0.70,0.95,"blNDC");
1114     text1->SetFillColor(0);
1115     sprintf (name, "%s (%s beam line)",text.c_str(),label[type].c_str());
1116     text1->AddText(name);
1117     for (int k=0; k<6; ++k) {
1118       TH1D* hist = (TH1D*)file->FindObjectAny(name1[k].c_str());
1119       //    std::cout << name1[k] << " read out at " << hist << std::endl;
1120       if (hist != 0) {
1121     sprintf (name,"%s%s",prefix.c_str(),name1[k].c_str());
1122     TCanvas *pad = new TCanvas(name,name,500,500);
1123     pad->SetRightMargin(0.10); pad->SetTopMargin(0.10);
1124     hist->GetYaxis()->SetTitle(ytitl1[k].c_str());
1125     hist->GetXaxis()->SetTitle(xtitl1[k].c_str());
1126     hist->SetTitle(title1[k].c_str()); hist->Rebin(type2[k]);
1127     hist->GetXaxis()->SetRangeUser(xrnglo1[k],xrnghi1[k]);
1128     hist->GetYaxis()->SetTitleOffset(1.4);
1129     if (type1[k] == 1) pad->SetLogy();
1130     hist->Draw();
1131     pad->Update();
1132     TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1133     if (st1 != NULL) {
1134       st1->SetFillColor(0);
1135       st1->SetY1NDC(0.81); st1->SetY2NDC(0.90);
1136       st1->SetX1NDC(0.65); st1->SetX2NDC(0.90);
1137     }
1138     pad->Modified();
1139     pad->Update();
1140     if (save) {
1141       sprintf (name, "c_%s%s.gif", prefix.c_str(), name1[k].c_str());
1142       pad->Print(name);
1143     }
1144       }
1145     }
1146     for (int k=0; k<5; ++k) {
1147       TProfile* hist = (TProfile*)file->FindObjectAny(name3[k].c_str());
1148       //    std::cout << name3[k] << " read out at " << hist << std::endl;
1149       if (hist != 0) {
1150     sprintf (name,"%s%s",prefix.c_str(),name3[k].c_str());
1151     TCanvas *pad = new TCanvas(name,name,500,500);
1152     pad->SetRightMargin(0.10); pad->SetTopMargin(0.10);
1153     hist->GetYaxis()->SetTitle(ytitl3[k].c_str());
1154     hist->GetXaxis()->SetTitle(xtitl3[k].c_str());
1155     hist->SetMarkerStyle(20); hist->SetMarkerSize(0.8);
1156     hist->SetTitle(title3[k].c_str()); 
1157     double xrnghi3 = (maxMod == 4) ? xrnghi30[k] : xrnghi31[k];
1158     hist->GetXaxis()->SetRangeUser(xrnglo3[k],xrnghi3);
1159     hist->GetYaxis()->SetTitleOffset(1.4);
1160     hist->Draw();
1161     pad->Update();
1162     TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
1163     if (st1 != NULL) {
1164       st1->SetFillColor(0);
1165       st1->SetY1NDC(0.75); st1->SetY2NDC(0.90);
1166       st1->SetX1NDC(0.65); st1->SetX2NDC(0.90);
1167     }
1168     text1->Draw("same");
1169     pad->Modified();
1170     pad->Update();
1171     if (save) {
1172       sprintf (name, "c_%s%s.gif", prefix.c_str(), name3[k].c_str());
1173       pad->Print(name);
1174     }
1175       }
1176     }
1177     for (int k=0; k<1; ++k) {
1178       TH2D* hist = (TH2D*)file->FindObjectAny(name2[k].c_str());
1179       //    std::cout << name2[k] << " read out at " << hist << std::endl;
1180       if (hist != 0) {
1181     sprintf (name,"%s%s",prefix.c_str(),name2[k].c_str());
1182     TCanvas *pad = new TCanvas(name,name,500,500);
1183     pad->SetRightMargin(0.10); pad->SetTopMargin(0.10);
1184     hist->GetYaxis()->SetTitle(ytitl2[k].c_str());
1185     hist->GetXaxis()->SetTitle(xtitl2[k].c_str());
1186     hist->GetXaxis()->SetRangeUser(xrnglo2[k],xrnghi2[k]);
1187     hist->GetYaxis()->SetRangeUser(xrnglo2[k],xrnghi2[k]);
1188     hist->SetMarkerStyle(20); hist->SetMarkerSize(0.2);
1189     hist->SetTitle(title2[k].c_str());
1190     hist->GetXaxis()->SetTitleOffset(1.4);
1191     hist->GetYaxis()->SetTitleOffset(1.4);
1192     hist->Draw();
1193     text1->Draw("same");
1194     if (save) {
1195       sprintf (name, "c_%s%s.gif", prefix.c_str(), name2[k].c_str());
1196       pad->Print(name);
1197     }
1198       }
1199     }
1200     for (int k=0; k<2; ++k) {
1201       TProfile2D* hist = (TProfile2D*)file->FindObjectAny(name4[k].c_str());
1202       //    std::cout << name4[k] << " read out at " << hist << std::endl;
1203       if (hist != 0) {
1204     sprintf (name,"%s%s",prefix.c_str(),name4[k].c_str());
1205     TCanvas *pad = new TCanvas(name,name,500,500);
1206     pad->SetRightMargin(0.10); pad->SetTopMargin(0.10);
1207     hist->GetYaxis()->SetTitle(ytitl4[k].c_str());
1208     hist->GetXaxis()->SetTitle(xtitl4[k].c_str());
1209     hist->GetXaxis()->SetRangeUser(xrnglo4[k],xrnghi4[k]);
1210     hist->GetYaxis()->SetRangeUser(xrnglo4[k],xrnghi4[k]);
1211     hist->SetMarkerStyle(20); hist->SetMarkerSize(0.2);
1212     hist->SetTitle(title4[k].c_str());
1213     hist->GetXaxis()->SetTitleOffset(1.4);
1214     hist->GetYaxis()->SetTitleOffset(1.4);
1215     hist->Draw();
1216     text1->Draw("same");
1217     if (save) {
1218       sprintf (name, "c_%s%s.gif", prefix.c_str(), name2[k].c_str());
1219       pad->Print(name);
1220     }
1221       }
1222     }
1223   }
1224 }
1225 
1226 void testTopology(bool fine) {
1227   HexTopology ht(fine);
1228   std::string ans;
1229   int cellmax = (fine) ? 239 : 132;
1230   int cell(-1), cellmin(0);
1231   std::cout << "Input Cell: ";
1232   std::cin >> cell;
1233   if (cell >= 0) cellmin = cellmax = cell;
1234   for (cell = cellmin; cell <= cellmax; ++cell) {
1235     for (int ext = 0; ext < 3; ++ext) {
1236       std::vector<int> listCell = ht.neighbours(cell,ext);
1237       std::cout << "Extent " << ext << " for Cell " << cell << " has "
1238         << listCell.size() << " cells" << std::endl;
1239       for (unsigned int k=0; k<listCell.size(); ++k) 
1240     std::cout << " " << listCell[k];
1241       std::cout << std::endl;
1242       std::cout << "Want to continue:";
1243       std::cin  >> ans;
1244       if (ans == "q" || ans == "Q" || ans == "s" || ans == "S") break;
1245     }
1246     if (ans == "q" || ans == "Q") break;
1247   }
1248 }
1249 
1250 bool inside(double x, double y) {
1251 
1252   const double size(6.185);
1253   const double tan30deg(0.5773502693);
1254   double absx = std::abs(x);
1255   double absy = std::abs(y);
1256   bool   flag(false);
1257 
1258   if (absx <= size && absy <= (2.*size*tan30deg)) {
1259     if (absy <= size*tan30deg || absx <= (2.*size - absy/tan30deg)) flag = true;
1260   }
1261   return flag;
1262 }
1263 
1264 void testGeometry() {
1265 
1266   HexCellGeometry geomc(false);
1267   for (int k = 0; k < 133; ++k) {
1268     std::pair<double,double> xy = geomc.position(k);
1269     std::cout << "Coarse Cell[" << k << "] " << xy.first << ":" << xy.second
1270           << std::endl;
1271   }
1272 
1273   HexCellGeometry geomf(true);
1274   for (int k = 0; k < 240; ++k) {
1275     std::pair<double,double> xy = geomf.position(k);
1276     std::cout << "Fine Cell[" << k << "] " << xy.first << ":" << xy.second
1277           << std::endl;
1278   }
1279 
1280   HexTBGeometry geomx;
1281   for (int i = 0; i < 7; ++i) {
1282     std::pair<double,double> xy = geomx.position_wafer(i);
1283       std::cout << "Wafer [" << i << "] " << xy.first << ":"  << xy.second 
1284         << std::endl;
1285     }
1286   for (int k = 0; k < 133; ++k) {
1287     std::pair<double,double> xy = geomx.position(0,k,true);
1288     std::cout << "EE Wafer [0] Cell[" << k << "] " << xy.first << ":" 
1289           << xy.second << std::endl;
1290   }
1291   for (int i = 0; i < 7; ++i) {
1292     for (int k = 0; k < 133; ++k) {
1293       std::pair<double,double> xy = geomx.position(i,k,false);
1294       std::cout << "FH Wafer [" << i << "] Cell[" << k << "] " << xy.first 
1295         << ":"  << xy.second << std::endl;
1296     }
1297   }
1298 }
[ 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!