Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:29:20

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //
0003 // Analysis script to compare energy distribution of TB06 data with MC
0004 //
0005 // TB06Analysis.C        Class to run over the tree created by TB06Analysis
0006 //                       within the framewok of CMSSW. This class produces
0007 //                       histogram file for trigger scintillators as well
0008 //                       as for energy in EB, HB and combined energy for
0009 //                       all tracks and tracks which are MIP in ECAL. To
0010 //                       run this, us the 2 following steps:
0011 // TB06Analysis t(infile, outfile, particle, energy, model, corrEB, corrHB);
0012 // t.Loop()
0013 //
0014 // where
0015 //   infile    string    Name of the input ROOT tree file 
0016 //   outfile   string    Name of the output ROOT histogram file
0017 //   particle  int       particle number (0: pi-, 1: pi+, 2: K-, 3: K+,
0018 //                                        4: p, 5: pbar)
0019 //   energy    int       signifies the beam momentum (a value between 0 and
0020 //                       14 to get momenta 2, 3, 4, 5, 6, 7, 8, 9, 20, 30,
0021 //                       50, 100, 150, 200, 300, 350)
0022 //   model     int       signifies Geant4 Physics list + Geant4 version 
0023 //                       (a value between 0 and 13 for the following lists
0024 //                        0 10.0.p2 QGSP-FTFP-BERT-EML                      1
0025 //                        1 10.0.p2 FTFP-BERT-EML                           2
0026 //                        2 10.2.p2 QGSP-FTFP-BERT-EMM                      4
0027 //                        3 10.2.p2 FTFP-BERT-EMM                           8
0028 //                        4 10.2.p2 FTFP-BERT-ATL-EMM                      16
0029 //            5 10.4.0.beta FTFP_BERT_EMM                      32
0030 //            6 10.3.ref08 FTFP_BERT_EMM                       64
0031 //            7 10.3.ref09 FTFP_BERT_EMM                      128
0032 //            8 10.3.ref10 FTFP_BERT_EMM                      256
0033 //            9 10.3.ref11 FTFP_BERT_EMM                      512
0034 //           10 10.3.p03 FTFP_BERT_EMM                       1024
0035 //                       11 10.4.cand01 FTFP_BERT_EMM                    2048
0036 //                       12 10.4.cand02 FTFP_BERT_EMM                    4096
0037 //                       13 10.4 FTFP_BERT_EMM                           8192
0038 //                       14 10.2.p02 FTFP_BERT_EMM 10.0.pre3            16384
0039 //                       15 10.4 FTFP_BERT_EMM 10.0.pre3                32768
0040 //                       16 10.4 FTFP_BERT_EMM VecGeom 10.0.pre3        65536
0041 //                       17 10.4 FTFP_BERT_EMM VecGeom+CLHEP 10.0.pre3 131072
0042 //                       18 10.4 FTFP_BERT_EMM CLHEP                   262144
0043 //                       19 10.4 FTFP_BERT_EMM VecGeom+CLHEP           524288
0044 //                       20 10.4.ref01 FTFP_BERT_EMM                  1048576
0045 //                       21 10.4.ref02 FTFP_BERT_EMM                  2097152
0046 //                       22 10.5.0.beta FTFP_BERT_EMM                 4194304
0047 //                       23 10.4.ref07 FTFP_BERT_EMM                  8388608
0048 //                       24 10.4.ref08 FTFP_BERT_EMM                 16777216
0049 //                       25 10.4.ref08 FTFP_BERT_EMM                 33554432
0050 //                       26 10.5.cand01 FTFP_BERT_EMM                67108864
0051 //                       27 10.4.ref00 FTFP_BERT_EMM                134217728
0052 //                       28 10.4.p03 FTFP_BERT_EMM                  268435456
0053 //                       29 10.5.ref00 FTFP_BERT_EMM                536870912
0054 //                       30 10.5.ref01 FTFP_BERT_EMM               1073741824
0055 //                       31 10.5.ref02 FTFP_BERT_EMM               2147483648
0056 //                       32 10.5.ref03 FTFP_BERT_EMM               4294967296
0057 //                       33 10.5.ref04 FTFP_BERT_EMM               8589935592
0058 //                       34 10.5.ref05 FTFP_BERT_EMM              17179871184
0059 //                       35 10.5.ref06 FTFP_BERT_EMM              34359742368
0060 //   corrEB    double     Correction to noise factor for EB (1.0)
0061 //   corrHB    double     Correction to noise factor for HB (1.0)
0062 //
0063 //  There are several auxiliary functions to make fits/plots:
0064 //
0065 //  DrawTrigger(std::string infile)  
0066 //                        Plots energy distribution in the trigger counters
0067 //                        present in the histogram file *infile*
0068 //  FitTrigger(std::string infile, int type)
0069 //                        Fits Landau function to the energy distribution in
0070 //                        the trigger counter indicated by *type* (0 S1, 1 S2,
0071 //                        2 S3, 3 S4) from the histogram file *infile*
0072 //  FitEBHB(std::string infile, int type)
0073 //                        Fits Crystal Ball (Gaussian) function to energy
0074 //                        deposit in EB (HB) as indicated by *type* : 0 (1)
0075 //                        from the histogram file *infile*
0076 //  DrawHistData(int particle, int energy, std::string dirName, bool rescale)
0077 //                        Makes plots of energy distribution (in ECAL, HCAL
0078 //                        and combined) for particle as given by *particle* 
0079 //                        with beam momentum as given by *energy* from a file 
0080 //                        in directory *dirName*. Possibility of rescaling
0081 //                        energy to match with the mean energy for the given
0082 //                        data point
0083 //  DrawHistDataAll(std::string dirName, bool rescale, bool save)
0084 //                        Plots energy distributions (in ECAL, HCAL, combined)
0085 //                        for all available particles and energies from files
0086 //                        in directory *dirName*. Possibility of rescaling
0087 //                        energy to match with the mean energy for the given
0088 //                        data point and saving the plots as pdf file
0089 //  plotDataMC(int particle, int energy, std::vector<int> models, int type,
0090 //             int rebin, int irtype, std::string prefix, bool approve,
0091 //             bool stat,  double xmin0, double xmax0, bool rescale)
0092 //                        Compares energy measured for all (MIPs in EB) 
0093 //                        tracks as indicated by *type* : 1 (2) for 
0094 //                        *particle* (a value between 0 and 5), *energy*
0095 //                        (a value between 0 and 14), *models* is a list
0096 //                        of models to be included or not. It assumes the 
0097 //                        data files are in directory TB06, MC files are in 
0098 //                        directory *prefix*X (X having a value between 0 
0099 //                        and 32). It produces raw comparison (irtype=1), 
0100 //                        ratio of MC/Data (irtype=2) or both (irtype=3) 
0101 //                        in ths same canvas. It can rebin the default 
0102 //                        histograms using *rebin* and ranges from *xmin0* 
0103 //                        to *xmax0* (if xmax0 is negative, the original 
0104 //                        range of histograms is kept). It can rescale the 
0105 //                        energy in the data histogram to match the mean. 
0106 //                        The legend and stat box are controlled using the 
0107 //                        flags *approve*, *stat*
0108 //  plotDataMCDist(int particle, int energy, std::vector<int> models, int rebin,
0109 //         std::string prefix, bool approve, bool stat, double xmin0,
0110 //                 double xmax0, int save)
0111 //                        Makes comparison plots for approval and save them
0112 //                        as pdf file if needed (*save* > 0). The definition
0113 //                        of other parameters as for plotDataMC
0114 //  plotDataMC(int particle, std::vector<int> models, bool ratio, 
0115 //             std::string dirName, bool approve)
0116 //                        Plots mean response for data and MC or the ratio of
0117 //                        response (MC/Data if *ratio* is true) as a function
0118 //                        of beam momentum for particle *particle*. The mean
0119 //                        values are read out from directory *dirName* and 
0120 //                        models to be plotted are given by *models* (the
0121 //                        usage is described earlier). *approve* decides the
0122 //                        content of legends in the plot
0123 //  plotDataMCResp(std::vector<int> models, std::string dirName, bool approve,
0124 //                 int save)
0125 //                        Plots mean response for data and MC and also the
0126 //                        ratio MC/Data for all available particles. The
0127 //                        legends are controlled by *approve* and saving
0128 //                        the plots as pdf file is done if *save* > 0
0129 ///////////////////////////////////////////////////////////////////////////////
0130 
0131 #include "TCanvas.h"
0132 #include "TChain.h"
0133 #include "TDirectory.h"
0134 #include "TF1.h"
0135 #include "TFile.h"
0136 #include "TFitResult.h"
0137 #include "TFitResultPtr.h"
0138 #include "TGraphAsymmErrors.h"
0139 #include "TH1D.h"
0140 #include "TH2.h"
0141 #include "THStack.h"
0142 #include "TLegend.h"
0143 #include "TMinuit.h"
0144 #include "TMath.h"
0145 #include "TPaveStats.h"
0146 #include "TPaveText.h"
0147 #include "TProfile.h"
0148 #include "TROOT.h"
0149 #include "TStyle.h"
0150 
0151 #include <cstdlib>
0152 #include <fstream>
0153 #include <iomanip>
0154 #include <iostream>
0155 #include <string>
0156 #include <vector>
0157 
0158 // Header file for the classes stored in the TTree if any.
0159 
0160 bool debug_ = false;
0161 
0162 class TB06Analysis {
0163 public :
0164   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
0165   Int_t           fCurrent; //!current Tree number in a TChain
0166 
0167   // Declaration of leaf types
0168   Double_t        eBeam_;
0169   Double_t        etaBeam_;
0170   Double_t        phiBeam_;
0171   Double_t        edepEC_;
0172   Double_t        edepHB_;
0173   Double_t        edepHO_;
0174   Double_t        noiseEC_;
0175   Double_t        noiseHB_;
0176   Double_t        noiseHO_;
0177   Double_t        edepS1_;
0178   Double_t        edepS2_;
0179   Double_t        edepS3_;
0180   Double_t        edepS4_;
0181   Double_t        edepVC_;
0182   Double_t        edepS7_;
0183   Double_t        edepS8_;
0184 
0185   // List of branches
0186   TBranch        *b_eBeam_;    //!
0187   TBranch        *b_etaBeam_;  //!
0188   TBranch        *b_phiBeam_;  //!
0189   TBranch        *b_edepEC_;   //!
0190   TBranch        *b_edepHB_;   //!
0191   TBranch        *b_edepHO_;   //!
0192   TBranch        *b_noiseEC_;  //!
0193   TBranch        *b_noiseHB_;  //!
0194   TBranch        *b_noiseHO_;  //!
0195   TBranch        *b_edepS1_;   //!
0196   TBranch        *b_edepS2_;   //!
0197   TBranch        *b_edepS3_;   //!
0198   TBranch        *b_edepS4_;   //!
0199   TBranch        *b_edepVC_;   //!
0200   TBranch        *b_edepS7_;   //!
0201   TBranch        *b_edepS8_;   //!
0202 
0203   TB06Analysis(std::string infname, std::string outfname, int ipar, int ien,
0204            int model, double corrEB=1, double corrHB=1);
0205   virtual ~TB06Analysis();
0206   virtual Int_t    Cut(Long64_t entry);
0207   virtual Int_t    GetEntry(Long64_t entry);
0208   virtual Long64_t LoadTree(Long64_t entry);
0209   virtual void     Init(TTree *tree);
0210   virtual void     Loop();
0211   virtual Bool_t   Notify();
0212   virtual void     Show(Long64_t entry = -1);
0213   std::string      outName_, modName_, partNam_, partFil_;
0214   double           scaleEB_, scaleHB_, mipCut_, cutS4_, cutVC_, cutS8_;
0215   double           xmin_, xmax_, delx_, corrEB_, corrHB_;
0216   int              model_, nbin_, iene_;
0217 };
0218 
0219 TB06Analysis::TB06Analysis(std::string inName, std::string outName, int ipar,
0220                int ien, int model, double corrEB,
0221                double corrHB) : fChain(0), outName_(outName),
0222                         corrEB_(corrEB), corrHB_(corrHB) {
0223 
0224   TFile      *file = new TFile(inName.c_str());
0225   TDirectory *dir  = (TDirectory*)file->FindObjectAny("testbeam");
0226   TTree      *tree = (TTree*)dir->Get("TB06Sim");
0227   Init(tree);
0228   const int nmodels=36;
0229   double scaleEB[nmodels] = { 1.010,  1.011,  1.011,  1.011,  1.011, 1.011, 1.011, 1.011, 1.011, 1.011, 1.011, 1.011, 1.010, 1.010, 1.015, 1.015, 1.015, 1.015, 1.015, 1.015, 1.010, 1.011, 1.015, 1.015, 1.015, 1.015, 1.015, 1.015, 1.015, 1.015, 1.015, 1.015, 1.015, 1.015, 1.015, 1.015};
0230   double scaleHB[nmodels] = {114.13, 114.29, 106.61, 106.54, 106.33, 106.43, 106.38, 106.50, 106.59, 107.37, 106.47, 107.14, 107.11, 107.11, 105.94, 106.59, 106.64, 106.64, 106.59, 106.64, 107.11, 107.09, 105.89, 105.93, 105.87, 105.93, 105.86, 106.58, 106.65, 105.86, 105.91, 105.91, 105.90, 105.86, 105.78, 105.81};
0231   double cutS4[nmodels]   = {0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028, 0.0028};
0232   double cutVC[nmodels]   = {0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014};
0233   double cutS8[nmodels]   = {0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014, 0.0014};
0234   std::string modelNames[nmodels] = {"10.0.p02 QGSP_FTFP_BERT_EML",
0235                      "10.0.p02 FTFP_BERT_EML",
0236                      "10.2.p02 QGSP_FTFP_BERT_EMM",
0237                      "10.2.p02 FTFP_BERT_EMM",
0238                      "10.2.p02 FTFP_BERT_ATL_EMM",
0239                      "10.4.0.beta FTFP_BERT_EMM",
0240                      "10.3.ref08 FTFP_BERT_EMM",
0241                      "10.3.ref09 FTFP_BERT_EMM",
0242                      "10.3.ref10 FTFP_BERT_EMM",
0243                      "10.3.ref11 FTFP_BERT_EMM",
0244                      "10.3.p03 FTFP_BERT_EMM",
0245                      "10.4.cand01 FTFP_BERT_EMM",
0246                      "10.4.cand02 FTFP_BERT_EMM",
0247                      "10.4 FTFP_BERT_EMM",
0248                      "10.2.p02 FTFP_BERT_EMM 10.0.pre3",
0249                      "10.4 FTFP_BERT_EMM 10.0.pre3",
0250                      "10.4 FTFP_BERT_EMM VecGeom 10.0.pre3",
0251                      "10.4 FTFP_BERT_EMM VecGeom+CLHEP 10.0.pre3",
0252                      "10.4 FTFP_BERT_EMM CLHEP",
0253                      "10.4 FTFP_BERT_EMM VecGeom+CLHEP",
0254                      "10.4.ref01 FTFP_BERT_EMM",
0255                      "10.4.ref02 FTFP_BERT_EMM",
0256                      "10.5.0.beta FTFP_BERT_EMM",
0257                      "10.4.ref07 FTFP_BERT_EMM",
0258                      "10.4.ref08 FTFP_BERT_EMM",
0259                      "10.4.ref08 FTFP_BERT_EMM",
0260                      "10.5.cand01 FTFP_BERT_EMM",
0261                      "10.4.ref00 FTFP_BERT_EMM",
0262                      "10.4.p03 FTFP_BERT_EMM",
0263                      "10.5.ref00 FTFP_BERT_EMM",
0264                      "10.5.ref01 FTFP_BERT_EMM",
0265                      "10.5.ref02 FTFP_BERT_EMM",
0266                      "10.5.ref03 FTFP_BERT_EMM",
0267                      "10.5.ref04 FTFP_BERT_EMM",
0268                      "10.5.ref05 FTFP_BERT_EMM"
0269                      "10.5.ref06 FTFP_BERT_EMM"
0270   };
0271   std::string partsF[7] = {"pi-","pi+","k-","k+","pro+","pro-","e-"};
0272   std::string partsN[7] = {"#pi^{-}","#pi^{+}","K^{-}","K^{+}","p","pbar","e"};
0273   int         iens[16]  = {2,3,4,5,6,7,8,9,20,30,50,100,150,200,300,350};
0274   /*
0275   double      xminh[16] = {-5,-5,-5,-5,-5,-5,-5,-5,5,10,10,10,10,10,10,10};
0276   double      xmaxh[16] = {12,20,20,25,25,25,30,35,50,60,100,200,250,350,450,500};
0277   double      delxs[16] = {.1,.1,.1,.1,.1,.25,.25,.25,.5,.5,1.,1.,2.,2.,2.,2.};
0278   */
0279   double      xminh[16] = {-2,-2,-1,-1,-1,-1,-1,-1,2,6,10,40,70,90,120,140};
0280   double      xmaxh[16] = {6,9,11,14,15,17,18,20,30,40,70,130,170,250,370,420};
0281   double      delxs[16] = {.2,.2,.2,.2,.2,.2,.2,.2,.5,.5,1.,1.,2.,2.,2.,2.};
0282 
0283   if (model < 0 || model >= nmodels)  model_ = 0;
0284   else                                model_ = model;
0285   if (ien < 0 || ien > 15)            ien    = 0;
0286   if (ipar < 0 || ipar > 6)           ipar   = 0;
0287   scaleEB_ = scaleEB[model_];
0288   scaleHB_ = scaleHB[model_];
0289   modName_ = modelNames[model_];
0290   mipCut_  = 1.20;
0291   cutS4_   = cutS4[model_];
0292   cutVC_   = cutVC[model_];
0293   cutS8_   = cutS8[model_];
0294   xmin_    = xminh[ien];
0295   xmax_    = xmaxh[ien];
0296   delx_    = delxs[ien];
0297   nbin_    = (int)((xmax_-xmin_)/delx_);
0298   iene_    = iens[ien];
0299   partFil_ = partsF[ipar];
0300   partNam_ = partsN[ipar];
0301   std::cout << "Model " << model_ << ":" << modName_ << " Particle " << partFil_
0302         << " Scale " << scaleEB_ << ":" << scaleHB_ << " Cuts " << cutS4_
0303         << ":" << cutVC_ << ":" << cutS8_ << std::endl;
0304 }
0305 
0306 TB06Analysis::~TB06Analysis() {
0307   if (!fChain) return;
0308   delete fChain->GetCurrentFile();
0309 }
0310 
0311 Int_t TB06Analysis::GetEntry(Long64_t entry) {
0312   // Read contents of entry.
0313   if (!fChain) return 0;
0314   return fChain->GetEntry(entry);
0315 }
0316 
0317 Long64_t TB06Analysis::LoadTree(Long64_t entry) {
0318   // Set the environment to read one entry
0319   if (!fChain) return -5;
0320   Long64_t centry = fChain->LoadTree(entry);
0321   if (centry < 0) return centry;
0322   if (fChain->GetTreeNumber() != fCurrent) {
0323     fCurrent = fChain->GetTreeNumber();
0324     Notify();
0325   }
0326   return centry;
0327 }
0328 
0329 void TB06Analysis::Init(TTree *tree) {
0330   // The Init() function is called when the selector needs to initialize
0331   // a new tree or chain. Typically here the branch addresses and branch
0332   // pointers of the tree will be set.
0333   // It is normally not necessary to make changes to the generated
0334   // code, but the routine can be extended by the user if needed.
0335   // Init() will be called many times when running on PROOF
0336   // (once per file to be processed).
0337   
0338   // Set branch addresses and branch pointers
0339   if (!tree) return;
0340   fChain = tree;
0341   fCurrent = -1;
0342   fChain->SetMakeClass(1);
0343 
0344   fChain->SetBranchAddress("eBeam_",   &eBeam_,   &b_eBeam_);
0345   fChain->SetBranchAddress("etaBeam_", &etaBeam_, &b_etaBeam_);
0346   fChain->SetBranchAddress("phiBeam_", &phiBeam_, &b_phiBeam_);
0347   fChain->SetBranchAddress("edepEC_",  &edepEC_,  &b_edepEC_);
0348   fChain->SetBranchAddress("edepHB_",  &edepHB_,  &b_edepHB_);
0349   fChain->SetBranchAddress("edepHO_",  &edepHO_,  &b_edepHO_);
0350   fChain->SetBranchAddress("noiseEC_", &noiseEC_, &b_noiseEC_);
0351   fChain->SetBranchAddress("noiseHB_", &noiseHB_, &b_noiseHB_);
0352   fChain->SetBranchAddress("noiseHO_", &noiseHO_, &b_noiseHO_);
0353   fChain->SetBranchAddress("edepS1_",  &edepS1_,  &b_edepS1_);
0354   fChain->SetBranchAddress("edepS2_",  &edepS2_,  &b_edepS2_);
0355   fChain->SetBranchAddress("edepS3_",  &edepS3_,  &b_edepS3_);
0356   fChain->SetBranchAddress("edepS4_",  &edepS4_,  &b_edepS4_);
0357   fChain->SetBranchAddress("edepVC_",  &edepVC_,  &b_edepVC_);
0358   fChain->SetBranchAddress("edepS7_",  &edepS7_,  &b_edepS7_);
0359   fChain->SetBranchAddress("edepS8_",  &edepS8_,  &b_edepS8_);
0360   Notify();
0361 
0362 }
0363 
0364 Bool_t TB06Analysis::Notify() {
0365   // The Notify() function is called when a new file is opened. This
0366   // can be either for a new TTree in a TChain or when when a new TTree
0367   // is started when using PROOF. It is normally not necessary to make changes
0368   // to the generated code, but the routine can be extended by the
0369   // user if needed. The return value is currently not used.
0370 
0371   return kTRUE;
0372 }
0373 
0374 void TB06Analysis::Show(Long64_t entry) {
0375   // Print contents of entry.
0376   // If entry is not specified, print current entry
0377   if (!fChain) return;
0378   fChain->Show(entry);
0379 }
0380 
0381 Int_t TB06Analysis::Cut(Long64_t ) {
0382   // This function may be called from Loop.
0383   // returns  1 if entry is accepted.
0384   // returns -1 otherwise.
0385   return 1;
0386 }
0387 
0388 void TB06Analysis::Loop() {
0389   //   In a ROOT session, you can do:
0390   //      root> .L TB06Analysis.C
0391   //      root> TB06Analysis t
0392   //      root> t.GetEntry(12); // Fill t data members with entry number 12
0393   //      root> t.Show();       // Show values of entry 12
0394   //      root> t.Show(16);     // Read and show values of entry 16
0395   //      root> t.Loop();       // Loop on all entries
0396   //
0397 
0398   //     This is the loop skeleton where:
0399   //    jentry is the global entry number in the chain
0400   //    ientry is the entry number in the current Tree
0401   //  Note that the argument to GetEntry must be:
0402   //    jentry for TChain::GetEntry
0403   //    ientry for TTree::GetEntry and TBranch::GetEntry
0404   //
0405   //       To read only selected branches, Insert statements like:
0406   // METHOD1:
0407   //    fChain->SetBranchStatus("*",0);  // disable all branches
0408   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0409   // METHOD2: replace line
0410   //    fChain->GetEntry(jentry);       //read all branches
0411   //by  b_branchname->GetEntry(ientry); //read only this branch
0412 
0413   if (fChain == 0) return;
0414   TFile *fout = new TFile(outName_.c_str(), "RECREATE");
0415   TH1D  *h_es[7], *h_EC, *h_HB, *h_EN1, *h_EN2;
0416   std::string names[7] = {"h1", "h2", "h3", "h4", "h5", "h6", "h7"};
0417   std::string title[7] = {"Energy in S1", "Energy in S2", "Energy in S3",
0418               "Energy in S4", "Energy in VC", "Energy in S7",
0419               "Energy in S8"};
0420   for (int i=0; i<7; ++i) {
0421     h_es[i] = new TH1D(names[i].c_str(), title[i].c_str(), 1000, 0, 0.02);
0422     h_es[i]->GetXaxis()->SetTitle("Energy (GeV)");
0423     h_es[i]->GetYaxis()->SetTitle("Events");
0424   }
0425   h_EC = new TH1D("EC", "Energy in EB", 2000, 0, 100.0);
0426   h_EC->GetXaxis()->SetTitle("Energy (GeV)");
0427   h_EC->GetYaxis()->SetTitle("Events");
0428   h_HB = new TH1D("HB", "Energy in HB", 2000, 0, 1.0);
0429   h_HB->GetXaxis()->SetTitle("Energy (GeV)");
0430   h_HB->GetYaxis()->SetTitle("Events");
0431   char nameh[100], titlh[100];
0432   sprintf (nameh, "EN1%s%d", partFil_.c_str(), model_);
0433   sprintf (titlh, "Total Energy (%d GeV/c %s)", iene_, partNam_.c_str());
0434   h_EN1 = new TH1D(nameh, titlh, nbin_, xmin_, xmax_);
0435   h_EN1->Sumw2();
0436   h_EN1->GetXaxis()->SetTitle("Energy (GeV)");
0437   h_EN1->GetYaxis()->SetTitle("Events");
0438   sprintf (nameh, "EN2%s%d", partFil_.c_str(), model_);
0439   sprintf (titlh, "Total Energy (MIP in ECAL %d GeV/c %s)", iene_, partNam_.c_str());
0440   h_EN2 = new TH1D(nameh, titlh, nbin_, xmin_, xmax_);
0441   h_EN2->Sumw2();
0442   h_EN2->GetXaxis()->SetTitle("Energy (GeV)");
0443   h_EN2->GetYaxis()->SetTitle("Events");
0444 
0445   Long64_t nentries = fChain->GetEntriesFast();
0446   Long64_t nbytes = 0, nb = 0;
0447   std::cout << "Correction Factors " << corrEB_ << ":" << corrHB_ << " and "
0448         << nentries << " entries" << std::endl;
0449 
0450   for (Long64_t jentry=0; jentry<nentries;jentry++) {
0451     Long64_t ientry = LoadTree(jentry);
0452     if (ientry < 0) break;
0453     nb = fChain->GetEntry(jentry);   nbytes += nb;
0454     // if (Cut(ientry) < 0) continue;
0455     h_es[0]->Fill(edepS1_);
0456     h_es[1]->Fill(edepS2_);
0457     h_es[2]->Fill(edepS3_);
0458     h_es[3]->Fill(edepS4_);
0459     h_es[4]->Fill(edepVC_);
0460     h_es[5]->Fill(edepS7_);
0461     h_es[6]->Fill(edepS8_);
0462     if (debug_) {
0463       std::cout << edepS4_ << ":" << edepVC_ << ":" << edepS8_ << " Cuts " 
0464         << cutS4_ << ":" << cutVC_ << ":" << cutS8_ << " Result " 
0465         << (edepS4_ < cutS4_ && edepVC_ < cutVC_ && edepS8_ < cutS8_) 
0466         << std::endl;
0467       std::cout << edepEC_ << ":" << edepHB_ << " Scale " << scaleEB_ << ":" 
0468         << scaleHB_ << " Noise " << corrEB_*noiseEC_ << ":" 
0469         << corrHB_*noiseHB_ << std::endl;
0470     }
0471     if (edepS4_ < cutS4_ && edepVC_ < cutVC_ && edepS8_ < cutS8_) {
0472       h_EC->Fill(edepEC_);
0473       h_HB->Fill(edepHB_);
0474       double enEB = scaleEB_*(edepEC_+corrEB_*noiseEC_);
0475       double enHB = scaleHB_*edepHB_+corrHB_*noiseHB_;
0476       double eTot = enEB+enHB;
0477       h_EN1->Fill(eTot);
0478       if (enEB < mipCut_) h_EN2->Fill(eTot);
0479     }
0480   }
0481   fout->cd(); fout->Write(); fout->Close();
0482 }
0483 
0484 TCanvas* DrawTrigger(std::string fileName) {
0485   std::string names[7]  = {"h1", "h2", "h3", "h4", "h5", "h6", "h7"};
0486   std::string title[7]  = {"S1", "S2", "S3", "S4", "Veto Counter", "MC1", "MC2"};
0487   int         colors[7] = {2,6,4,1,7,9,3};
0488   int         mtype[7]  = {20,21,22,23,24,33,25};
0489 
0490   TCanvas* pad(0);
0491   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0492   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
0493   gStyle->SetOptTitle(0);
0494   gStyle->SetOptStat(1110);       gStyle->SetOptFit(0);
0495   TFile* file = new TFile(fileName.c_str());
0496   if (file) {  
0497     bool first(true);
0498     double dy  = 0.06;
0499     double yh  = 0.90;
0500     double yh1 = yh-5*dy - 0.01;
0501     pad = new TCanvas("Trigger","TriggerCounter", 700, 500);
0502     pad->SetLogy();
0503     TLegend *legend = new TLegend(0.65, yh1-0.20, 0.90, yh1);
0504     legend->SetFillColor(kWhite);
0505     for (int i=0; i<7; ++i) {
0506       TH1D *hist = (TH1D*)file->FindObjectAny(names[i].c_str());
0507       if (hist) {
0508     hist->SetLineColor(colors[i]);
0509     hist->SetMarkerColor(colors[i]);
0510     hist->SetMarkerStyle(mtype[i]);
0511     hist->GetXaxis()->SetTitle("Energy (GeV)");
0512     hist->GetYaxis()->SetTitle("Events");
0513     hist->GetYaxis()->SetLabelOffset(0.005);
0514     hist->GetYaxis()->SetTitleOffset(1.20);
0515     if (first) hist->Draw("");
0516     else       hist->Draw("sames");
0517     first = false;
0518     pad->Update();
0519     TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
0520     if (st1 != NULL) {
0521       st1->SetFillColor(kWhite);
0522       st1->SetLineColor(colors[i]);
0523       st1->SetTextColor(colors[i]);
0524       st1->SetY1NDC(yh-dy); st1->SetY2NDC(yh);
0525       st1->SetX1NDC(0.70); st1->SetX2NDC(0.90);
0526       yh -= dy;
0527     }
0528     legend->AddEntry(hist,title[i].c_str(),"lp");
0529       }
0530     }
0531     legend->Draw("same");
0532     pad->Update();
0533   }
0534   return pad;
0535 }
0536 
0537 Double_t langaufun(Double_t *x, Double_t *par) {
0538 
0539   //Fit parameters:
0540   //par[0]=Width (scale) parameter of Landau density
0541   //par[1]=Most Probable (MP, location) parameter of Landau density
0542   //par[2]=Total area (integral -inf to inf, normalization constant)
0543   //par[3]=Width (sigma) of convoluted Gaussian function
0544   //
0545   //In the Landau distribution (represented by the CERNLIB approximation), 
0546   //the maximum is located at x=-0.22278298 with the location parameter=0.
0547   //This shift is corrected within this function, so that the actual
0548   //maximum is identical to the MP parameter.
0549 
0550   // Numeric constants
0551   Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
0552   Double_t mpshift  = -0.22278298;       // Landau maximum location
0553   
0554   // Control constants
0555   Double_t np = 100.0;      // number of convolution steps
0556   Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
0557 
0558   // Variables
0559   Double_t xx;
0560   Double_t mpc;
0561   Double_t fland;
0562   Double_t sum = 0.0;
0563   Double_t xlow,xupp;
0564   Double_t step;
0565   Double_t i;
0566   Double_t par0=0.2;
0567 
0568   // MP shift correction
0569   mpc = par[1] - mpshift * par0 * par[1]; 
0570 
0571   // Range of convolution integral
0572   xlow = x[0] - sc * par[3];
0573   xupp = x[0] + sc * par[3];
0574 
0575   step = (xupp-xlow) / np;
0576 
0577   // Convolution integral of Landau and Gaussian by sum
0578   for(i=1.0; i<=np/2; i++) {
0579     xx = xlow + (i-.5) * step;
0580     fland = TMath::Landau(xx,mpc,par0*par[1], kTRUE); // / par[0];
0581     sum += fland * TMath::Gaus(x[0],xx,par[3]);
0582 
0583     xx = xupp - (i-.5) * step;
0584     fland = TMath::Landau(xx,mpc,par0*par[1], kTRUE); // / par[0];
0585     sum += fland * TMath::Gaus(x[0],xx,par[3]);
0586   }
0587 
0588   return (par[2] * step * sum * invsq2pi / par[3]);
0589 }
0590 
0591 Double_t crystalfun(Double_t *x, Double_t *par) {
0592 
0593   // Variables
0594   Double_t xx    = x[0] - par[1];
0595   Double_t sigma = par[2];
0596   Double_t an    = par[3]; // 6.45
0597   Double_t alpha = par[4]; // 0.91
0598   
0599   double crystal = 0;
0600   if (xx > -alpha*sigma) {
0601     crystal = par[0]*std::exp(-0.5*(xx/sigma)*(xx/sigma));
0602   } else {
0603     Double_t den = (1.0-(alpha/an)*(xx/sigma) - (alpha*alpha)/an);
0604     if (den > 0) {
0605       Double_t den1 = an*std::log(den);
0606       crystal       = par[0]*std::exp(-0.5*alpha*alpha)/std::exp(den1);
0607     }
0608   }
0609   return crystal;
0610 }
0611 
0612 std::pair<TF1*,TFitResultPtr> functionFit(TH1D *his, double *fitrange, 
0613                       double *startvalues, 
0614                       double *parlimitslo, 
0615                       double *parlimitshi, int mode) {
0616 
0617   char FunName[100];
0618   std::string fname("LanGau");
0619   if (mode < 0) fname = "CrysBall";
0620   sprintf(FunName,"%s_%s",fname.c_str(), his->GetName());
0621   
0622   TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
0623   if (ffitold) delete ffitold;
0624 
0625   TF1 *ffit;
0626   int npar=4;
0627   if (mode >=0) {
0628     ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],npar);
0629   } else {
0630     npar = 5;
0631     ffit = new TF1(FunName,crystalfun,fitrange[0],fitrange[1],npar);
0632   }
0633   ffit->SetParameters(startvalues);
0634   if (mode >=0) {
0635     ffit->SetParNames("Width","MP","Area","GSigma");
0636   } else {
0637     ffit->SetParNames("Area","Mean","Width","AN","Alpha");
0638   }
0639   if (mode == 0) {
0640     for (int i=0; i<npar; i++) 
0641       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0642   }
0643   // fit within specified range, use ParLimits, do not plot
0644   TFitResultPtr fptr = his->Fit(FunName,"wwqRB0S"); 
0645   return std::pair<TF1*,TFitResultPtr>(ffit,fptr); 
0646 }
0647 
0648 TCanvas* FitTrigger(std::string fileName, int type) {
0649   std::string names[7]  = {"h1", "h2", "h3", "h4", "h5", "h6", "h7"};
0650   std::string title[7]  = {"S1", "S2", "S3", "S4", "Veto Counter", "MC1", "MC2"};
0651   int         colors[7] = {2,6,4,1,7,9,3};
0652   int         mtype[7]  = {20,21,22,23,24,33,25};
0653 
0654   TCanvas* pad(0);
0655   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0656   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
0657   gStyle->SetOptTitle(0);
0658   gStyle->SetOptStat(1110);       gStyle->SetOptFit(10);
0659   TFile* file = new TFile(fileName.c_str());
0660   if (file) {
0661     double dy  = 0.15;
0662     double yh  = 0.90;
0663     double yh1 = yh-dy-0.01;
0664     TH1D *hist = (TH1D*)file->FindObjectAny(names[type].c_str());
0665     if (hist) {
0666       char name[100];
0667       sprintf(name,"TriggerFit%d",type);
0668       pad = new TCanvas(name,name,700,500);
0669       TLegend *legend = new TLegend(0.65, yh1-0.05, 0.90, yh1);
0670       legend->SetFillColor(kWhite);
0671       hist->SetLineColor(colors[type]);
0672       hist->SetMarkerColor(colors[type]);
0673       hist->SetMarkerStyle(mtype[type]);
0674       hist->GetXaxis()->SetTitle("Energy (GeV)");
0675       hist->GetYaxis()->SetTitle("Events");
0676       hist->GetYaxis()->SetLabelOffset(0.005);
0677       hist->GetYaxis()->SetTitleOffset(1.20);
0678       double LowEdge(0.001);
0679       double HighEdge(0.005);
0680       TF1 *Fitfun = new TF1(name,"landau",LowEdge,HighEdge);
0681       TFitResultPtr Fit = hist->Fit(Fitfun, "+0wwqR");
0682       hist->GetXaxis()->SetRangeUser(0, 0.005);
0683       hist->Draw("");
0684       Fitfun->Draw("same");
0685       //      hist->Fit("landau", "wwqRS", "", LowEdge, HighEdge);
0686       pad->Update();
0687       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
0688       if (st1 != NULL) {
0689     st1->SetFillColor(kWhite);
0690     st1->SetLineColor(colors[type]);
0691     st1->SetTextColor(colors[type]);
0692     st1->SetY1NDC(yh-dy); st1->SetY2NDC(yh);
0693     st1->SetX1NDC(0.60); st1->SetX2NDC(0.90);
0694       }
0695       legend->AddEntry(hist,title[type].c_str(),"lp");
0696       legend->Draw("same");
0697       pad->Update();
0698     }
0699   }
0700   return pad;
0701 }
0702 
0703 TCanvas* FitEBHB(std::string fileName, int type, int model=0, bool save=false) {
0704   std::string names[2]  = {"EC", "HB"};
0705   std::string title[2]  = {"Energy in EB", "Energy in HB"};
0706   int         colors[6] = {2,6,4,1,7,9};
0707   int         mtype[6]  = {20,21,22,23,24,33};
0708 
0709   TCanvas* pad(0);
0710   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0711   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
0712   gStyle->SetOptTitle(0);
0713   gStyle->SetOptStat(1110);       gStyle->SetOptFit(10);
0714   char name[100];
0715   sprintf (name, "model%d/%s", model, fileName.c_str());
0716   TFile* file = new TFile(name);
0717   if (file) {
0718     double dy  = 0.20;
0719     double yh  = 0.90;
0720     double yh1 = yh-dy-0.01;
0721     double xh  = (type == 0) ? 0.40 : 0.90;
0722     TH1D *hist = (TH1D*)file->FindObjectAny(names[type].c_str());
0723     if (hist) {
0724       sprintf(name,"%sFit%d",names[type].c_str(),model);
0725       pad = new TCanvas(name,name,700,500);
0726       hist->SetLineColor(colors[type]);
0727       hist->SetLineWidth(2);
0728       hist->SetMarkerColor(colors[type]);
0729       hist->SetMarkerStyle(mtype[type]);
0730       hist->GetXaxis()->SetTitle("Energy (GeV)");
0731       hist->GetYaxis()->SetTitle("Events");
0732       hist->GetYaxis()->SetLabelOffset(0.005);
0733       hist->GetYaxis()->SetTitleOffset(1.20);
0734       double LowEdge  = hist->GetMean() - 0.8*hist->GetRMS();
0735       double HighEdge = hist->GetMean() + 2.*hist->GetRMS();
0736       TFitResultPtr Fit = hist->Fit("gaus", "+0wwqRS", "", LowEdge, HighEdge);
0737       std::cout << Fit->Value(0) << " 1 " << Fit->Value(1) << " 2 " << Fit->Value(2) << std::endl;
0738       double startvalues[5], fitrange[2], lowValue[5], highValue[5];
0739       startvalues[0] = Fit->Value(0); lowValue[0] = 0.0; highValue[0] = 10.0*startvalues[0];
0740       startvalues[1] = Fit->Value(1); lowValue[1] = 0.5*startvalues[1]; highValue[1] = 2.*startvalues[1];
0741       startvalues[2] = Fit->Value(2); lowValue[2] = 0.5*startvalues[2]; highValue[2] = 2.*startvalues[2];
0742       startvalues[3] = 6.45;              lowValue[3] = 0.5*startvalues[3]; highValue[3] = 2.*startvalues[3];
0743       startvalues[4] = 0.91;              lowValue[4] = 0.5*startvalues[4]; highValue[4] = 2.*startvalues[4];
0744       TF1 *Fitfun;
0745       if (type == 0) {
0746     fitrange[0] = startvalues[1]-10.*startvalues[2];
0747     fitrange[1] = startvalues[1]+2.*startvalues[2];
0748     std::pair<TF1*,TFitResultPtr> ffit = functionFit(hist, fitrange, startvalues, lowValue,highValue,-1);
0749     Fitfun = ffit.first;
0750     Fit    = ffit.second;
0751       } else {
0752     fitrange[0] = startvalues[1]-10.*startvalues[2];
0753     fitrange[1] = startvalues[1]+10.*startvalues[2];
0754     Fitfun = new TF1(name,"gaus",fitrange[0],fitrange[1]);
0755     Fit = hist->Fit(Fitfun, "+0wwqRS");
0756       }
0757       std::cout << LowEdge << " " << startvalues[1]-2.*startvalues[3] << " " 
0758         << HighEdge << " " << fitrange[1] << " Fit " << Fit->Value(1) 
0759         << ":" << (50.0/Fit->Value(1)) << std::endl;
0760       double low  = fitrange[0]-2.*startvalues[2]; 
0761       double high = fitrange[1]+2.*startvalues[2];
0762       hist->GetXaxis()->SetRangeUser(low, high);
0763       hist->Draw("");
0764       Fitfun->Draw("same");
0765       pad->Update();
0766       
0767       TPaveStats* st1 = (TPaveStats*)hist->GetListOfFunctions()->FindObject("stats");
0768       if (st1 != NULL) {
0769     st1->SetFillColor(kWhite);
0770     st1->SetLineColor(colors[type]);
0771     st1->SetTextColor(colors[type]);
0772     st1->SetY1NDC(yh-dy);  st1->SetY2NDC(yh);
0773     st1->SetX1NDC(xh-0.3); st1->SetX2NDC(xh);
0774       }
0775       TLegend *legend = new TLegend(xh-0.25, yh1-0.05, xh, yh1);
0776       legend->SetFillColor(kWhite);
0777       legend->AddEntry(hist,title[type].c_str(),"lp");
0778       legend->Draw("same");
0779       pad->Update();
0780       if (save) {
0781     sprintf (name, "%s.pdf", pad->GetName());
0782     pad->Print(name);
0783       }
0784     }
0785   }
0786   return pad;
0787 }
0788 
0789 double GetScaleFactor(int type, int ie, std::string dirName="RespAll") {
0790   char infil1[200], infil2[200];
0791   std::string partsF[6] = {"pim","pip","km","kp","prop","prom"};
0792   int         iens[16]  = {2,3,4,5,6,7,8,9,20,30,50,100,150,200,300,350};
0793   if (dirName == "") {
0794     sprintf(infil1,"%s.txt",  partsF[type].c_str());
0795     sprintf(infil2,"%sn.txt", partsF[type].c_str());
0796   } else {
0797     sprintf(infil1,"%s/%s.txt",  dirName.c_str(), partsF[type].c_str());
0798     sprintf(infil2,"%s/%sn.txt", dirName.c_str(), partsF[type].c_str());
0799   }
0800   double pbeam, resp, reso, scale(1.0), resp1(0), resp2(0);
0801   bool   ok(false);
0802   std::ifstream fInput1(infil1);
0803   if (!fInput1.good()) {
0804     std::cout << "Cannot open file " << infil1 << std::endl;
0805   } else {
0806     while (1) {
0807       fInput1 >> pbeam >> resp >> reso;
0808       if (!fInput1.good()) break;
0809       int ip = (int)(pbeam+0.01);
0810       if (ip == iens[ie]) {
0811     resp1 = resp; ok = true; break;
0812       }
0813     }
0814     fInput1.close();
0815   }
0816   if (ok) {
0817     std::ifstream fInput2(infil2);
0818     ok = false;
0819     if (!fInput2.good()) {
0820       std::cout << "Cannot open file " << infil2 << std::endl;
0821     } else {
0822       while (1) {
0823     fInput2 >> pbeam >> resp >> reso;
0824     if (!fInput2.good()) break;
0825     int ip = (int)(pbeam+0.01);
0826     if (ip == iens[ie]) {
0827       resp2 = resp; ok = true; break;
0828     }
0829       }
0830       fInput2.close();
0831     }
0832   }
0833   if (ok) scale = resp1/resp2;
0834   if (debug_) 
0835     std::cout << type << ":" << ie << ":" << iens[ie] << " response " << resp1 
0836           << ":" << resp2 << ":" << scale << std::endl;
0837   return scale;
0838 }
0839 
0840 bool FillHistData(char* infile, TH1D* h_eEB, TH1D* h_eHB, TH1D* h_eTot,
0841           double scale=1.0) {
0842 
0843   bool flag(false);
0844   std::ifstream fInput(infile);
0845   if (!fInput.good()) {
0846     std::cout << "Cannot open file " << infile << std::endl;
0847   } else {
0848     char buffer [1024];
0849     unsigned int all(1), good(0);
0850     fInput.getline(buffer, 1024);
0851     if (debug_) std::cout << buffer << std::endl;
0852     double eEB, eHB, eHO;
0853     double eEBl(9999), eEBh(-9999), eHBl(9999), eHBh(-9999), eTotl(9999), eToth(-9999);
0854     while (1) {
0855       fInput >> eEB >> eHB >> eHO;
0856       if (!fInput.good()) break;
0857       all++; good++;
0858       double eTot = eEB+eHB;
0859       if (h_eEB  != 0) h_eEB->Fill(scale*eEB);
0860       if (h_eHB  != 0) h_eHB->Fill(scale*eHB);
0861       if (h_eTot != 0) h_eTot->Fill(scale*eTot);
0862       if (eEB  < eEBl)  eEBl  = eEB;  if (eEB  > eEBh)  eEBh  = eEB;
0863       if (eHB  < eHBl)  eHBl  = eHB;  if (eHB  > eHBh)  eHBh  = eHB;
0864       if (eTot < eTotl) eTotl = eTot; if (eTot > eToth) eToth = eTot;
0865     }
0866     if (debug_) {
0867       std::cout << "Reads " << all << " (" << good << ") records from "
0868         << infile << std::endl << "Minimum/maximum for EB " << eEBl
0869         << ":" << eEBh << "  HB " << eHBl << ":" << eHBh << "  Total "
0870         << eTotl << ":" << eToth << std::endl;
0871     }
0872     fInput.close();
0873     flag = (good>0);
0874   }
0875   return flag;
0876 }
0877 
0878 std::vector<TCanvas*> DrawHistData(int type, int ie, std::string dirName="TB06",
0879                    bool rescale=false) {
0880   std::string partsF[6] = {"pi-","pi+","k-","k+","pro+","pro-"};
0881   std::string partsN[6] = {"#pi^{-}","#pi^{+}","K^{-}","K^{+}","p","pbar"};
0882   std::string names[3]  = {"ECAL", "HCAL", "Total"};
0883   int         iens[16]  = {2,3,4,5,6,7,8,9,20,30,50,100,150,200,300,350};
0884   double      xmine(-2);
0885   double      xmaxe[16] = {3,4,4,6,6,8,8,10,45,45,100,200,200,350,350,500};
0886   /*
0887   double      xminh(-5);
0888   double      xmaxh[16] = {12,20,20,25,25,25,30,35,50,60,100,200,250,350,450,500};
0889   double      delxs[16] = {.1,.1,.1,.1,.1,.25,.25,.25,.5,.5,1.,1.,2.,2.,2.,2.};
0890   */
0891   double      xminh[16] = {-2,-2,-1,-1,-1,-1,-1,-1,2,6,10,40,70,90,120,140};
0892   double      xmaxh[16] = {6,9,11,14,15,17,18,20,30,40,70,130,170,250,370,420};
0893   double      delxs[16] = {.2,.2,.2,.2,.2,.2,.2,.2,.5,.5,1.,1.,2.,2.,2.,2.};
0894 
0895   std::vector<TCanvas*> pads;
0896   char infile[200], title[50], name[50];
0897   if (dirName == "") sprintf(infile,"%s%d.txt", partsF[type].c_str(), iens[ie]);
0898   else               sprintf(infile,"%s/%s%d.txt", dirName.c_str(), partsF[type].c_str(), iens[ie]);
0899   sprintf (title, "%d GeV %s", iens[ie], partsN[type].c_str());
0900   TH1D* hist[3];
0901   double xmin(xmine), xmax, delx;
0902   int    nbin(0);
0903   xmax = xmaxe[ie];
0904   delx = delxs[ie];
0905   nbin = (int)((xmax-xmin)/delx);
0906   sprintf (name, "%s%s%d", names[0].c_str(), partsF[type].c_str(), iens[ie]);
0907   hist[0] = new TH1D(name, title, nbin, xmin, xmax);
0908   hist[0]->Sumw2();
0909   xmin = xminh[ie];
0910   xmax = xmaxh[ie];
0911   nbin = (int)((xmax-xmin)/delx);
0912   sprintf (name, "%s%s%d", names[1].c_str(), partsF[type].c_str(), iens[ie]);
0913   hist[1] = new TH1D(name, title, nbin, xmin, xmax);
0914   hist[1]->Sumw2();
0915   sprintf (name, "%s%s%d", names[2].c_str(), partsF[type].c_str(), iens[ie]);
0916   hist[2] = new TH1D(name, title, nbin, xmin, xmax);
0917   hist[2]->Sumw2();
0918   double scale = rescale ? GetScaleFactor(type, ie) : 1.0;
0919   bool ok = FillHistData(infile, hist[0], hist[1], hist[2], scale);
0920 
0921   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
0922   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
0923   gStyle->SetOptTitle(0);
0924   gStyle->SetOptStat(1110);       gStyle->SetOptFit(0);
0925   if (ok) {
0926     double mean[3], rms[3];
0927     for (int k=0; k<3; ++k) {
0928       sprintf(name, "%s%s%d", names[k].c_str(), partsF[type].c_str(), iens[ie]);
0929       TCanvas* pad = new TCanvas(name, name, 700, 500);
0930       pad->SetRightMargin(0.10); pad->SetTopMargin(0.10);
0931       sprintf(name, "%s Energy (GeV)", names[k].c_str());
0932       hist[k]->GetXaxis()->SetTitle(name); 
0933       hist[k]->GetYaxis()->SetTitle("Events");
0934       hist[k]->GetYaxis()->SetTitleOffset(1.2);
0935       hist[k]->Draw(); pad->Update();
0936       TPaveStats* st1 = (TPaveStats*)hist[k]->GetListOfFunctions()->FindObject("stats");
0937       if (st1 != NULL) {
0938     st1->SetFillColor(0);
0939     st1->SetY1NDC(0.75); st1->SetY2NDC(0.90);
0940     st1->SetX1NDC(0.70); st1->SetX2NDC(0.90);
0941       }
0942       mean[k] = hist[k]->GetMean(); rms[k] = hist[k]->GetRMS();
0943       TPaveText *text = new TPaveText(0.70,0.70,0.90,0.745,"blNDC");
0944       text->SetFillColor(0); text->AddText(title);   text->Draw("same");
0945       pad->Modified();       pad->Update();
0946       pads.push_back(pad);
0947     }
0948     std::cout << partsF[type] << " " << iens[ie] << " GeV";
0949     for (int k=0; k<3; ++k) std::cout << " " << mean[k] << " " << rms[k];
0950     std::cout << " " << mean[2]/iens[ie] << " " << rms[2]/mean[2] << std::endl;
0951   }
0952   return pads;
0953 }
0954 
0955 void DrawHistDataAll(std::string dirName="TB06", bool rescale=false,
0956              bool save=false) {
0957   char filename[100];
0958   for (int k1=0; k1<6; ++k1) {
0959     for (int k2=0; k2<16; ++k2) {
0960       std::vector<TCanvas*> pads = DrawHistData(k1,k2,dirName,rescale);
0961       if (save) {
0962     for (unsigned int k=0; k<pads.size(); ++k) {
0963       sprintf (filename, "%s.pdf", pads[k]->GetName());
0964       pads[k]->Print(filename);
0965     }
0966       }
0967     }
0968   }
0969 }
0970 
0971 
0972 TCanvas* plotDataMC(int ipar, int ien, std::vector<int> models, int type,
0973             int rebin=1, int irtype=1, std::string prefix="model", 
0974             bool approve=false, bool stat=true, double xmin0=-1, 
0975             double xmax0=-1, bool rescale=true) {
0976   static const int nmodels=36;
0977   std::string modelNames[nmodels] = {"G4 10.0.p02 QGSP_FTFP_BERT_EML",
0978                      "G4 10.0.p02 FTFP_BERT_EML",
0979                      "G4 10.2.p02 QGSP_FTFP_BERT_EMM",
0980                      "G4 10.2.p02 FTFP_BERT_EMM",
0981                      "G4 10.2.p02 FTFP_BERT_ATL_EMM",
0982                      "G4 10.4.beta FTFP_BERT_EMM",
0983                      "G4 10.3.ref08 FTFP_BERT_EMM",
0984                      "G4 10.3.ref09 FTFP_BERT_EMM",
0985                      "G4 10.3.ref10 FTFP_BERT_EMM",
0986                      "G4 10.3.ref11 FTFP_BERT_EMM",
0987                      "G4 10.3.p03 FTFP_BERT_EMM",
0988                      "G4 10.4.cand01 FTFP_BERT_EMM",
0989                      "G4 10.4.cand02 FTFP_BERT_EMM",
0990                      "G4 10.4 FTFP_BERT_EMM",
0991                      "G4 10.2.p02 FTFP_BERT_EMM 10.0.pre3",
0992                      "G4 10.4 FTFP_BERT_EMM 10.0.pre3",
0993                      "G4 10.4 FTFP_BERT_EMM VecGeom 10.0.pre3",
0994                      "G4 10.4 FTFP_BERT_EMM VecGeom+CLHEP 10.0.pre3",
0995                      "G4 10.4 FTFP_BERT_EMM CLHEP",
0996                      "G4 10.4 FTFP_BERT_EMM VecGeom+CLHEP",
0997                      "G4 10.4.ref01 FTFP_BERT_EMM",
0998                      "G4 10.4.ref02 FTFP_BERT_EMM",
0999                      "G4 10.5.0.beta FTFP_BERT_EMM",
1000                      "G4 10.4.ref07 FTFP_BERT_EMM",
1001                      "G4 10.4.ref08 FTFP_BERT_EMM",
1002                      "G4 10.4.ref09 FTFP_BERT_EMM",
1003                      "10.5.cand01 FTFP_BERT_EMM",
1004                      "10.4.ref00 FTFP_BERT_EMM",
1005                      "10.4.p03 FTFP_BERT_EMM",
1006                      "10.5.ref00 FTFP_BERT_EMM",
1007                      "10.5.ref01 FTFP_BERT_EMM",
1008                      "10.5.ref02 FTFP_BERT_EMM",
1009                      "10.5.ref03 FTFP_BERT_EMM",
1010                      "10.5.ref04 FTFP_BERT_EMM",
1011                      "10.5.ref05 FTFP_BERT_EMM",
1012                      "10.5.ref06 FTFP_BERT_EMM"
1013   };
1014   std::string partsF[6] = {"pi-","pi+","k-","k+","pro+","pro-"};
1015   std::string partsM[6] = {"pim","pip","km","kp","prop","prom"};
1016   std::string partsN[6] = {"#pi^{-}","#pi^{+}","K^{-}","K^{+}","proton","antiproton"};
1017   int         iens[16]  = {2,3,4,5,6,7,8,9,20,30,50,100,150,200,300,350};
1018   int         types[16] = {0,0,0,0,0,0,0,0, 1, 1, 1,  1,  1,  2,  2,  2};
1019   int         colors[nmodels]= { 2, 7, 6, 4, 9, 8,40, 7, 6, 9,46,48,30, 2, 6,
1020                  2, 7, 4, 2, 6, 4, 1, 7, 6, 7, 9, 8, 6, 2, 9,
1021                  1, 7, 7, 4, 7, 4};
1022   int         mtype[nmodels] = {21,22,23,24,25,26,27,22,23,25,29,33,42,21,21,
1023                 22,23,24,21,22,23,24,21,22,23,33,26,22,21,23,
1024                 26,21,21,26,21,26};
1025   int         colorD(1), mtypeD(20);
1026   std::string titlty[2] = {"All events", "MIP in ECAL"};
1027 
1028   if (debug_) {
1029     std::cout << "plotDataMC " << ipar << ", " << ien << ", " << models.size()
1030           << " Models:";
1031     for (unsigned int k=0; k<models.size(); ++k)
1032       std::cout << " " << models[k] << ",";
1033     std::cout << " Type... " << type << ", " << rebin << ", " << irtype << ", " 
1034           << prefix << ", " << approve << ", " << xmin0 << ", " << xmax0 
1035           << ", " << rescale << std::endl;
1036   }
1037   TCanvas* pad(0);
1038   gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
1039   gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
1040   gStyle->SetOptTitle(0);         gStyle->SetOptFit(10);
1041   if (stat)   gStyle->SetOptStat(1110);
1042   else        gStyle->SetOptStat(0);
1043   if (type != 2) type = 1;
1044   std::vector<int>   icol;
1045   std::vector<TH1D*> hists;
1046   if (irtype < 0 || irtype > 3) irtype = 1;
1047 
1048   //Get some model information
1049   char     infile[120], title[100], name[120];
1050   double   xmin(-5.0), xmax(500.0), ymax(0);
1051   int      nbin(0), nmod(0);
1052   for (unsigned int k=0; k<models.size(); ++k) {
1053     if (models[k] <= nmodels) {
1054       int i = models[k];
1055       sprintf (infile,"%s%d/%s%d.root", prefix.c_str(), i, partsM[ipar].c_str(),
1056            iens[ien]);
1057       TFile *file = new TFile(infile);
1058       if (debug_) 
1059     std::cout << "File " << i << " " << infile << "  " << file << std::endl;
1060       if (file) {
1061     sprintf (name,"EN%d%s%d",type,partsF[ipar].c_str(),i);
1062     TH1D* h1 = (TH1D*)file->FindObjectAny(name);
1063     if (debug_) 
1064       std::cout << "Hist " << i << "  " << name << "   " << h1 << std::endl;
1065     if (h1) {
1066       nmod++;
1067       if (nbin == 0) {
1068         nbin = h1->GetNbinsX();
1069         xmin = h1->GetBinLowEdge(1);
1070         xmax = h1->GetBinLowEdge(nbin) + h1->GetBinWidth(nbin);
1071       }
1072     }
1073     file->Close();
1074       }
1075     }
1076   }
1077   if (xmax0 > 0) {
1078     nbin = (int)((xmax0-xmin0)/((xmax-xmin)*rebin));
1079     xmax = xmax0;
1080     xmin = xmin0;
1081   } else {
1082     nbin /= rebin;
1083   }
1084 
1085   //Data first
1086   TH1D *hist(0);
1087   sprintf (infile, "TB06/%s%d.txt", partsF[ipar].c_str(), iens[ien]);
1088   sprintf (title, "%d GeV %s", iens[ien], partsN[ipar].c_str());
1089   sprintf (name, "Data%d%d%s%d", type, irtype, partsF[ipar].c_str(), iens[ien]);
1090   double scalef = rescale ? GetScaleFactor(ipar, ien) : 1.0;
1091   std::ifstream fInput(infile);
1092   if (!fInput.good()) {
1093     std::cout << "Cannot open file " << infile << std::endl;
1094   } else {
1095     hist = new TH1D(name, title, nbin, xmin, xmax);
1096     hist->Sumw2();
1097     hist->GetXaxis()->SetTitle("Energy (GeV)");
1098     hist->GetYaxis()->SetTitle("Events");
1099     hist->GetYaxis()->SetTitleOffset(1.3);
1100     hist->SetMarkerColor(colorD);
1101     hist->SetMarkerStyle(mtypeD);
1102     hist->SetMarkerSize(1.2);
1103     hist->SetLineColor(colorD);
1104     hist->SetLineWidth(2);
1105     char buffer [1024];
1106     fInput.getline(buffer, 1024);
1107     double eEB, eHB, eHO;
1108     while (1) {
1109       fInput >> eEB >> eHB >> eHO;
1110       if (!fInput.good()) break;
1111       eEB *= scalef; eHB *= scalef;
1112       double eTot = (eEB+eHB);
1113       if      (type == 1) hist->Fill(eTot);
1114       else if (eEB < 1.2) hist->Fill(eTot);
1115     }
1116     fInput.close();
1117     for (int k=1; k<=nbin; ++k) {
1118       if (hist->GetBinContent(k)+hist->GetBinError(k) > ymax) {
1119     ymax = hist->GetBinContent(k)+hist->GetBinError(k);
1120       }
1121     }
1122     hists.push_back(hist);
1123     icol.push_back(colorD);
1124   }
1125 
1126   //Now loop of MC files
1127   if (hist) {
1128     double total(0);
1129     double legsize(0.04);
1130     for (int k=1; k<=nbin; ++k) total += hist->GetBinContent(k);
1131     double yext(0);
1132     if (irtype == 1) {
1133       yext = 500;
1134       sprintf (name,"EN%d%s%dH",type,partsF[ipar].c_str(),iens[ien]);
1135     } else if (irtype == 2) {
1136       yext = 300;
1137       sprintf (name,"EN%d%s%dR",type,partsF[ipar].c_str(),iens[ien]);
1138     } else {
1139       yext = 700;
1140       sprintf (name,"EN%d%s%d",type,partsF[ipar].c_str(),iens[ien]);
1141     }
1142     pad = new TCanvas(name, name, 700, yext);
1143     double dy  = 0.08;
1144     double xh(0.89), yh(0.89);
1145     double yh1 = yh-dy*(nmod+1)- 0.02;
1146     TLegend *legend(0);
1147     if (!stat) {
1148       if (approve) legend = new TLegend(xh-0.5, yh-legsize*(nmod+1)-0.01, xh, yh-0.01);
1149       else legend = new TLegend(xh-0.35, yh-legsize*(nmod+1)-0.01, xh, yh-0.01);
1150     } else if (types[ien] == 2) {
1151       legend = new TLegend(xh-0.35, yh-legsize*(nmod+1)-0.01, xh, yh-0.01);
1152       xh     = 0.4;
1153     } else if (types[ien] == 1) {
1154       legend = new TLegend(0.11, yh-legsize*(nmod+1), 0.45, yh-0.01);
1155     } else {
1156       legend = new TLegend(xh-0.35, yh1-legsize*(nmod+1), xh, yh1);
1157     }
1158     legend->SetFillColor(kWhite); legend->SetBorderSize(0); 
1159     legend->SetMargin(0.2);
1160 
1161     sprintf (name,"ENStack1%d%s%d",type,partsF[ipar].c_str(),iens[ien]);
1162     THStack *Hs = new THStack(name,"");
1163     Hs->Add(hist,"pe sames");
1164     sprintf (title, "%d GeV %s (%s)", iens[ien], partsN[ipar].c_str(), titlty[type-1].c_str());
1165     legend->AddEntry(hist,title,"lp");
1166 
1167     for (unsigned int k=0; k<models.size(); ++k) {
1168       if (models[k] <= nmodels) {
1169     int i = models[k];
1170     sprintf (infile,"%s%d/%s%d.root",prefix.c_str(),i,partsM[ipar].c_str(),iens[ien]);
1171     TFile *file = new TFile(infile);
1172     if (file) {
1173       sprintf (name,"EN%d%s%d",type,partsF[ipar].c_str(),i);
1174       TH1D* h1 = (TH1D*)file->FindObjectAny(name);
1175       if (h1) {
1176         sprintf (name,"EN%d%s%d1",type,partsF[ipar].c_str(),i);
1177         TH1D *h2 = new TH1D(name, "", nbin, xmin, xmax);
1178         h2->Sumw2();
1179         double totm(0);
1180         for (int k=1; k<=h1->GetNbinsX(); ++k) totm += h1->GetBinContent(k);
1181         double scale = total/totm;
1182         int ibin = (int)((xmin-h1->GetBinLowEdge(1))/h1->GetBinWidth(1));
1183         for (int k=1; k<=nbin; ++k) {
1184           double cont(0);
1185           for (int k1=0; k1<rebin; ++k1) {
1186         ++ibin;
1187         double cv = h1->GetBinContent(ibin);
1188         cont += scale*cv;
1189           }
1190           h2->SetBinContent(k,cont);
1191           if (cont > ymax) ymax = cont;
1192           
1193         }
1194         h2->SetMarkerColor(colors[i]);
1195         h2->SetMarkerStyle(mtype[i]);
1196         h2->SetMarkerSize(1.2);
1197         h2->SetLineColor(colors[i]);
1198         h2->SetLineWidth(2);
1199         Hs->Add(h2,"hist sames");
1200         legend->AddEntry(h2,modelNames[i].c_str(),"lp");
1201         icol.push_back(colors[i]);
1202         hists.push_back(h2);
1203       }
1204     }
1205       }
1206     }
1207 
1208     int    imax = (ymax > 100) ? (int)(0.01*ymax) : (int)(0.1*ymax);
1209     double ymx  = (ymax > 100) ? 100*(imax+2) : 10*(imax+2);
1210     Hs->SetMinimum(0.0); Hs->SetMaximum(ymx);
1211     if (irtype != 2) {
1212       TPad *pad1(pad);
1213       if (irtype != 1) {
1214     sprintf (name,"ENPad1%d%s%d",type,partsF[ipar].c_str(),iens[ien]);
1215     pad1 = new TPad(name,"pad1",0,0.3,1,1);
1216     pad1->SetBottomMargin(0.01);
1217       } 
1218       pad1->SetTopMargin(0.10); pad1->SetRightMargin(0.10);
1219       pad1->Draw(); pad1->cd();
1220       Hs->Draw("nostack");
1221       Hs->GetHistogram()->GetYaxis()->SetTitle("Events");
1222       if (irtype == 1) {
1223     Hs->GetHistogram()->GetXaxis()->SetTitle("Energy (GeV)");
1224     Hs->GetHistogram()->GetXaxis()->SetTitleOffset(0.90);
1225     Hs->GetHistogram()->GetYaxis()->SetTitleOffset(1.20);
1226       }
1227       pad1->Update(); pad1->Modified();
1228       if (stat) {
1229     for (unsigned int k=0; k<hists.size(); ++k) {
1230       TPaveStats* st1 = (TPaveStats*)hists[k]->GetListOfFunctions()->FindObject("stats");
1231       if (st1 != NULL) {
1232         st1->SetFillColor(kWhite);  
1233         st1->SetLineColor(icol[k]); st1->SetTextColor(icol[k]);
1234         st1->SetY1NDC(yh-dy);  st1->SetY2NDC(yh);
1235         st1->SetX1NDC(xh-0.3); st1->SetX2NDC(xh);
1236         yh -= dy;
1237       }
1238     }
1239       }
1240       legend->Draw("same");
1241       if (approve) {
1242     double yht = ((types[ien] == 2) || (!stat)) ? (0.88-legsize*(nmod+1)) : 
1243       ((types[ien] == 1) ? (yh1-0.02) : (yh1-legsize*(nmod+1)-0.02));
1244     TPaveText* text = new TPaveText(0.60, yht-0.04, 0.89, yht, "brNDC");
1245     text->AddText("CMS Preliminary");
1246     text->Draw("same");
1247 //  TPaveText* txt2 = new TPaveText(0.60, yht-0.08, 0.89,yht-0.04, "brNDC");
1248     TPaveText* txt2 = new TPaveText(0.11, 0.85, 0.35, 0.89, "brNDC");
1249     txt2->AddText("2006 Test Beam Data");
1250     txt2->Draw("same");
1251       }
1252     }
1253 
1254     pad->cd();
1255     if (irtype != 1) {
1256       sprintf (name,"ENStack2%d%s%d",type,partsF[ipar].c_str(),iens[ien]);
1257       THStack *Hsr = new THStack(name,"");
1258       Hsr->SetMinimum(0.0); Hsr->SetMaximum(3.99);
1259       TPad *pad2(pad);
1260       if (irtype != 2) {
1261     sprintf (name,"ENPad2%d%s%d",type,partsF[ipar].c_str(),iens[ien]);
1262     pad2  = new TPad(name,"pad2",0,0,1,0.3);
1263     pad2->SetTopMargin(0.005);
1264       }
1265       pad2->SetBottomMargin(0.20); pad2->SetRightMargin(0.10);
1266       pad2->Draw(); pad2->cd();
1267       double xh2 = (types[ien] == 2) ? 0.90 : xh;
1268       double yh2 = (irtype == 2) ? 0.88 : 0.98;
1269       double xh1 = (approve) ? (xh2-0.40) : ((irtype == 2) ? (xh2-0.60) : (xh2-0.30));
1270       double yh1 = (irtype == 2) ? (yh2-0.05*nmod) : (yh2-0.10*nmod);
1271       TLegend *legend1 = new TLegend(xh1, yh1, xh2, yh2);
1272       legend1->SetFillColor(kWhite); legend1->SetBorderSize(0); 
1273       legend1->SetMargin(0.2);
1274       TH1D *h_ref =  (TH1D*)hists[0];
1275       int imx1 = (ymax > 100) ? 10 : 1;
1276       int imx2 = (ymax > 100) ? 100 : 5;
1277       for (unsigned int i=1; i<hists.size(); i++) {
1278     sprintf (name,"ENRatio%d%s%d",i,partsF[ipar].c_str(),iens[ien]);
1279     TH1D *ratio = new TH1D(name, "Ratio", nbin, xmin, xmax);
1280     double sumNum(0), sumDen(0);
1281     for (int k=1; k<=nbin; ++k) {
1282       if (h_ref->GetBinContent(k) > imx1 && hists[i]->GetBinContent(k) > imx1) {
1283         double rat = hists[i]->GetBinContent(k)/h_ref->GetBinContent(k);
1284         double drt = rat*(h_ref->GetBinError(k)/h_ref->GetBinContent(k));
1285         ratio->SetBinContent(k,rat); ratio->SetBinError(k,drt);
1286         if (h_ref->GetBinContent(k) > imx2) {
1287           if (rat > 1.) {
1288         rat = 1./rat; drt *= (rat*rat);
1289           }
1290           sumNum += (fabs(1.0-rat)/(drt*drt));
1291           sumDen += (1.0/(drt*drt));
1292         }
1293       }
1294     }
1295     double mean  = (sumDen>0) ? (sumNum/sumDen) : 0;
1296     double error = (sumDen>0) ? 1.0/sqrt(sumDen) : 0;
1297     std::cout << "Model " << i << " Delta " << mean << " +- " << error 
1298           << std::endl;
1299     if (approve) {
1300       sprintf (name, "%s",modelNames[i-1].c_str());
1301     } else if (irtype == 2) {
1302       sprintf (name, "#delta = %6.3f #pm %6.3f  %s",mean,error,modelNames[i-1].c_str());
1303     } else {
1304       sprintf (name, "#delta = %6.3f #pm %6.3f",mean,error);
1305     }
1306     legend1->AddEntry(ratio,name,"lp");
1307     ratio->SetStats(0);
1308     ratio->SetLineColor(icol[i]);
1309     ratio->SetLineWidth(2);
1310     ratio->SetMarkerColor(icol[i]);
1311     ratio->SetMarkerStyle(mtype[i]);
1312     ratio->SetMarkerSize(1.2);
1313     ratio->SetLineWidth(2);
1314     Hsr->Add(ratio,"pe same");
1315       }
1316       Hsr->Draw("nostack");
1317       pad2->Update();
1318       Hsr->GetHistogram()->SetStats(0);
1319       Hsr->GetHistogram()->GetYaxis()->SetTitle("MC/Data");
1320       Hsr->GetHistogram()->GetXaxis()->SetTitle("Energy (GeV)");
1321       if (irtype == 2) {
1322     Hsr->GetHistogram()->GetXaxis()->SetTitleSize(0.06);
1323     Hsr->GetHistogram()->GetYaxis()->SetTitleSize(0.06);
1324     Hsr->GetHistogram()->GetXaxis()->SetTitleOffset(1.00);
1325     Hsr->GetHistogram()->GetYaxis()->SetTitleOffset(0.60);
1326     Hsr->GetHistogram()->GetXaxis()->SetLabelSize(0.07);
1327     Hsr->GetHistogram()->GetYaxis()->SetLabelSize(0.07);
1328       } else {
1329     Hsr->GetHistogram()->GetXaxis()->SetTitleSize(0.10);
1330     Hsr->GetHistogram()->GetYaxis()->SetTitleSize(0.10);
1331     Hsr->GetHistogram()->GetXaxis()->SetTitleOffset(0.90);
1332     Hsr->GetHistogram()->GetYaxis()->SetTitleOffset(0.35);
1333     Hsr->GetHistogram()->GetXaxis()->SetLabelSize(0.10);
1334     Hsr->GetHistogram()->GetYaxis()->SetLabelSize(0.10);
1335       }
1336       TLine *line = new TLine(xmin,1.0,xmax,1.0);
1337       line->SetLineStyle(2); line->SetLineWidth(2);
1338       line->SetLineColor(kBlack); line->Draw();
1339       if (approve) {
1340     double dyt = (irtype == 3) ? 0.08 : 0.08;
1341     TPaveText* text = new TPaveText(0.11,yh2-dyt,xh2-0.29,yh2,"brNDC");
1342     sprintf (name, "2006 Test Beam Data (%d GeV %s %s)", iens[ien], 
1343          partsN[ipar].c_str(), titlty[type-1].c_str());
1344     text->AddText(name);
1345     text->Draw("same");
1346     TPaveText* txt2 = new TPaveText(xh2-0.25, yh2-dyt, xh2, yh2, "brNDC");
1347     txt2->AddText("CMS Preliminary");
1348     txt2->Draw("same");
1349       } else {
1350     legend1->Draw("same");
1351       }
1352       pad2->Update();
1353     }
1354     pad->Update();
1355   }
1356   return pad;
1357 }
1358 
1359 void plotDataMCDist(int ipar, int ien, std::vector<int> models, int rebin=2,
1360             std::string prefix="model", bool approve=true, 
1361             bool stat=true, double xmin0=-1, double xmax0=-1, 
1362             int save=0) {
1363 
1364   for (int type=1; type<3; ++type) {
1365     char filename[100];
1366     TCanvas* c1 = plotDataMC(ipar, ien, models, type, rebin, 1, prefix, approve, stat, xmin0, xmax0, true);
1367     if (save != 0) {
1368       if (save < 0) sprintf (filename, "%s.jpg", c1->GetName());
1369       else          sprintf (filename, "%s.pdf", c1->GetName());
1370       c1->Print(filename);
1371     }
1372     TCanvas* c2 = plotDataMC(ipar, ien, models, type, rebin, 2, prefix, approve, stat, xmin0, xmax0, true);
1373     if (save != 0) {
1374       if (save < 0) sprintf (filename, "%s.jpg", c2->GetName());
1375       else          sprintf (filename, "%s.pdf", c2->GetName());
1376       c2->Print(filename);
1377     }
1378   }
1379 }
1380 
1381 void plotDataMCDist(std::string prefix, int type, int irtype,
1382             std::vector<int> models, int rebin=2, bool approve=true, 
1383             bool stat=true, double xmin0=-1, double xmax0=-1,
1384             int save=0) {
1385   int ienMin[6] = {0,0,2,3,0,0};
1386   int ienMax[6] = {8,8,6,6,8,8};
1387 
1388   char filename[100];
1389   for (int ipar=0; ipar<6; ++ipar) {
1390     for (int ien=ienMin[ipar]; ien<ienMax[ipar]; ++ien) {
1391       TCanvas* c1 = plotDataMC(ipar, ien, models, type, rebin, irtype, prefix, 
1392                    approve, stat, xmin0, xmax0, true);
1393       if (save != 0) {
1394     if (save < 0) sprintf (filename, "%s.jpg", c1->GetName());
1395     else          sprintf (filename, "%s.pdf", c1->GetName());
1396     c1->Print(filename);
1397       }
1398     }
1399   }
1400 }
1401 
1402 TCanvas* plotDataMC(int ipar, std::vector<int> models, bool ratio=false,
1403             std::string dirName="RespAll", bool approve=false) {
1404 
1405   static const int nmodels=37;
1406   std::string names[nmodels] = {"Test Beam Data",
1407                 "G4 10.0.p02 QGSP_FTFP_BERT_EML",
1408                 "G4 10.0.p02 FTFP_BERT_EML",
1409                 "G4 10.2.p02 QGSP_FTFP_BERT_EMM",
1410                 "G4 10.2.p02 FTFP_BERT_EMM",
1411                 "G4 10.2.p02 FTFP_BERT_ATL_EMM",
1412                 "G4 10.4.beta FTFP_BERT_EMM",
1413                 "G4 10.3.ref08 FTFP_BERT_EMM",
1414                 "G4 10.3.ref09 FTFP_BERT_EMM",
1415                 "G4 10.3.ref10 FTFP_BERT_EMM",
1416                 "G4 10.3.ref11 FTFP_BERT_EMM",
1417                 "G4 10.3.p03 FTFP_BERT_EMM",
1418                 "G4 10.4.cand01 FTFP_BERT_EMM",
1419                 "G4 10.4.cand02 FTFP_BERT_EMM",
1420                 "G4 10.4 FTFP_BERT_EMM",
1421                 "G4 10.2.p02 FTFP_BERT_EMM 10.0.pre3",
1422                 "G4 10.4 FTFP_BERT_EMM 10.0.pre3",
1423                 "G4 10.4 FTFP_BERT_EMM VecGeom 10.0.pre3",
1424                 "G4 10.4 FTFP_BERT_EMM VecGeom+CLHEP 10.0.pre3",
1425                 "G4 10.4 FTFP_BERT_EMM CLHEP",
1426                 "G4 10.4 FTFP_BERT_EMM VecGeom+CLHEP",
1427                 "G4 10.4.ref01 FTFP_BERT_EMM",
1428                 "G4 10.4.ref02 FTFP_BERT_EMM",
1429                 "G4 10.5.0.beta FTFP_BERT_EMM",
1430                 "G4 10.4.ref07 FTFP_BERT_EMM",
1431                 "G4 10.4.ref08 FTFP_BERT_EMM",
1432                 "G4 10.4.ref09 FTFP_BERT_EMM",
1433                 "G4 10.5.cand01 FTFP_BERT_EMM",
1434                 "G4 10.4.ref00 FTFP_BERT_EMM",
1435                 "G4 10.4.p03 FTFP_BERT_EMM",
1436                 "G4 10.5.ref00 FTFP_BERT_EMM",
1437                 "G4 10.5.ref01 FTFP_BERT_EMM",
1438                 "G4 10.5.ref02 FTFP_BERT_EMM",
1439                 "G4 10.5.ref03 FTFP_BERT_EMM",
1440                 "G4 10.5.ref04 FTFP_BERT_EMM",
1441                 "G4 10.5.ref05 FTFP_BERT_EMM",
1442                 "G4 10.5.ref06 FTFP_BERT_EMM"
1443   };
1444   std::string partsM[6] = {"pim","pip","km","kp","prop","prom"};
1445   std::string partsN[6] = {"#pi^{-}","#pi^{+}","K^{-}","K^{+}","proton","antiproton"};
1446   double      xmax[6]   = {400.0,30.0,15.0,15.0,400.0,15.0};
1447   double      ylow[6]   = {0.4,0.4,0.2,0.2,0.2,0.6};
1448   double      ylowr[6]  = {0.9,0.9,0.7,0.5,0.9,0.9};
1449   double      ymaxr[6]  = {1.2,1.2,1.2,1.4,1.2,1.2};
1450   int         colors[nmodels+1]= { 1, 2, 7, 6, 4, 9, 8,40, 7, 6, 9,46,48,30, 2,
1451                    6, 2, 7, 4, 2, 6, 4, 1, 7, 6, 7, 9, 8, 2, 9,
1452                    6, 7,30,30, 7,30, 7};
1453   int         mtype[nmodels+1] = {20,21,22,23,24,25,26,27,22,23,25,29,33,42,21,
1454                   21,22,23,24,21,22,23,24,21,22,23,33,26,22,21,
1455                   23,33,21,21,33,21,33};
1456 
1457   TCanvas* canvas(0);
1458   char     cname[100];
1459   int      nm(0);
1460   for (unsigned int i=0; i<models.size(); ++i) {
1461     if (models[i] < nmodels) ++nm;
1462   }
1463   double ymax = 0.948;
1464   double ymin = ymax-0.04*nm;
1465   TLegend*  legend(0);
1466   if (ratio) {
1467     legend = new TLegend(0.175, ymin, 0.970, ymax);
1468   } else {
1469     legend = new TLegend(0.175, ymin, 0.973, ymax);
1470   }
1471   legend->SetBorderSize(0); legend->SetFillColor(kWhite);
1472   legend->SetMargin(0.2);
1473   std::vector<TGraphAsymmErrors*> graphs;
1474   const int NPT=20;
1475   double    mom[NPT], dmom[NPT], momd[NPT], meand[NPT], dmeand[NPT], mean[NPT], dmean[NPT];
1476 
1477   // First data
1478   char                infile[100];
1479   if (dirName == "") {
1480     sprintf(infile,"%s.txt",  partsM[ipar].c_str());
1481   } else {
1482     sprintf(infile,"%s/%s.txt",  dirName.c_str(), partsM[ipar].c_str());
1483   }
1484   double    pb, resp, errsp;
1485   bool      ok(false);
1486   int       nptd(0);
1487   //First data
1488   std::ifstream fInput1(infile);
1489   if (!fInput1.good()) {
1490     std::cout << "Cannot open file " << infile << std::endl;
1491   } else {
1492     ok = true;
1493     while (1) {
1494       fInput1 >> pb >> resp >> errsp;
1495       if (!fInput1.good()) break;
1496       momd[nptd] = pb; dmom[nptd] = 0; meand[nptd] = resp; dmeand[nptd] = errsp; ++nptd;
1497     }
1498     fInput1.close();
1499     ok = (nptd > 0);
1500     if (debug_) {
1501       std::cout << "Reads " << nptd << " points from " << infile << std::endl;
1502       for (int k=0; k<nptd; ++k)
1503     std::cout << "[" << k << "] " << momd[k] << " +- " << dmom[k] << "   "
1504           << meand[k] << " +- " << dmeand[k] << std::endl;
1505     }
1506   }
1507   if (ok) {
1508     if (!ratio) {
1509       TGraphAsymmErrors *graph = new TGraphAsymmErrors(nptd,momd,meand,dmom,dmom,dmeand,dmeand);
1510       graph->SetMarkerStyle(mtype[0]);
1511       graph->SetMarkerColor(colors[0]);
1512       graph->SetMarkerSize(1.4);
1513       graph->SetLineColor(colors[0]);
1514       graph->SetLineWidth(2);
1515       graphs.push_back(graph);
1516       sprintf(cname, "2006 %s (%s)", names[0].c_str(), partsN[ipar].c_str());
1517       legend->AddEntry(graph, cname, "lp");
1518     }
1519     for (unsigned int k=0; k<models.size(); ++k) {
1520       int npt(0);
1521       double sumNum(0), sumDen(0);
1522       if (models[k] < nmodels) {
1523     int i = models[k];
1524     if (dirName == "") {
1525       sprintf(infile,"%sm%d.txt",partsM[ipar].c_str(),i);
1526     } else {
1527       sprintf(infile,"%s/%sm%d.txt",dirName.c_str(),partsM[ipar].c_str(),i);
1528     }
1529     std::ifstream fInput2(infile);
1530     if (!fInput2.good()) {
1531       std::cout << "Cannot open file " << infile << std::endl;
1532     } else {
1533       while (1) {
1534         fInput2 >> pb >> resp;
1535         if (!fInput2.good()) break;
1536         if (ratio) {
1537           for (int k=0; k<nptd; ++k) {
1538         if (std::abs(momd[k]-pb) < 0.1) {
1539           double rat =  resp/meand[k]; double drt = dmeand[k]/meand[k];
1540           mom[npt] = pb; dmom[npt] = 0; mean[npt] = rat; dmean[npt] = drt;
1541           if (rat > 1.) {
1542             rat = 1./rat; drt *= (rat*rat);
1543           }
1544           sumNum += (fabs(1.0-rat)/(drt*drt));
1545           sumDen += (1.0/(drt*drt));
1546           ++npt; break;
1547         }
1548           }
1549         } else {
1550           mom[npt] = pb; dmom[npt] = 0; mean[npt] = resp; dmean[npt] = 0.001; ++npt;
1551         }
1552       }
1553     }
1554     fInput2.close();
1555     if (npt > 0) {
1556       if (debug_) {
1557         std::cout << "Reads " << npt << " points from " << infile 
1558               << std::endl;
1559         for (int k=0; k<npt; ++k)
1560           std::cout << "[" << k << "] " << mom[k] << " +- " << dmom[k]
1561             << "   " << mean[k] << " +- " << dmean[k] << std::endl;
1562       }
1563       TGraphAsymmErrors *graph = new TGraphAsymmErrors(npt,mom,mean,dmom,dmom,dmean,dmean);
1564       graph->SetMarkerStyle(mtype[i+1]);
1565       graph->SetMarkerColor(colors[i+1]);
1566       graph->SetMarkerSize(1.2);
1567       graph->SetLineColor(colors[i+1]);
1568       graph->SetLineWidth(2);
1569       graphs.push_back(graph);
1570       if (ratio) {
1571         double rat = (sumDen>0) ? (sumNum/sumDen) : 0;
1572         double err = (sumDen>0) ? sqrt(1./sumDen) : 0;
1573         std::cout << "Particle " << partsM[ipar] << " Model "
1574               << names[i+1] << " ratio " << rat << " +- " 
1575               << err << std::endl;
1576         if (approve) sprintf (cname, "%s", names[i+1].c_str());
1577         else         sprintf (cname, "(#delta = %5.3f) %s", rat, names[i+1].c_str());
1578       } else {
1579         sprintf (cname, "%s", names[i+1].c_str());
1580       }
1581       legend->AddEntry(graph, cname, "lp");
1582     }
1583       }
1584     }
1585     // Now prepare plot
1586     gStyle->SetCanvasBorderMode(0); gStyle->SetCanvasColor(kWhite);
1587     gStyle->SetPadColor(kWhite);    gStyle->SetFillColor(kWhite);
1588     gStyle->SetOptTitle(kFALSE);    gStyle->SetPadBorderMode(0);
1589     gStyle->SetCanvasBorderMode(0); gStyle->SetOptStat(0);
1590     if (ratio) {
1591       sprintf(cname, "c_%sRespR", partsM[ipar].c_str());
1592       canvas = new TCanvas(cname, cname, 500, 300);
1593     } else {
1594       sprintf(cname, "c_%sResp", partsM[ipar].c_str());
1595       canvas = new TCanvas(cname, cname, 500, 500);
1596     }
1597     canvas->SetTopMargin(0.05);       canvas->SetLogx();
1598     canvas->SetLeftMargin(0.15);      canvas->SetRightMargin(0.025);
1599     canvas->SetBottomMargin(0.14);
1600     TH1F *vFrame = canvas->DrawFrame(1.0, 0.01, xmax[ipar], 2.0);
1601     vFrame->GetXaxis()->SetRangeUser(1.0,xmax[ipar]);
1602     if (ratio) {
1603       if (approve) vFrame->GetYaxis()->SetRangeUser(ylowr[ipar],ymaxr[ipar]);
1604       else         vFrame->GetYaxis()->SetRangeUser(ylowr[ipar],1.4);
1605     } else {
1606       vFrame->GetYaxis()->SetRangeUser(ylow[ipar],1.0);
1607     }
1608     if (ratio) {
1609       vFrame->GetXaxis()->SetLabelSize(0.06);
1610       vFrame->GetYaxis()->SetLabelSize(0.06);
1611       vFrame->GetXaxis()->SetTitleSize(0.06);
1612       vFrame->GetYaxis()->SetTitleSize(0.06);
1613       vFrame->GetXaxis()->SetTitleOffset(0.9);
1614       vFrame->GetYaxis()->SetTitleOffset(1.0);
1615     } else {
1616       vFrame->GetXaxis()->SetLabelSize(0.04);
1617       vFrame->GetYaxis()->SetLabelSize(0.04);
1618       vFrame->GetXaxis()->SetTitleSize(0.04);
1619       vFrame->GetYaxis()->SetTitleSize(0.04);
1620       vFrame->GetXaxis()->SetTitleOffset(1.2);
1621       vFrame->GetYaxis()->SetTitleOffset(1.6);
1622     }
1623     vFrame->GetXaxis()->SetLabelOffset(0.0);
1624     vFrame->GetXaxis()->SetTitle("p_{Beam} (GeV/c)");
1625     if (ratio) {
1626       vFrame->GetYaxis()->SetTitle("MC/Data (Response)");
1627     } else {
1628       if (approve) vFrame->GetYaxis()->SetTitle("Mean of E_{Measured}/p_{Beam}");
1629       else         vFrame->GetYaxis()->SetTitle("<E_{Measured}>/p_{Beam}");
1630     }
1631     for (unsigned int ii=0; ii<graphs.size(); ++ii) {
1632       if (ii == 0 && !ratio) graphs[ii]->Draw("P");
1633       else                   graphs[ii]->Draw("LP");
1634     }
1635     if ((!approve) || (!ratio)) legend->Draw();
1636     double yvl = (approve && ratio) ? ymax : ymin;
1637     if ((!approve) || ratio) {
1638       TPaveText* text = new TPaveText(0.16,yvl-0.08, 0.70,yvl-0.01,"brNDC");
1639       sprintf(cname, "2006 %s (%s)", names[0].c_str(), partsN[ipar].c_str());
1640       text->AddText(cname);
1641       text->Draw("same");
1642     }
1643     TPaveText* txt2 = new TPaveText(0.70, yvl-0.08, 0.948, yvl-0.01, "brNDC");
1644     sprintf (cname, "CMS Preliminary");
1645     txt2->AddText(cname);
1646     txt2->Draw("same");
1647     if (ratio) {
1648       TLine *line = new TLine(1.0,1.0,xmax[ipar],1.0);
1649       line->SetLineStyle(2); line->SetLineWidth(2);
1650       line->SetLineColor(kBlack); line->Draw();
1651     }
1652   }
1653   return canvas;
1654 
1655 }
1656 
1657 void plotDataMCResp(std::vector<int> models, std::string dirName="RespAll", 
1658             bool approve=true, int save=0) {
1659 
1660   for (int k=0; k<6; ++k) {
1661     char filename[100];
1662     TCanvas* c1 = plotDataMC(k, models, false, dirName, approve);
1663     if (save != 0) {
1664       if (save < 0) sprintf (filename, "%s.jpg", c1->GetName());
1665       else          sprintf (filename, "%s.pdf", c1->GetName());
1666       c1->Print(filename);
1667     }
1668     TCanvas* c2 = plotDataMC(k, models, true, dirName, approve);
1669     if (save != 0) {
1670       if (save < 0) sprintf (filename, "%s.jpg", c2->GetName());
1671       else          sprintf (filename, "%s.pdf", c2->GetName());
1672       c2->Print(filename);
1673     }
1674   }
1675 }
1676 
1677 void convertel(std::string indir="nmodel0", std::string outdir="model0", int md=0) {
1678 
1679   TB06Analysis p01((indir+"/el50e.root"),  (outdir+"/el50e.root"), 6,10, md, 1.0, 1.0);
1680   p01.Loop();
1681   TB06Analysis p02((indir+"/el50h.root"),  (outdir+"/el50h.root"), 6,10, md, 1.0, 1.0);
1682   p02.Loop();
1683 }
1684 
1685 void convert(std::string indir="nmodel0", std::string outdir="model0", int md=0) {
1686 
1687   TB06Analysis p01((indir+"/km4.root"),  (outdir+"/km4.root"),    2, 3, md, 1.0, 1.0);
1688   p01.Loop();
1689   TB06Analysis p02((indir+"/km5.root"),  (outdir+"/km5.root"),    2, 4, md, 1.0, 1.0);
1690   p02.Loop();
1691   TB06Analysis p03((indir+"/km6.root"),  (outdir+"/km6.root"),    2, 5, md, 1.0, 1.0);
1692   p03.Loop();
1693   TB06Analysis p04((indir+"/km7.root"),  (outdir+"/km7.root"),    2, 6, md, 1.0, 1.0);
1694   p04.Loop();
1695   TB06Analysis p05((indir+"/kp5.root"),  (outdir+"/kp5.root"),    3, 4, md, 1.0, 1.0);
1696   p05.Loop();
1697   TB06Analysis p06((indir+"/kp6.root"),  (outdir+"/kp6.root"),    3, 5, md, 1.0, 1.0);
1698   p06.Loop();
1699   TB06Analysis p07((indir+"/kp7.root"),  (outdir+"/kp7.root"),    3, 6, md, 1.0, 1.0);
1700   p07.Loop();
1701   TB06Analysis p08((indir+"/pim2.root"),  (outdir+"/pim2.root"),   0, 0, md, 1.0, 1.0);
1702   p08.Loop();
1703   TB06Analysis p09((indir+"/pim3.root"),  (outdir+"/pim3.root"),   0, 1, md, 1.0, 1.0);
1704   p09.Loop();
1705   TB06Analysis p10((indir+"/pim4.root"),  (outdir+"/pim4.root"),   0, 2, md, 1.0, 1.0);
1706   p10.Loop();
1707   TB06Analysis p11((indir+"/pim5.root"),  (outdir+"/pim5.root"),   0, 3, md, 1.0, 1.0);
1708   p11.Loop();
1709   TB06Analysis p12((indir+"/pim6.root"),  (outdir+"/pim6.root"),   0, 4, md, 1.0, 1.0);
1710   p12.Loop();
1711   TB06Analysis p13((indir+"/pim7.root"),  (outdir+"/pim7.root"),   0, 5, md, 1.0, 1.0);
1712   p13.Loop();
1713   TB06Analysis p14((indir+"/pim8.root"),  (outdir+"/pim8.root"),   0, 6, md, 1.0, 1.0);
1714   p14.Loop();
1715   TB06Analysis p15((indir+"/pim9.root"),  (outdir+"/pim9.root"),   0, 7, md, 1.0, 1.0);
1716   p15.Loop();
1717   TB06Analysis p16((indir+"/pim20.root"),  (outdir+"/pim20.root"),  0, 8, md, 1.0, 1.0);
1718   p16.Loop();
1719   TB06Analysis p17((indir+"/pim30.root"),  (outdir+"/pim30.root"),  0, 9, md, 1.0, 1.0);
1720   p17.Loop();
1721   TB06Analysis p18((indir+"/pim50.root"),  (outdir+"/pim50.root"),  0,10, md, 1.0, 1.0);
1722   p18.Loop();
1723   TB06Analysis p19((indir+"/pim100.root"), (outdir+"/pim100.root"), 0,11, md, 1.0, 1.0);
1724   p19.Loop();
1725   TB06Analysis p20((indir+"/pim150.root"), (outdir+"/pim150.root"), 0,12, md, 1.0, 1.0);
1726   p20.Loop();
1727   TB06Analysis p21((indir+"/pim200.root"), (outdir+"/pim200.root"), 0,13, md, 1.0, 1.0);
1728   p21.Loop();
1729   TB06Analysis p22((indir+"/pim300.root"), (outdir+"/pim300.root"), 0,14, md, 1.0, 1.0);
1730   p22.Loop();
1731   TB06Analysis p23((indir+"/pip2.root"),  (outdir+"/pip2.root"),   1, 0, md, 1.0, 1.0);
1732   p23.Loop();
1733   TB06Analysis p24((indir+"/pip3.root"),  (outdir+"/pip3.root"),   1, 1, md, 1.0, 1.0);
1734   p24.Loop();
1735   TB06Analysis p25((indir+"/pip4.root"),  (outdir+"/pip4.root"),   1, 2, md, 1.0, 1.0);
1736   p25.Loop();
1737   TB06Analysis p26((indir+"/pip5.root"),  (outdir+"/pip5.root"),   1, 3, md, 1.0, 1.0);
1738   p26.Loop();
1739   TB06Analysis p27((indir+"/pip6.root"),  (outdir+"/pip6.root"),   1, 4, md, 1.0, 1.0);
1740   p27.Loop();
1741   TB06Analysis p28((indir+"/pip7.root"),  (outdir+"/pip7.root"),   1, 5, md, 1.0, 1.0);
1742   p28.Loop();
1743   TB06Analysis p29((indir+"/pip8.root"),  (outdir+"/pip8.root"),   1, 6, md, 1.0, 1.0);
1744   p29.Loop();
1745   TB06Analysis p30((indir+"/pip9.root"),  (outdir+"/pip9.root"),   1, 7, md, 1.0, 1.0);
1746   p30.Loop();
1747   TB06Analysis p31((indir+"/prom2.root"),  (outdir+"/prom2.root"),  5, 0, md, 1.0, 1.0);
1748   p31.Loop();
1749   TB06Analysis p32((indir+"/prom3.root"),  (outdir+"/prom3.root"),  5, 1, md, 1.0, 1.0);
1750   p32.Loop();
1751   TB06Analysis p33((indir+"/prom4.root"),  (outdir+"/prom4.root"),  5, 2, md, 1.0, 1.0);
1752   p33.Loop();
1753   TB06Analysis p34((indir+"/prom5.root"),  (outdir+"/prom5.root"),  5, 3, md, 1.0, 1.0);
1754   p34.Loop();
1755   TB06Analysis p35((indir+"/prom6.root"),  (outdir+"/prom6.root"),  5, 4, md, 1.0, 1.0);
1756   p35.Loop();
1757   TB06Analysis p36((indir+"/prom7.root"),  (outdir+"/prom7.root"),  5, 5, md, 1.0, 1.0);
1758   p36.Loop();
1759   TB06Analysis p37((indir+"/prom8.root"),  (outdir+"/prom8.root"),  5, 6, md, 1.0, 1.0);
1760   p37.Loop();
1761   TB06Analysis p38((indir+"/prom9.root"),  (outdir+"/prom9.root"),  5, 7, md, 1.0, 1.0);
1762   p38.Loop();
1763   TB06Analysis p39((indir+"/prop2.root"),  (outdir+"/prop2.root"),  4, 0, md, 1.0, 1.0);
1764   p39.Loop();
1765   TB06Analysis p40((indir+"/prop3.root"),  (outdir+"/prop3.root"),  4, 1, md, 1.0, 1.0);
1766   p40.Loop();
1767   TB06Analysis p41((indir+"/prop4.root"),  (outdir+"/prop4.root"),  4, 2, md, 1.0, 1.0);
1768   p41.Loop();
1769   TB06Analysis p42((indir+"/prop5.root"),  (outdir+"/prop5.root"),  4, 3, md, 1.0, 1.0);
1770   p42.Loop();
1771   TB06Analysis p43((indir+"/prop6.root"),  (outdir+"/prop6.root"),  4, 4, md, 1.0, 1.0);
1772   p43.Loop();
1773   TB06Analysis p44((indir+"/prop7.root"),  (outdir+"/prop7.root"),  4, 5, md, 1.0, 1.0);
1774   p44.Loop();
1775   TB06Analysis p45((indir+"/prop8.root"),  (outdir+"/prop8.root"),  4, 6, md, 1.0, 1.0);
1776   p45.Loop();
1777   TB06Analysis p46((indir+"/prop9.root"),  (outdir+"/prop9.root"),  4, 7, md, 1.0, 1.0);
1778   p46.Loop();
1779   TB06Analysis p47((indir+"/prop20.root"), (outdir+"/prop20.root"), 4, 8, md, 1.0, 1.0);
1780   p47.Loop();
1781   TB06Analysis p48((indir+"/prop30.root"), (outdir+"/prop30.root"), 4, 9, md, 1.0, 1.0);
1782   p48.Loop();
1783   TB06Analysis p49((indir+"/prop350.root"),(outdir+"/prop350.root"),4,15, md, 1.0, 1.0);
1784   p49.Loop();
1785 }
1786 
1787 void makePlots(int model1, int model2, int model3=-1, int model4=-1,
1788            int model5=-1, int save=0) {
1789   std::vector<int> models;
1790   if (model1 >= 0) models.push_back(model1);
1791   if (model2 >= 0) models.push_back(model2);
1792   if (model3 >= 0) models.push_back(model3);
1793   if (model4 >= 0) models.push_back(model4);
1794   if (model5 >= 0) models.push_back(model5);
1795   for (int ien=0; ien<8; ++ien)
1796    plotDataMCDist(0,ien,models,2,"model",true,true,-1,-1,save);
1797    for (int ien=0; ien<8; ++ien)
1798     plotDataMCDist(1,ien,models,2,"model",true,true,-1,-1,save);
1799   /*
1800   for (int ien=2; ien<6; ++ien)
1801     plotDataMCDist(2,ien,models,2,"model",true,true,-1,-1,save);
1802   for (int ien=3; ien<6; ++ien)
1803     plotDataMCDist(3,ien,models,2,"model",true,true,-1,-1,save);
1804   */
1805   for (int ien=0; ien<8; ++ien)
1806     plotDataMCDist(4,ien,models,2,"model",true,true,-1,-1,save);
1807   /*
1808   for (int ien=0; ien<8; ++ien)
1809     plotDataMCDist(5,ien,models,2,"model",true,true,-1,-1,save);
1810   */
1811 }
1812 
1813 void PrintMeanRatio(const char* infile="ResMean.out", int maxset=-1) {
1814   std::ifstream fInput(infile);
1815   if (!fInput.good()) {
1816     std::cout << "Cannot open file " << infile << std::endl;
1817   } else {
1818     int nset(0);
1819     fInput >> nset;
1820     if (fInput.good()) {
1821       int nsets = ((maxset < 0) || (maxset >= nset)) ? nset-1 : maxset;
1822       std::cout << "Will look into " << nsets << ":" << nset << " sets\n";
1823       for (int iset=0; iset<=nsets; ++iset) {
1824     char buff[200];
1825     int  npt(0);
1826     fInput >> npt >> buff;
1827     std::cout << "Data set " << iset << " " << buff << " with " << npt 
1828           << " points" << std::endl;
1829     double sumNum(0), sumDen(0);
1830     for (int ip=0; ip<npt; ++ip) {
1831       double pbeam, rat, drt;
1832       fInput >> pbeam >> rat >> drt;
1833       if (rat > 1.) {
1834         rat = 1./rat; drt *= (rat*rat);
1835       }
1836       sumNum += (fabs(1.0-rat)/(drt*drt));
1837       sumDen += (1.0/(drt*drt));
1838     }
1839     double mean  = (sumDen>0) ? (sumNum/sumDen) : 0;
1840     double error = (sumDen>0) ? 1.0/sqrt(sumDen) : 0;
1841     std::cout << "Set " << iset << " Delta " << mean << " +- " << error 
1842           << std::endl;
1843       }
1844     }
1845     fInput.close();
1846   }
1847 }