HcalCalibAlgos/macros/RecJet.C

0001 //////////////////////////////////////////////////////////
0002 // This class has been automatically generated on
0003 // Thu Jan 29 14:18:27 2015 by ROOT version 5.34/19
0004 // from TTree RecJet/RecJet Tree
0005 // found on file: True.root
0006 //////////////////////////////////////////////////////////
0007
0008 #include <TCanvas.h>
0009 #include <TChain.h>
0010 #include <TFile.h>
0011 #include <TGraph.h>
0012 #include <TH1D.h>
0013 #include <TH2D.h>
0014 #include <TPaveStats.h>
0015 #include <TROOT.h>
0016 #include <TStyle.h>
0017 #include <TTree.h>
0018 #include <iostream>
0019 #include <string>
0020 #include <cmath>
0021 #include <map>
0022 #include <fstream>
0023 #include <TF1.h>
0024
0025 #include <iomanip>
0026 #include <sstream>
0027
0028 // Header file for the classes stored in the TTree if any.
0029
0030 // Fixed size dimensions of array or collections stored in the TTree if any.
0031
0032 class RecJet {
0033 public :
0034   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
0035   Int_t           fCurrent; //!current Tree number in a TChain
0036
0037   // Declaration of leaf types
0038   Int_t           cells;
0039   Int_t           mysubd;
0040   Int_t           depth;
0041   Int_t           ieta;
0042   Int_t           iphi;
0043   Float_t         mom0_MB;
0044   Float_t         mom1_MB;
0045   Float_t         mom2_MB;
0046   Float_t         mom3_MB;
0047   Float_t         mom4_MB;
0048   Float_t         mom0_Diff;
0049   Float_t         mom1_Diff;
0050   Float_t         mom2_Diff;
0051   Float_t         mom3_Diff;
0052   Float_t         mom4_Diff;
0053   Int_t           trigbit;
0054   Double_t        rnnumber;
0055
0056   // List of branches
0057   TBranch        *b_cells;     //!
0058   TBranch        *b_mysubd;    //!
0059   TBranch        *b_depth;     //!
0060   TBranch        *b_ieta;      //!
0061   TBranch        *b_iphi;      //!
0062   TBranch        *b_mom0_MB;   //!
0063   TBranch        *b_mom1_MB;   //!
0064   TBranch        *b_mom2_MB;   //!
0065   TBranch        *b_mom3_MB;   //!
0066   TBranch        *b_mom4_MB;   //!
0067   TBranch        *b_mom0_Diff;   //!
0068   TBranch        *b_mom1_Diff;   //!
0069   TBranch        *b_mom2_Diff;   //!
0070   TBranch        *b_mom3_Diff;   //!
0071   TBranch        *b_mom4_Diff;
0072   TBranch        *b_trigbit;   //!
0073   TBranch        *b_rnnumber;  //!
0074
0075   struct MyInfo {
0076     double Mom0, Mom1, Mom2, Mom3, Mom4;
0077     MyInfo() {
0078       Mom0 = Mom1 = Mom2 = Mom3 = Mom4 = 0.;
0079     }
0080   };
0081
0082   struct CFactors {
0083     double cfac1, efac1, cfac2, efac2;
0084     CFactors() {
0085       cfac1 = cfac2 = 1;
0086       efac1 = efac2 = 0;
0087     }
0088   };
0089
0090   struct Hists {
0091     TH1D *h1, *h2, *h3, *h4, *h5, *h6, *h7, *h8, *h9;
0092     Hists() {
0093       h1 = h2 = h3 = h4 = h5= h6= h7= h8= h9=0;
0094     }
0095   };
0096   TFile*                                           file;
0097   std::string                                      detector;
0098   int                                              mode_;
0099   bool                                             loadTrig_;
0100   std::map<unsigned int,MyInfo>                    mNoise_;
0101   std::map<std::pair<unsigned int,int>,MyInfo>     mTrig_;
0102   std::vector<unsigned int>                        dets_;
0103   double                                           factor_;
0104   double                                           err_mean, err_var;
0105   std::map<unsigned int,CFactors>                  corrFactor_;
0106
0107   RecJet(std::string fname, int mode=0);
0108   virtual ~RecJet();
0109   bool             OpenFile(std::string fname);
0110   virtual Int_t    GetEntry(Long64_t entry);
0111   virtual Long64_t LoadTree(Long64_t entry);
0112   virtual void     Init(TTree *tree);
0113   virtual Bool_t   Notify();
0114   virtual void     Show(Long64_t entry = -1);
0115   virtual Int_t    Cut(Long64_t entry);
0116   std::map<unsigned int,RecJet::MyInfo> LoopMap();
0117   virtual void     LoopNoise();
0118   virtual void     Loop(int subdet, std::string indx, bool clear);
0119   virtual void     LoopIter(int subdet, std::string indx, bool clear, int maxIter=100);
0120   virtual void     LoopIterate(int subdet, std::string indx, double emin, double emax, int maxIter=100);
0121     std::map<unsigned int,RecJet::Hists> MakeCorr(std::map<unsigned int,RecJet::MyInfo>&, int, int);
0122   virtual void     LoopMean(int subdet, std::string indx);
0123   virtual void     eta_distribution(std::string fname="Eta.root", bool var=true);
0124   virtual void     Alleta_distribution(std::string fname="AllEta.root", bool var=true, bool HBHE=true);
0125   virtual void     det_distribution(std::string fname="Det.root");
0126   virtual void     Fit(std::string rootfile="HCALNZS2015_Final_Pedestal_magnet_test.root", std::string textfile="fit.txt");
0127   std::pair<double,double> SubtractNoise(unsigned int, double, double, bool);
0128   virtual void     Disturb(int subd);
0129   std::map<unsigned int,RecJet::MyInfo> MultiplyEnergy(std::map<unsigned int,RecJet::MyInfo>&,int);
0130   virtual void     ChangeMoments(std::map<unsigned int,RecJet::MyInfo>::iterator&, double);
0131   std::map<unsigned int,RecJet::MyInfo> LoadNoise();
0132   virtual void     MeanVariance(std::map<unsigned int,RecJet::MyInfo>::iterator &mitr, std::pair<double,double>& mean, std::pair<double,double>& variance);
0133   std::vector<std::pair<double,double> > Staggered_CorrFactor(int subd, std::map<unsigned int,RecJet::Hists>::iterator &hitr, bool varmethod);
0134   std::vector<std::pair<double,double> > CorrFactor(int subd, std::map<unsigned int,RecJet::Hists>::iterator &hitr, bool varmethod);
0135   void StoreCorrFactor(const std::map<unsigned int,RecJet::Hists> &, std::map<unsigned int,RecJet::CFactors>&, bool clear=true);
0136   std::pair<double,double> MeanDeviation(int, std::map<unsigned int,RecJet::CFactors>&);
0137   void             WriteCorrFactor(std::string& outFile, std::string& infile, bool var=false);
0138   virtual void     Error_study(std::string fname, std::string det, int depth);
0139   virtual void     LoopMapTrig();
0140   virtual void     LoopTrig(int sdet, int eta, int phi, int dep, std::string indx);
0141   std::map<unsigned int,RecJet::MyInfo> BuildMap(int sdet, bool cfirst, double emin, double emax);
0142 };
0143
0144 RecJet::RecJet(std::string fname, int mode) : fChain(0), file(0), mode_(mode),
0145                           loadTrig_(false), factor_(1.) {
0146   // if parameter tree is not specified (or zero), connect the file
0147   // used to generate this class and read the Tree.
0148   std::cout << "Open file " << fname << std::endl;
0149   //  std::cout << "hi.." << std::endl;
0150   file = new TFile(fname.c_str());
0151   //file->cd("minbiasana");
0152   //  std::cout << "hi.." << std::endl;
0153   TTree *tree;
0154   //gDirectory->GetObject("RecJet", tree);
0155   file->GetObject("RecJet", tree);
0156   Init(tree);
0157 }
0158
0159 RecJet::~RecJet() {
0160   if (!fChain) return;
0161   delete fChain->GetCurrentFile();
0162 }
0163
0164 bool RecJet::OpenFile(std::string fname) {
0165   if (fChain) {
0166     std::cout << "Close File " << fChain->GetCurrentFile()->GetName() << std::endl;
0167     delete fChain->GetCurrentFile();
0168   }
0169   file = new TFile(fname.c_str());
0170   if (file) {
0171     file->cd("minbiasana");
0172     TTree *tree;
0173     gDirectory->GetObject("RecJet", tree);
0174     Init(tree);
0175     return true;
0176   } else {
0177     return false;
0178   }
0179 }
0180
0181 Int_t RecJet::GetEntry(Long64_t entry) {
0182   // Read contents of entry.
0183   if (!fChain) return 0;
0184   return fChain->GetEntry(entry);
0185 }
0186
0187 Long64_t RecJet::LoadTree(Long64_t entry) {
0188   // Set the environment to read one entry
0189   if (!fChain) return -5;
0190   Long64_t centry = fChain->LoadTree(entry);
0191   if (centry < 0) return centry;
0192   if (fChain->GetTreeNumber() != fCurrent) {
0193     fCurrent = fChain->GetTreeNumber();
0194     Notify();
0195   }
0196   return centry;
0197 }
0198
0199 void RecJet::Init(TTree *tree) {
0200   // The Init() function is called when the selector needs to initialize
0201   // a new tree or chain. Typically here the branch addresses and branch
0202   // pointers of the tree will be set.
0203   // It is normally not necessary to make changes to the generated
0204   // code, but the routine can be extended by the user if needed.
0205   // Init() will be called many times when running on PROOF
0206   // (once per file to be processed).
0207
0208   // Set branch addresses and branch pointers
0209   if (!tree) return;
0210   fChain = tree;
0211   fCurrent = -1;
0212   fChain->SetMakeClass(1);
0213
0214   fChain->SetBranchAddress("cells", &cells, &b_cells);
0215   fChain->SetBranchAddress("mysubd", &mysubd, &b_mysubd);
0216   fChain->SetBranchAddress("depth", &depth, &b_depth);
0217   fChain->SetBranchAddress("ieta", &ieta, &b_ieta);
0218   fChain->SetBranchAddress("iphi", &iphi, &b_iphi);
0219   fChain->SetBranchAddress("mom0_MB", &mom0_MB, &b_mom0_MB);
0220   fChain->SetBranchAddress("mom1_MB", &mom1_MB, &b_mom1_MB);
0221   fChain->SetBranchAddress("mom2_MB", &mom2_MB, &b_mom2_MB);
0222   fChain->SetBranchAddress("mom3_MB", &mom3_MB, &b_mom3_MB);
0223   fChain->SetBranchAddress("mom4_MB", &mom4_MB, &b_mom4_MB);
0224   if (mode_ == 1) {
0225     fChain->SetBranchAddress("mom0_Diff", &mom0_Diff, &b_mom0_Diff);
0226     fChain->SetBranchAddress("mom1_Diff", &mom1_Diff, &b_mom1_Diff);
0227     fChain->SetBranchAddress("mom2_Diff", &mom2_Diff, &b_mom2_Diff);
0228     fChain->SetBranchAddress("mom3_Diff", &mom3_Diff, &b_mom3_Diff);
0229     fChain->SetBranchAddress("mom4_Diff", &mom4_Diff, &b_mom4_Diff);
0230   }
0231   fChain->SetBranchAddress("trigbit", &trigbit, &b_trigbit);
0232   fChain->SetBranchAddress("rnnumber", &rnnumber, &b_rnnumber);
0233   Notify();
0234 }
0235
0236 Bool_t RecJet::Notify() {
0237   // The Notify() function is called when a new file is opened. This
0238   // can be either for a new TTree in a TChain or when when a new TTree
0239   // is started when using PROOF. It is normally not necessary to make changes
0240   // to the generated code, but the routine can be extended by the
0241   // user if needed. The return value is currently not used.
0242
0243   return kTRUE;
0244 }
0245
0246 void RecJet::Show(Long64_t entry) {
0247   // Print contents of entry.
0248   // If entry is not specified, print current entry
0249   if (!fChain) return;
0250   fChain->Show(entry);
0251 }
0252
0253 Int_t RecJet::Cut(Long64_t ) {
0254   // This function may be called from Loop.
0255   // returns  1 if entry is accepted.
0256   // returns -1 otherwise.
0257   return 1;
0258 }
0259
0260 std::map<unsigned int,RecJet::MyInfo> RecJet::LoopMap() {
0261   //   In a ROOT session, you can do:
0262   //      Root > .L RecJet.C
0263   //      Root > RecJet t
0264   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
0265   //      Root > t.Show();       // Show values of entry 12
0266   //      Root > t.Show(16);     // Read and show values of entry 16
0267   //      Root > t.Loop();       // Loop on all entries
0268   //
0269   //     This is the loop skeleton where:
0270   //    jentry is the global entry number in the chain
0271   //    ientry is the entry number in the current Tree
0272   //  Note that the argument to GetEntry must be:
0273   //    jentry for TChain::GetEntry
0274   //    ientry for TTree::GetEntry and TBranch::GetEntry
0275   //
0276   //       To read only selected branches, Insert statements like:
0277   // METHOD1:
0278   //    fChain->SetBranchStatus("*",0);  // disable all branches
0279   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
0280   // METHOD2: replace line
0281   //    fChain->GetEntry(jentry);       //read all branches
0282   //by  b_branchname->GetEntry(ientry); //read only this branch
0283
0284   std::map<unsigned int,RecJet::MyInfo> m_;
0285   if (fChain != 0) {
0286     Long64_t nentries = fChain->GetEntriesFast();
0287     Long64_t nbytes = 0, nb = 0;
0288
0289     for (Long64_t jentry=0; jentry<nentries;jentry++) {
0290       Long64_t ientry = LoadTree(jentry);
0291       if (ientry < 0) break;
0292
0293       nb = fChain->GetEntry(jentry);   nbytes += nb;
0294       unsigned int detId1 = ((mysubd<<20) | ((depth&0x1f)<<14) |
0295                  ((ieta>0)?(0x2000|(ieta<<7)):((-ieta)<<7)) |
0296                  (iphi&0x7f));
0297       std::map<unsigned int, RecJet::MyInfo>::iterator mitr = m_.find(detId1);
0298       if (mitr == m_.end()) {
0299     RecJet::MyInfo info;
0300     m_[detId1] = info ;
0301     mitr = m_.find(detId1) ;
0302       }
0303       if (mode_ == 1) {
0304     mitr->second.Mom0 += mom0_Diff;
0305     mitr->second.Mom1 += mom1_Diff;
0306     mitr->second.Mom2 += mom2_Diff;
0307     mitr->second.Mom3 += mom3_Diff;
0308     mitr->second.Mom4 += mom4_Diff;
0309       } else {
0310     mitr->second.Mom0 += mom0_MB;
0311     mitr->second.Mom1 += mom1_MB;
0312     mitr->second.Mom2 += mom2_MB;
0313     mitr->second.Mom3 += mom3_MB;
0314     mitr->second.Mom4 += mom4_MB;
0315       }
0316     }
0317   }
0318   return m_;
0319 }
0320
0321 void RecJet::LoopNoise() {
0322
0323   mNoise_ = LoopMap();
0324   if (mNoise_.size() == 0) return;
0325   TH1D* h[4];
0326   std::string det[4] = {"HB", "HE", "HO", "HF"};
0327   char fname[80], name[80], title[80];
0328   for (int subd=0; subd<4; ++subd) {
0329     sprintf (fname, "Pedestal_2015D_%s.root", det[subd].c_str());
0330     TFile f(fname, "recreate");
0331     detector = det[subd];
0332     ofstream myfile;
0333     sprintf (fname,"%s_Noise.txt",detector.c_str());
0334     myfile.open(fname);
0335     sprintf(name,"%s",det[subd].c_str());
0336     sprintf(title,"Energy Distribution for %s",det[subd].c_str());
0337     if (subd !=3) h[subd] = new TH1D(name,title,100,-1.,1.);
0338     else          h[subd] = new TH1D(name,title,100,0.,5.);
0339     std::map<unsigned int,RecJet::Hists> hh_ ;
0340     for (std::map<unsigned int,RecJet::MyInfo>::iterator mitr = mNoise_.begin();
0341      mitr !=mNoise_.end(); mitr++) {
0342       unsigned int sdet    = ((mitr->first)>>20)&0x7;
0343       if (sdet == (unsigned int)(subd+1)) {
0344     unsigned int detId2  = (mitr->first)&0x7fff80;
0345     std::map<unsigned int,RecJet::Hists>::iterator hitr = hh_.find(detId2);
0346     int keta = (((mitr->first)&0x2000) ? (((mitr->first)>>7)&0x3f) :
0347             -(((mitr->first)>>7)&0x3f));
0348     int dept = (((mitr->first)>>14)&0x1f);
0349     int kphi = ((mitr->first)&0x7f);
0350     if (hitr == hh_.end()) {
0351       RecJet::Hists hh;
0352       hh_[detId2]  = hh;
0353       hitr = hh_.find(detId2);
0354       sprintf (name,  "hist%s%d%d", detector.c_str(), keta, dept);
0355       sprintf (title, "#phi (%s #eta %d depth %d)", detector.c_str(), keta, dept);
0356       hitr->second.h1 = new TH1D(name, title, 72, 0, 72.);
0357       hitr->second.h1->GetXaxis()->SetTitle("i#phi");
0358       hitr->second.h1->GetXaxis()->CenterTitle();
0359       hitr->second.h1->GetYaxis()->SetTitle("Mean Energy (GeV)");
0360       hitr->second.h1->GetYaxis()->CenterTitle();
0361       sprintf (name,  "histv%s%d%d", detector.c_str(), keta, dept);
0362       hitr->second.h2 = new TH1D(name, title, 72, 0, 72.);
0363       hitr->second.h2->GetXaxis()->SetTitle("i#phi");
0364       hitr->second.h2->GetXaxis()->CenterTitle();
0365       hitr->second.h2->GetYaxis()->SetTitle("Variance (GeV^2)");
0366       hitr->second.h2->GetYaxis()->CenterTitle();
0367     }
0368     std::pair<double,double> mean, variance;
0369     MeanVariance(mitr,mean,variance);
0370         h[subd]->Fill(mean.first);
0371     hitr->second.h1->SetBinContent(kphi,mean.first);
0372     hitr->second.h1->SetBinError(kphi,mean.second);
0373     hitr->second.h2->SetBinContent(kphi,variance.first);
0374     hitr->second.h2->SetBinError(kphi,variance.second);
0375     myfile << mitr->first << "\t" << mitr->second.Mom0 << "\t"
0376            << mitr->second.Mom1 << "\t" << mitr->second.Mom2 << "\t"
0377            << mitr->second.Mom3 << "\t" << mitr->second.Mom4 << std::endl;
0378       }
0379     }
0380     h[subd]->Write();
0381     myfile.close();
0382     for (std::map<unsigned int,RecJet::Hists>::iterator hitr = hh_.begin();
0383      hitr != hh_.end(); ++hitr) {
0384       int keta = (((hitr->first)&0x2000) ? (((hitr->first)>>7)&0x3f) :
0385             -(((hitr->first)>>7)&0x3f));
0386       int dept = (((hitr->first)>>14)&0x1f);
0387       sprintf (fname, "correction_mean%s%d%d", detector.c_str(), keta, dept);
0388       sprintf (title, "Correction Factor (Mean) for (%s #eta %d depth %d)", detector.c_str(), keta, dept);
0389       hitr->second.h3 = new TH1D(fname, title, 72, 0., 72.);
0390       hitr->second.h3->GetXaxis()->SetTitle("i#phi");
0391       hitr->second.h3->GetXaxis()->CenterTitle();
0392       hitr->second.h3->GetYaxis()->SetTitle("Correction Factor");
0393       hitr->second.h3->GetYaxis()->CenterTitle();
0394       sprintf (fname, "correction_variance%s%d%d", detector.c_str(), keta,dept);
0395       sprintf (title, "Correction Factor (Variance) for (%s #eta %d depth %d)", detector.c_str(), keta, dept);
0396       hitr->second.h4 = new TH1D(fname, title, 72, 0., 72.);
0397       hitr->second.h4->GetXaxis()->SetTitle("i#phi");
0398       hitr->second.h4->GetXaxis()->CenterTitle();
0399       hitr->second.h4->GetYaxis()->SetTitle("Correction Factor");
0400       hitr->second.h4->GetYaxis()->CenterTitle();
0401       sprintf (fname, "correction_variance_mine%s%d%d", detector.c_str(), keta,dept);
0402       sprintf (title, "Correction Factor (Variance) for (%s #eta %d depth %d)", detector.c_str(), keta, dept);
0403       hitr->second.h8 = new TH1D(fname, title, 72, 0., 72.);
0404       hitr->second.h8->GetXaxis()->SetTitle("i#phi");
0405       hitr->second.h8->GetXaxis()->CenterTitle();
0406       hitr->second.h8->GetYaxis()->SetTitle("Corrected Variance");
0407       hitr->second.h8->GetYaxis()->CenterTitle();
0408       sprintf (fname, "correction_mean_mine%s%d%d", detector.c_str(), keta,dept);
0409       sprintf (title, "Correction Factor (Mean) for (%s #eta %d depth %d)", detector.c_str(), keta, dept);
0410       hitr->second.h9 = new TH1D(fname, title, 72, 0., 72.);
0411       hitr->second.h9->GetXaxis()->SetTitle("i#phi");
0412       hitr->second.h9->GetXaxis()->CenterTitle();
0413       hitr->second.h9->GetYaxis()->SetTitle("Corrected Mean");
0414       hitr->second.h9->GetYaxis()->CenterTitle();
0415       std::vector<std::pair<double,double> > corr = CorrFactor(subd, hitr, false);
0416       for (unsigned int i=0; i<corr.size(); ++i) {
0417     if (corr[i].first > 0){
0418       hitr->second.h3->SetBinContent(i+1, corr[i].first);
0419       hitr->second.h3->SetBinError(i+1,   corr[i].second);
0420     }
0421       }
0422       corr = CorrFactor(subd, hitr, true);
0423       for (unsigned int i=0; i<corr.size(); ++i) {
0424     if (corr[i].first > 0){
0425       hitr->second.h4->SetBinContent(i+1, corr[i].first);
0426       hitr->second.h4->SetBinError(i+1,   corr[i].second);
0427     }
0428       }
0429       hitr->second.h1->Write();
0430       hitr->second.h2->Write();
0431       hitr->second.h3->Write();
0432       hitr->second.h4->Write();
0433       hitr->second.h8->Write();
0434       hitr->second.h9->Write();
0435
0436     }
0437     f.Close();
0438   }
0439 }
0440
0441 void RecJet::Loop(int sdet, std::string indx, bool clear) {
0442
0443   std::map<unsigned int,RecJet::MyInfo> m_ = LoopMap();
0444
0445   if (m_.size() == 0) return;
0446
0447   char fname[80], name[80], title[80];
0448   m_ = MultiplyEnergy(m_,0);
0449
0450   int subd = ((sdet>=1 && sdet<=4) ? (sdet-1) : 0);
0451   std::string det[4] = {"HB", "HE", "HO", "HF"};
0452   sprintf (fname, "HCALNZS2015%s_phiRecoDistribution_Final%s.root", indx.c_str(), det[subd].c_str());
0453   TFile f(fname, "recreate");
0454   detector = det[subd];
0455   sprintf(name,"%s", det[subd].c_str());
0456   sprintf(title,"Energy Distribution for %s", det[subd].c_str());
0457   std::map<unsigned int,RecJet::Hists> hh_ = MakeCorr(m_, subd, -1);
0458
0459   StoreCorrFactor(hh_, corrFactor_, clear);
0460   std::pair<double,double> meandev = MeanDeviation(sdet, corrFactor_);
0461   std::cout << "Mean deviation from 1: " << meandev.first << " and "
0462         << meandev.second << std::endl;
0463   sprintf(name, "Comparison%s", detector.c_str());
0464   sprintf(title,"Comparison of Correction Factor from Variance & Mean (%s)", detector.c_str());
0465   TH2D* h = new TH2D(name, title, 100, 0.5, 2.0, 100., 0.5, 2.0);
0466   for (std::map<unsigned int,RecJet::Hists>::const_iterator hitr = hh_.begin();
0467        hitr != hh_.end(); ++hitr) {
0468     if (hitr->second.h3 != 0 && hitr->second.h4 != 0) {
0469       for (int i=1; i<=hitr->second.h3->GetNbinsX(); ++i) {
0470     double x1 = hitr->second.h3->GetBinContent(i);
0471     double x2 = hitr->second.h4->GetBinContent(i);
0472         if (x1 != 0 && x2 != 0) h->Fill(x1,x2);
0473       }
0474     }
0475     if (hitr->second.h1 != 0) hitr->second.h1->Write();
0476     if (hitr->second.h2 != 0) hitr->second.h2->Write();
0477     if (hitr->second.h3 != 0) hitr->second.h3->Write();
0478     if (hitr->second.h4 != 0) hitr->second.h4->Write();
0479     if (hitr->second.h5 != 0) hitr->second.h5->Write();
0480     if (hitr->second.h6 != 0) hitr->second.h6->Write();
0481     if (hitr->second.h7 != 0) hitr->second.h7->Write();
0482     if (hitr->second.h8 != 0) hitr->second.h8->Write();
0483     if (hitr->second.h9 != 0) hitr->second.h9->Write();
0484
0485
0486   }
0487   h->Write();
0488   f.Close();
0489 }
0490
0491 void RecJet::LoopIter(int sdet, std::string indx, bool clear, int maxIter) {
0492
0493   std::map<unsigned int,RecJet::MyInfo> m_ = LoopMap();
0494   if (m_.size() == 0) return;
0495
0496   int subd = ((sdet>=1 && sdet<=4) ? (sdet-1) : 0);
0497   std::string det[4] = {"HB", "HE", "HO", "HF"};
0498   char fname[80], name[80], title[80];
0499   sprintf (fname, "Pileup%s_phiRecoDistribution_%s.root", indx.c_str(), det[subd].c_str());
0500   TFile f(fname, "recreate");
0501   detector = det[subd];
0502   sprintf (name, "Conv1%s", detector.c_str());
0503   sprintf (title, "Mean Deviation vs Iteration (%s)", detector.c_str());
0504   TH1D* h1 = new TH1D(name, title, maxIter+1, 0, maxIter+1);
0505   sprintf (name, "Conv2%s", detector.c_str());
0506   TH1D* h2 = new TH1D(name, title, maxIter+1, 0, maxIter+1);
0507   // 0th iteration and initiate correction factor
0508   std::map<unsigned int,RecJet::Hists> hh_ = MakeCorr(m_, subd, 0);
0509   StoreCorrFactor(hh_, corrFactor_, clear);
0510   std::pair<double,double> meandev = MeanDeviation(sdet, corrFactor_);
0511   h1->SetBinContent(0,meandev.first);
0512   h2->SetBinContent(0,meandev.second);
0513   std::cout << "Mean deviation from 1 in iteration 0: " << meandev.first
0514         << " and " << meandev.second << std::endl;
0515   for (std::map<unsigned int,RecJet::Hists>::const_iterator hitr = hh_.begin();
0516        hitr != hh_.end(); ++hitr) {
0517     if (hitr->second.h1 != 0) hitr->second.h1->Write();
0518     if (hitr->second.h2 != 0) hitr->second.h2->Write();
0519     if (hitr->second.h3 != 0) hitr->second.h3->Write();
0520     if (hitr->second.h4 != 0) hitr->second.h4->Write();
0521     if (hitr->second.h5 != 0) hitr->second.h5->Write();
0522     if (hitr->second.h6 != 0) hitr->second.h6->Write();
0523     if (hitr->second.h7 != 0) hitr->second.h7->Write();
0524     if (hitr->second.h8 != 0) hitr->second.h8->Write();
0525     if (hitr->second.h9 != 0) hitr->second.h9->Write();
0526
0527   }
0528   //Now iterate
0529   for (int kit = 1; kit <= maxIter; ++kit) {
0530     std::map<unsigned int,RecJet::MyInfo>   m2_  = MultiplyEnergy(m_,1);
0531     std::map<unsigned int,RecJet::Hists>    hh2_ = MakeCorr(m_, subd, kit);
0532     std::map<unsigned int,RecJet::CFactors> cF2_;
0533     StoreCorrFactor(hh2_, cF2_, true);
0534     meandev = MeanDeviation(sdet, cF2_);
0535     h1->SetBinContent(kit,meandev.first);
0536     h2->SetBinContent(kit,meandev.second);
0537     std::cout << "Mean deviation from 1 in iteration " << kit << ": "
0538           << meandev.first << " and " << meandev.second << std::endl;
0539     //Multiply cF2 to CorrFactor and make exit route if meandev is small
0540     for (std::map<unsigned int,RecJet::CFactors>::iterator it = corrFactor_.begin();
0541        it != corrFactor_.end(); ++it) {
0542       std::map<unsigned int,RecJet::CFactors>::iterator jt = cF2_.find(it->first);
0543       if (jt != cF2_.end()) {
0544     it->second.cfac1 *= (jt->second.cfac1);
0545     it->second.efac1 *= (jt->second.cfac1);
0546     it->second.cfac2 *= (jt->second.cfac2);
0547     it->second.efac2 *= (jt->second.cfac2);
0548       }
0549     }
0550     if (meandev.first < 0.0001 || kit==maxIter) {
0551       for (std::map<unsigned int,RecJet::Hists>::const_iterator hitr = hh2_.begin();
0552        hitr != hh2_.end(); ++hitr) {
0553     if (hitr->second.h1 != 0) hitr->second.h1->Write();
0554     if (hitr->second.h2 != 0) hitr->second.h2->Write();
0555     if (hitr->second.h3 != 0) hitr->second.h3->Write();
0556     if (hitr->second.h4 != 0) hitr->second.h4->Write();
0557     if (hitr->second.h5 != 0) hitr->second.h5->Write();
0558     if (hitr->second.h6 != 0) hitr->second.h6->Write();
0559     if (hitr->second.h7 != 0) hitr->second.h7->Write();
0560     if (hitr->second.h8 != 0) hitr->second.h8->Write();
0561     if (hitr->second.h9 != 0) hitr->second.h9->Write();
0562
0563       }
0564       break;
0565     }
0566   }
0567
0568   sprintf (name, "Comparison%s", detector.c_str());
0569   sprintf (title,"Comparison of Correction Factor from Variance & Mean (%s)", detector.c_str());
0570   TH2D* h = new TH2D(name, title, 100, 0.5, 2.0, 100, 0.5, 2.0);
0571   for (std::map<unsigned int,RecJet::CFactors>::iterator itr = corrFactor_.begin();
0572        itr != corrFactor_.end(); ++itr) {
0573     double x1 = itr->second.cfac1;
0574     double x2 = itr->second.cfac2;
0575     if (x1 != 0 && x2 != 0) h->Fill(x1,x2);
0576   }
0577   h->Write();
0578   h1->Write();
0579   h2->Write();
0580   f.Close();
0581 }
0582
0583     /*h1->distribution of mean energy
0584       h2->distribution of variance
0585       h3->correction factor from variance method
0586       h4->correction factor from mean method
0587       h5->corrected mean by variance method
0588       h6->corrected mean from mean method
0589       h7->corrected variance
0590       h8->corrected variance directly from the mean of variance
0591       h9->corrected mean directly from the mean of mean energy*/
0592
0593
0594 void RecJet::LoopIterate(int sdet, std::string indx, double emin, double emax, int maxIter) {
0595
0596   int subd = ((sdet>=1 && sdet<=4) ? (sdet-1) : 0);
0597   std::string det[4] = {"HB", "HE", "HO", "HF"};
0598   char fname[80], name[80], title[80];
0599   sprintf (fname, "Pileup%s_phiRecoDistribution_%s.root", indx.c_str(), det[subd].c_str());
0600   TFile f(fname, "recreate");
0601   detector = det[subd];
0602   sprintf (name, "Conv1%s", detector.c_str());
0603   sprintf (title, "Mean Deviation vs Iteration (%s)", detector.c_str());
0604   TH1D* h1 = new TH1D(name, title, maxIter+1, 0, maxIter+1);
0605   sprintf (name, "Conv2%s", detector.c_str());
0606   TH1D* h2 = new TH1D(name, title, maxIter+1, 0, maxIter+1);
0607   // 0th iteration and initiate correction factor
0608   std::map<unsigned int,RecJet::MyInfo> m_ = BuildMap(sdet, true, emin, emax);
0609   std::map<unsigned int,RecJet::Hists> hh_ = MakeCorr(m_, subd, 0);
0610   StoreCorrFactor(hh_, corrFactor_, true);
0611   std::pair<double,double> meandev = MeanDeviation(sdet, corrFactor_);
0612   h1->SetBinContent(0,meandev.first);
0613   h2->SetBinContent(0,meandev.second);
0614   std::cout << "Mean deviation from 1 in iteration 0: " << meandev.first
0615         << " and " << meandev.second << std::endl;
0616   for (std::map<unsigned int,RecJet::Hists>::const_iterator hitr = hh_.begin();
0617        hitr != hh_.end(); ++hitr) {
0618     if (hitr->second.h1 != 0) hitr->second.h1->Write();
0619     if (hitr->second.h2 != 0) hitr->second.h2->Write();
0620     if (hitr->second.h3 != 0) hitr->second.h3->Write();
0621     if (hitr->second.h4 != 0) hitr->second.h4->Write();
0622     if (hitr->second.h5 != 0) hitr->second.h5->Write();
0623     if (hitr->second.h6 != 0) hitr->second.h6->Write();
0624     if (hitr->second.h7 != 0) hitr->second.h7->Write();
0625     if (hitr->second.h8 != 0) hitr->second.h8->Write();
0626     if (hitr->second.h9 != 0) hitr->second.h9->Write();
0627
0628   }
0629   //Now iterate
0630   for (int kit = 1; kit <= maxIter; ++kit) {
0631     std::map<unsigned int,RecJet::MyInfo>   m2_  = BuildMap(sdet, false, emin, emax);
0632     std::map<unsigned int,RecJet::Hists>    hh2_ = MakeCorr(m_, subd, kit);
0633     std::map<unsigned int,RecJet::CFactors> cF2_;
0634     StoreCorrFactor(hh2_, cF2_, true);
0635     meandev = MeanDeviation(sdet, cF2_);
0636     h1->SetBinContent(kit,meandev.first);
0637     h2->SetBinContent(kit,meandev.second);
0638     std::cout << "Mean deviation from 1 in iteration " << kit << ": "
0639           << meandev.first << " and " << meandev.second << std::endl;
0640     //Multiply cF2 to CorrFactor and make exit route if meandev is small
0641     for (std::map<unsigned int,RecJet::CFactors>::iterator it = corrFactor_.begin();
0642        it != corrFactor_.end(); ++it) {
0643       std::map<unsigned int,RecJet::CFactors>::iterator jt = cF2_.find(it->first);
0644       if (jt != cF2_.end()) {
0645     it->second.cfac1 *= (jt->second.cfac1);
0646     it->second.efac1 *= (jt->second.cfac1);
0647     it->second.cfac2 *= (jt->second.cfac2);
0648     it->second.efac2 *= (jt->second.cfac2);
0649       }
0650     }
0651     if (meandev.first < 0.0001 || kit==maxIter) {
0652       for (std::map<unsigned int,RecJet::Hists>::const_iterator hitr = hh2_.begin();
0653        hitr != hh2_.end(); ++hitr) {
0654     if (hitr->second.h1 != 0) hitr->second.h1->Write();
0655     if (hitr->second.h2 != 0) hitr->second.h2->Write();
0656     if (hitr->second.h3 != 0) hitr->second.h3->Write();
0657     if (hitr->second.h4 != 0) hitr->second.h4->Write();
0658     if (hitr->second.h5 != 0) hitr->second.h5->Write();
0659     if (hitr->second.h6 != 0) hitr->second.h6->Write();
0660     if (hitr->second.h7 != 0) hitr->second.h7->Write();
0661     if (hitr->second.h8 != 0) hitr->second.h8->Write();
0662     if (hitr->second.h9 != 0) hitr->second.h9->Write();
0663
0664       }
0665       break;
0666     }
0667   }
0668
0669   sprintf (name, "Comparison%s", detector.c_str());
0670   sprintf (title,"Comparison of Correction Factor from Variance & Mean (%s)", detector.c_str());
0671   TH2D* h = new TH2D(name, title, 100, 0.5, 2.0, 100, 0.5, 2.0);
0672   for (std::map<unsigned int,RecJet::CFactors>::iterator itr = corrFactor_.begin();
0673        itr != corrFactor_.end(); ++itr) {
0674     double x1 = itr->second.cfac1;
0675     double x2 = itr->second.cfac2;
0676     if (x1 != 0 && x2 != 0) h->Fill(x1,x2);
0677   }
0678   h->Write();
0679   h1->Write();
0680   h2->Write();
0681   f.Close();
0682 }
0683
0684 std::map<unsigned int, RecJet::Hists>
0685 RecJet::MakeCorr(std::map<unsigned int,RecJet::MyInfo>& m_, int subd,
0686          int itr) {
0687
0688   char name[80], title[80];
0689   std::map<unsigned int, RecJet::Hists> hh_ ;
0690   for (std::map<unsigned int,RecJet::MyInfo>::iterator mitr = m_.begin();
0691        mitr !=m_.end(); mitr++) {
0692     unsigned int sdet0  = ((mitr->first)>>20)&0x7;
0693     if (sdet0 == (unsigned int)(subd+1)) {
0694       unsigned int detId2  = (mitr->first)&0x7fff80 ;
0695       std::map<unsigned int,RecJet::Hists>::iterator hitr = hh_.find(detId2);
0696       if (hitr == hh_.end()) {
0697     RecJet::Hists hh;
0698     hh_[detId2] = hh;
0699     hitr = hh_.find(detId2);
0700     int keta = (((mitr->first)&0x2000) ? (((mitr->first)>>7)&0x3f) :
0701             -(((mitr->first)>>7)&0x3f));
0702     int dept = (((mitr->first)>>14)&0x1f);
0703     if (itr >= 0) {
0704       sprintf (name,  "hist%s%d %d%d", detector.c_str(), keta, dept, itr);
0705       sprintf (title, "#phi (%s #eta %d depth %d iter %d)", detector.c_str(), keta, dept, itr);
0706     } else {
0707       sprintf (name,  "hist%s%d%d", detector.c_str(), keta, dept);
0708       sprintf (title, "#phi (%s #eta %d depth %d)", detector.c_str(), keta, dept);
0709     }
0710     hitr->second.h1 = new TH1D(name, title, 72, 0, 72.);
0711         hitr->second.h1->GetXaxis()->SetTitle("i#phi");
0712     hitr->second.h1->GetXaxis()->CenterTitle();
0713         hitr->second.h1->GetYaxis()->SetTitle("Mean Energy(GeV)");
0714     hitr->second.h1->GetYaxis()->CenterTitle();
0715     if (itr >= 0) {
0716       sprintf (name,  "histv%s%d %d%d", detector.c_str(), keta, dept, itr);
0717     } else {
0718       sprintf (name,  "histv%s%d%d", detector.c_str(), keta, dept);
0719     }
0720     hitr->second.h2 = new TH1D(name, title, 72, 0, 72.);
0721     hitr->second.h2->GetXaxis()->SetTitle("i#phi");
0722     hitr->second.h2->GetXaxis()->CenterTitle();
0723     hitr->second.h2->GetYaxis()->SetTitle("Variance (GeV^2)");
0724     hitr->second.h2->GetYaxis()->CenterTitle();
0725       }
0726       std::pair<double,double> mean, variance;
0727       MeanVariance(mitr,mean,variance);
0728       std::pair<double,double> nsub = SubtractNoise(mitr->first,mean.first,
0729                             mean.second,false);
0730       int kphi = ((mitr->first)&0x7f);
0731       hitr->second.h1->SetBinContent(kphi,nsub.first);
0732       hitr->second.h1->SetBinError(kphi,nsub.second);
0733       nsub = SubtractNoise(mitr->first,variance.first,
0734                variance.second,true);
0735       hitr->second.h2->SetBinContent(kphi,variance.first);
0736       hitr->second.h2->SetBinError(kphi,variance.second);
0737     }
0738   }
0739   for (std::map<unsigned int, RecJet::Hists>::iterator hitr = hh_.begin();
0740        hitr != hh_.end(); ++hitr) {
0741     int keta = (((hitr->first)&0x2000) ? (((hitr->first)>>7)&0x3f) :
0742         -(((hitr->first)>>7)&0x3f));
0743     int dept = (((hitr->first)>>14)&0x1f);
0744     if (itr >= 0) {
0745       sprintf (name,  "correctionMean%s%d %d%d", detector.c_str(), keta, dept, itr);
0746       sprintf (title, "Correction Factor (Mean) for (%s #eta %d depth %d iter %d)", detector.c_str(), keta, dept, itr);
0747     } else {
0748       sprintf (name,  "correctionMean%s%d%d", detector.c_str(), keta, dept);
0749       sprintf (title, "Correction Factor (Mean) for (%s #eta %d depth %d)", detector.c_str(), keta, dept);
0750     }
0751     hitr->second.h3 = new TH1D(name, title, 72, 0., 72.);
0752     hitr->second.h3->GetXaxis()->SetTitle("i#phi");
0753     hitr->second.h3->GetXaxis()->CenterTitle();
0754     hitr->second.h3->GetYaxis()->SetTitle("Correction Factor");
0755     hitr->second.h3->GetYaxis()->CenterTitle();
0756     if (itr >= 0) {
0757       sprintf (name,  "correctionVariance%s%d %d%d", detector.c_str(), keta, dept, itr);
0758       sprintf (title, "Correction Factor (Variance) for (%s #eta %d depth %d iter %d)", detector.c_str(), keta, dept, itr);
0759     } else {
0760       sprintf (name,  "correctionVariance%s%d%d", detector.c_str(), keta, dept);
0761       sprintf (title, "Correction Factor (Variance) for (%s #eta %d depth %d)", detector.c_str(), keta, dept);
0762     }
0763     hitr->second.h4 = new TH1D(name, title, 72, 0., 72.);
0764     hitr->second.h4->GetXaxis()->SetTitle("i#phi");
0765     hitr->second.h4->GetXaxis()->CenterTitle();
0766     hitr->second.h4->GetYaxis()->SetTitle("Correction Factor");
0767     hitr->second.h4->GetYaxis()->CenterTitle();
0768     sprintf (name,  "Corrected_Mean_var%s%d %d%d", detector.c_str(), keta, dept,itr);
0769     sprintf (title, "Corrected_Mean (by Variance) for (%s #eta %d depth %d) %d", detector.c_str(), keta, dept,itr);
0770     hitr->second.h5 = new TH1D(name, title, 72, 0., 72.);
0771     hitr->second.h5->GetXaxis()->SetTitle("i#phi");
0772     hitr->second.h5->GetXaxis()->CenterTitle();
0773     hitr->second.h5->GetYaxis()->SetTitle("Corrected_Mean");
0774     hitr->second.h5->GetYaxis()->CenterTitle();
0775     sprintf (name,  "Corrected_Mean_mean%s%d %d%d", detector.c_str(), keta, dept,itr);
0776     sprintf (title, "Corrected_Mean (by Mean) for (%s #eta %d depth %d) %d", detector.c_str(), keta, dept,itr);
0777     hitr->second.h6 = new TH1D(name, title, 72, 0., 72.);
0778     hitr->second.h6->GetXaxis()->SetTitle("i#phi");
0779     hitr->second.h6->GetXaxis()->CenterTitle();
0780     hitr->second.h6->GetYaxis()->SetTitle("Corrected_Mean");
0781     hitr->second.h6->GetYaxis()->CenterTitle();
0782     sprintf (name,  "Corrected_Variance%s%d %d%d", detector.c_str(), keta, dept,itr);
0783     sprintf (title, "Corrected_Variance (by Variance) for (%s #eta %d depth %d) %d", detector.c_str(), keta, dept,itr);
0784     hitr->second.h7 = new TH1D(name, title, 72, 0., 72.);
0785     hitr->second.h7->GetXaxis()->SetTitle("i#phi");
0786     hitr->second.h7->GetXaxis()->CenterTitle();
0787     hitr->second.h7->GetYaxis()->SetTitle("Corrected_Var");
0788     hitr->second.h7->GetYaxis()->CenterTitle();
0789     sprintf (name,  "Corrected_Variance_(my method)%s%d %d%d", detector.c_str(), keta, dept,itr);
0790     sprintf (title, "Corrected_Variance (by Variance) for (%s #eta %d depth %d) %d", detector.c_str(), keta, dept,itr);
0791     hitr->second.h8 = new TH1D(name, title, 72, 0., 72.);
0792     hitr->second.h8->GetXaxis()->SetTitle("i#phi");
0793     hitr->second.h8->GetXaxis()->CenterTitle();
0794     hitr->second.h8->GetYaxis()->SetTitle("Corrected_Var");
0795     hitr->second.h8->GetYaxis()->CenterTitle();
0796     sprintf (name,  "Corrected_Mean_(my method)%s%d %d%d", detector.c_str(), keta, dept,itr);
0797     sprintf (title, "Corrected_Variance (by Variance) for (%s #eta %d depth %d) %d", detector.c_str(), keta, dept,itr);
0798     hitr->second.h9 = new TH1D(name, title, 72, 0., 72.);
0799     hitr->second.h9->GetXaxis()->SetTitle("i#phi");
0800     hitr->second.h9->GetXaxis()->CenterTitle();
0801     hitr->second.h9->GetYaxis()->SetTitle("Corrected_Mean");
0802     hitr->second.h9->GetYaxis()->CenterTitle();
0803
0804     std::vector<std::pair<double,double> > corr1 = CorrFactor(subd,hitr,false);
0805     for (unsigned int i=0; i<corr1.size(); ++i) {
0806       if (corr1[i].first != 0 && corr1[i].second !=0){
0807     hitr->second.h3->SetBinContent(i+1, corr1[i].first);
0808     hitr->second.h3->SetBinError(i+1,   corr1[i].second);
0809     double bin_cont = hitr->second.h1->GetBinContent(i+1) * corr1[i].first;
0810         double efac1 = (hitr->second.h1->GetBinError(i+1))*(corr1[i].first);
0811         double efac2 = (hitr->second.h1->GetBinContent(i+1))*(corr1[i].second);
0812         double bin_err = std::sqrt(efac1*efac1+efac2*efac2);
0813         if (bin_cont !=0) {
0814       hitr->second.h6->SetBinContent(i+1, bin_cont);
0815       hitr->second.h6->SetBinError(i+1, bin_err);
0816     }
0817       }
0818     }
0819     std::vector<std::pair<double,double> > corr2 = CorrFactor(subd,hitr, true);
0820     for (unsigned int i=0; i<corr2.size(); ++i) {
0821       if (corr2[i].first != 0 && corr2[i].second !=0) {
0822     double var = hitr->second.h2->GetBinContent(i+1)* corr2[i].first * corr2[i].first;
0823     hitr->second.h7->SetBinContent(i+1,var);
0824     double var_err1 = corr2[i].first * corr2[i].first * hitr->second.h2->GetBinError(i+1);
0825     double var_err2 = 2. * corr2[i].first * corr2[i].second * hitr->second.h2->GetBinContent(i+1);
0826     double var_err = std::sqrt(var_err1*var_err1 + var_err2*var_err2);
0827     hitr->second.h7->SetBinError(i+1,var_err);
0828     hitr->second.h4->SetBinContent(i+1, corr2[i].first);
0829     hitr->second.h4->SetBinError(i+1,   corr2[i].second);
0830     double bin_cont = hitr->second.h1->GetBinContent(i+1) * corr2[i].first;
0831         if(bin_cont !=0) hitr->second.h5->SetBinContent(i+1, bin_cont);
0832     double efac1 = (hitr->second.h1->GetBinError(i+1))*(corr2[i].first);
0833     double efac2 = (hitr->second.h1->GetBinContent(i+1))*(corr2[i].second);
0834     double bin_err = std::sqrt(efac1*efac1+efac2*efac2);
0835         if(bin_cont !=0) hitr->second.h5->SetBinError(i+1, bin_err);
0836       }
0837     }
0838   }
0839   return hh_;
0840 }
0841
0842 void RecJet::LoopMean(int sdet, std::string indx) {
0843   std::map<unsigned int,RecJet::MyInfo> m_ = LoopMap();
0844   if (m_.size() == 0) return;
0845   int subd = ((sdet>=1 && sdet<=4) ? (sdet-1) : 0);
0846   std::string det[4] = {"HB", "HE", "HO", "HF"};
0847   char fname[80],name[80],title[80];
0848   sprintf (fname, "True_etaRecoDistribution%s_%s.root", indx.c_str(), det[subd].c_str());
0849   TFile f(fname, "recreate");
0850   detector = det[subd];
0851
0852   std::map<unsigned int,std::pair<TH1D*,TH1D*> > myMap;
0853   for (std::map<unsigned int,RecJet::MyInfo>::iterator mitr = m_.begin();
0854        mitr !=m_.end(); mitr++) {
0855     unsigned int sdet0  = ((mitr->first)>>20)&0x7;
0856     if (sdet0 == (unsigned int)(subd+1)) {
0857       unsigned int detId2 =  (mitr->first)&0x7C000 ;
0858       std::map<unsigned int, std::pair<TH1D*,TH1D*> >::iterator hitr = myMap.find(detId2);
0859       int b2 = (subd == 3) ? 42 : ((subd == 1) ? 30 : 17);
0860       if (hitr == myMap.end()) {
0861     int dept = (((mitr->first)>>14)&0x1f);
0862     sprintf (name,  "hist%s%d", detector.c_str(), dept);
0863     sprintf (title, "Mean Response (%s depth %d)", detector.c_str(), dept);
0864     TH1D* h1 = new TH1D(name, title, 2*b2, -(double)(b2), (double)(b2));
0865     h1->GetXaxis()->SetTitle("i#eta"); h1->GetXaxis()->CenterTitle();
0866     sprintf (name,  "histv%s%d", detector.c_str(), dept);
0867     TH1D* h2 = new TH1D(name, title, 2*b2, -(double)(b2), (double)(b2));
0868     h2->GetXaxis()->SetTitle("i#eta"); h2->GetXaxis()->CenterTitle();
0869     if (subd == 3) {
0870       h1->GetYaxis()->SetTitle("Mean # PE"); h2->GetYaxis()->SetTitle("Variace(# PE)");
0871     } else {
0872       h1->GetYaxis()->SetTitle("Mean Energy (GeV)"); h2->GetYaxis()->SetTitle("Variance (GeV^2)");
0873     }
0874     myMap[detId2] = std::pair<TH1D*,TH1D*>(h1,h2);
0875     hitr = myMap.find(detId2);
0876       }
0877       int keta        = (((mitr->first)&0x2000) ? (((mitr->first)>>7)&0x3f) :
0878              -(((mitr->first)>>7)&0x3f));
0879       int bin         = keta + b2 + 1;
0880       std::pair<double,double> mean, variance;
0881       MeanVariance(mitr,mean,variance);
0882       std::pair<double,double> nsub = SubtractNoise(mitr->first,mean.first,mean.second,false);
0883       hitr->second.first->SetBinContent(bin,nsub.first);
0884       hitr->second.first->SetBinError(bin,nsub.second);
0885       nsub = SubtractNoise(mitr->first,variance.first,variance.second,true);
0886       hitr->second.second->SetBinContent(bin,nsub.first);
0887       hitr->second.second->SetBinError(bin,nsub.second);
0888     }
0889   }
0890   for (std::map<unsigned int, std::pair<TH1D*,TH1D*> >::iterator hitr = myMap.begin();
0891        hitr != myMap.end(); ++hitr) {
0892     hitr->second.first->Write(); hitr->second.second->Write();
0893   }
0894   f.Close();
0895 }
0896
0897 void RecJet::Alleta_distribution(std::string fname, bool var, bool HBHE) {
0898   std::map<unsigned int,RecJet::MyInfo> m_ = LoopMap();
0899   TFile f(fname.c_str(), "recreate");
0900   if (m_.size() == 0) return;
0901   TH1D* h;
0902   char name[50], title[50];
0903   if (var) {
0904     sprintf(name,"Alleta");
0905     sprintf(title,"Variance distribution w.r.t.eta");
0906   } else {
0907     sprintf(name,"Alleta");
0908     sprintf(title,"Energy distribution w.r.t.eta");
0909   }
0910   h = new TH1D(name,title,82, -41., 41.);
0911   if (var) {
0912     h->GetYaxis()->SetTitle("Variance in Gev^2"); h->GetYaxis()->CenterTitle();
0913   } else {
0914     h->GetYaxis()->SetTitle("Mean Energy in Gev"); h->GetYaxis()->CenterTitle();
0915   }
0916   h->GetXaxis()->SetTitle("i#eta"); h->GetXaxis()->CenterTitle();
0917
0918
0919   std::map<int, std::pair<double,double> > h_eta;
0920   for (std::map<unsigned int,RecJet::MyInfo>::iterator mitr = m_.begin();
0921        mitr !=m_.end(); mitr++) {
0922     int eta = (mitr->first)&0x2000 ? ((mitr->first)>>7)&0x3f : -(((mitr->first)>>7)&0x3f);
0923     std::map<int,std::pair<double,double> >::iterator hitr = h_eta.find(eta);
0924     if (hitr == h_eta.end()) {
0925       double energy=0.;
0926       double error =0.;
0927       h_eta[eta] = std::pair<double,double>(energy,error);
0928       hitr          = h_eta.find(eta);
0929     }
0930     std::pair<double,double> mean, variance;
0931     MeanVariance(mitr,mean,variance);
0932     std::pair<double, double> nsub;
0933     if (var) {
0934       nsub = SubtractNoise(mitr->first,variance.first,variance.second,true);
0935     } else {
0936      nsub = SubtractNoise(mitr->first,mean.first,mean.second,false);
0937     }
0938     hitr->second.first +=nsub.first;
0939     hitr->second.second +=pow(nsub.second,2);
0940   }
0941   for(std::map<int,std::pair<double,double> >::iterator hitr = h_eta.begin();
0942       hitr !=h_eta.end(); hitr++) {
0943     int eta = hitr->first;
0944     if(HBHE && fabs(eta)>29) continue;
0945     if(!HBHE && fabs(eta) <=29) continue;
0946     h->SetBinContent(40+eta+1,hitr->second.first);
0947     h->SetBinError(40+eta+1,std::sqrt(hitr->second.second));
0948   }
0949   h->Write();
0950   f.Close();
0951 }
0952
0953
0954 void RecJet::eta_distribution(std::string fname, bool var) {
0955   std::map<unsigned int,RecJet::MyInfo> m_ = LoopMap();
0956   char name[50], title[50];
0957   TFile f(fname.c_str(), "recreate");
0958   if (m_.size() == 0) return;
0959   std::map<unsigned int, TH1D*> h;
0960   std::string det[4] = {"HB", "HE", "HO", "HF"};
0961   for(int idet=1; idet<=4; idet++){
0962     for(int depth=1; depth<=3; depth++){
0963       unsigned int detId = (((idet&0x7)<<5)|(depth&0x1f));
0964       if (var) {
0965     sprintf(name,"%s%d_eta",det[idet-1].c_str(),depth);
0966     sprintf(title,"Variance distribution w.r.t.eta for %s depth %d", det[idet-1].c_str(),depth);
0967       } else {
0968     sprintf(name,"%s%d_eta",det[idet-1].c_str(),depth);
0969         sprintf(title,"Energy distribution w.r.t.eta for %s depth %d", det[idet-1].c_str(),depth);
0970       }
0971       TH1D* hist = new TH1D(name,title,82, -41., 41.);
0972       if (var) {
0973     hist->GetYaxis()->SetTitle("Variance in Gev^2"); hist->GetYaxis()->CenterTitle();
0974       } else {
0975     hist->GetYaxis()->SetTitle("Mean Energy in Gev"); hist->GetYaxis()->CenterTitle();
0976       }
0977       hist->GetXaxis()->SetTitle("i#eta"); hist->GetXaxis()->CenterTitle();
0978       h[detId] = hist;
0979     }
0980   }
0981
0982   std::map<unsigned int, std::pair<double,double> > h_eta;
0983   for (std::map<unsigned int,RecJet::MyInfo>::iterator mitr = m_.begin();
0984        mitr !=m_.end(); mitr++) {
0985     int subd = ((mitr->first)>>20)&0x7;
0986     int eta = (mitr->first)&0x2000 ? ((mitr->first)>>7)&0x3f : -(((mitr->first)>>7)&0x3f);
0987     int depth = ((mitr->first)>>14)&0x1f;
0988     unsigned int detId = (((subd&0x7)<<12) | ((depth&0x1f)<<7) | ((eta>0) ? (0x40|eta) : -eta));
0989     std::map<unsigned int,std::pair<double,double> >::iterator hitr = h_eta.find(detId);
0990     if (hitr == h_eta.end()) {
0991       double energy=0.;
0992       double error =0.;
0993       h_eta[detId] = std::pair<double,double>(energy,error);
0994       hitr          = h_eta.find(detId);
0995       }
0996     std::pair<double,double> mean, variance;
0997     MeanVariance(mitr,mean,variance);
0998     std::pair<double,double> nsub;
0999     if (var){
1000       nsub = SubtractNoise(mitr->first,variance.first,variance.second,true);
1001     } else {
1002       nsub = SubtractNoise(mitr->first,mean.first,mean.second,false);
1003     }
1004     hitr->second.first +=nsub.first;
1005     hitr->second.second +=pow(nsub.second,2);
1006   }
1007   for(std::map<unsigned int,std::pair<double,double> >::iterator hitr = h_eta.begin();
1008       hitr !=h_eta.end(); hitr++) {
1009     unsigned int detId = (((hitr->first)>>7)&0xff);
1010     int eta = (((hitr->first)&0x40) ? ((hitr->first)&0x3f) : -((hitr->first)&0x3f));
1011     for(std::map<unsigned int,TH1D* >::iterator itr = h.begin();
1012     itr !=h.end(); itr++){
1013       if(itr->first != detId)continue;
1014       itr->second->SetBinContent(40+eta+1,hitr->second.first);
1015       itr->second->SetBinError(40+eta+1,std::sqrt(hitr->second.second));
1016     }
1017   }
1018   for(std::map<unsigned int,TH1D* >::iterator itr = h.begin();
1019       itr !=h.end(); itr++) if(itr->second->GetEntries() !=0) itr->second->Write();
1020   f.Close();
1021 }
1022
1023 void RecJet::det_distribution(std::string fname) {
1024   TFile f(fname.c_str(),"recreate");
1025   std::map<unsigned int,RecJet::MyInfo> m_ = LoopMap();
1026   if (m_.size() == 0) return;
1027   TH1D* h[4];
1028   std::string det[4] = {"HB", "HE","HO", "HF"};
1029   char        name[700], title[700];
1030   for (int idet=1; idet<=4; idet++) {
1031     sprintf(name, "%s", det[idet-1].c_str());
1032     sprintf (title, "Noise distribution for %s", det[idet-1].c_str());
1033     h[idet-1] = new TH1D(name,title,32,-4., 4.);
1034   }
1035   for (std::map<unsigned int,RecJet::MyInfo>::iterator mitr = m_.begin();
1036        mitr !=m_.end(); mitr++) {
1037     int subd = ((mitr->first)>>20)&0x7;
1038     std::pair<double,double> mean, variance;
1039     MeanVariance(mitr,mean,variance);
1040     h[subd-1]->Fill(mean.first);
1041   }
1042   for (int i=0; i<4; i++) h[i]->Write();
1043 }
1044
1045 void RecJet::Fit(std::string rootfile, std::string textfile) {
1046   TFile filex(rootfile.c_str());
1047   ofstream myfile;
1048   myfile.open(textfile.c_str());
1049   myfile << std::setw(25) << "Histo" << std::setw(8) << "Entries"
1050      << std::setw(10) << "Mean" << std::setw(8) << "RMS" << std::setw(20)
1051      << "Const_fit/Error" << std::setw(20) << "Mean_fit/Error"
1052      << std::setw(20) << "Sigma_fit/Error" << std::setw(10)
1053          << "Chi2/ndf" << "\n" <<std::endl;
1054   std::string det[3] = {"HB", "HE", "HF"};
1055   for(int idet=0; idet<=2; idet++){
1056     for(int ietax=-41; ietax<=41; ietax++){
1057       for(int iphix=1; iphix<=72; iphix++){
1058         for(int dep=1; dep<=3; dep++){
1059       std::ostringstream s;
1060           s << det[idet].c_str() << "eta" << ietax << "phi" << iphix << "dep" << dep;
1061       std::string d = s.str();
1062           TH1D* h; filex.GetObject(d.c_str(),h);
1063           if(h && h->GetEntries()) {
1064             myfile << std::setw(25)  << h->GetTitle() << std::setw(8)
1065            << h->GetEntries() << std::setw(10) << std::setprecision(4)
1066            << h->GetMean() << std::setw(8) << std::setprecision(4)
1067            << h->GetRMS();
1068             TF1 *f1 = new TF1("f1","[0]*TMath::Exp((-0.5/pow([2],2))*pow((x[0] -[1]),2))",-3,3);
1069             f1->SetParameters(500,h->GetMean(),h->GetRMS());
1070             f1->SetParNames ("Constant","Mean_value","Sigma");
1071             h->Fit("f1","R");
1072             double par[3];
1073             f1->GetParameters(&par[0]);
1074             myfile << std::setw(10) << std::setprecision(4) << par[0]
1075            << std::setw(9) << std::setprecision(4) << f1->GetParError(0)
1076                    << std::setw(10) << std::setprecision(4) << par[1]
1077            << std::setw(9) << std::setprecision(4) << f1->GetParError(1)
1078                    << std::setw(10) << std::setprecision(4) << par[2]
1079            << std::setw(9) << std::setprecision(4) <<f1->GetParError(2)
1080                    <<std::setw(10) << f1->GetChisquare()/f1->GetNDF() << endl;
1081           }
1082         }
1083       }
1084     }
1085   }
1086 }
1087
1088
1089 void RecJet::Disturb(int subd) {
1090   std::cout << "Factor of multiplication: ";
1091   std::cin >> factor_;
1092   dets_.clear();
1093   std::cout << "Enter detid for subd=" << subd << std::endl;
1094   int detId[3];
1095   std::string detid[3] = {"depth", "ieta", "iphi"};
1096   while (1) {
1097     bool ok(true);
1098     for (int i=0; i<3; i++) {
1099       std::cout << "Enter the value of " << detid[i].c_str() << std::endl;
1100       std::cin >> detId[i];
1101       if (i != 1 && detId[i] <= 0) {
1102     ok = false; break;
1103       }
1104     }
1105     if (ok) {
1106       unsigned int detId1 = ((subd<<20) | ((detId[0]&0x1f)<<14) |
1107                  ((detId[1]>0) ? (0x2000|(detId[1]<<7)) :
1108                   ((-detId[1])<<7)) | (detId[2]&0x7f));
1109       dets_.push_back(detId1);
1110     } else {
1111       break;
1112     }
1113   }
1114 }
1115
1116 std::map<unsigned int,RecJet::MyInfo> RecJet::MultiplyEnergy(std::map<unsigned int, RecJet::MyInfo>& m_, int corrfac) {
1117
1118   if (corrfac > 0) {
1119     for (std::map<unsigned int,RecJet::MyInfo>::iterator mitr = m_.begin();
1120      mitr !=m_.end(); mitr++) {
1121       std::map<unsigned int,RecJet::CFactors>::iterator itr = corrFactor_.find(mitr->first);
1122       double fac = ((itr == corrFactor_.end()) ? 1 :
1123             ((corrfac == 1) ? itr->second.cfac1 : itr->second.cfac2));
1124       ChangeMoments(mitr,fac);
1125     }
1126   } else {
1127     for (unsigned k=0; k<dets_.size(); ++k) {
1128       std::map<unsigned int,RecJet::MyInfo>::iterator mitr = m_.find(dets_[k]);
1129       if (mitr != m_.end()) ChangeMoments(mitr,factor_);
1130     }
1131   }
1132   return m_;
1133 }
1134
1135 void RecJet::ChangeMoments(std::map<unsigned int,RecJet::MyInfo>::iterator& mitr, double fac) {
1136
1137   mitr->second.Mom1 *= fac;
1138   mitr->second.Mom2 *= (fac*fac);
1139   mitr->second.Mom3 *= (fac*fac*fac);
1140   mitr->second.Mom4 *= (fac*fac*fac*fac);
1141 }
1142
1143 std::pair<double,double> RecJet::SubtractNoise(unsigned int detId, double mean,
1144                            double error, bool ifVar) {
1145
1146   std::map<unsigned int,RecJet::MyInfo>::iterator mitr = mNoise_.find(detId);
1147   if (mitr != mNoise_.end()) {
1148     std::cout << "//////////////////enter" << std::endl;
1149     double m2(0), e2(0);
1150     if (ifVar) {
1151       double mean1    = mitr->second.Mom1 / mitr->second.Mom0;
1152       double mean2    = mitr->second.Mom2 / mitr->second.Mom0;
1153       double mean3    = mitr->second.Mom3 / mitr->second.Mom0;
1154       double mean4    = mitr->second.Mom4 / mitr->second.Mom0;
1155       m2 = (mean2 - mean1*mean1);
1156       e2 = ((mean4+8*mean1*mean1*mean2-4*mean1*mean3-4*mean1*mean1*mean1*mean1-
1157          mean2*mean2)/mitr->second.Mom0);
1158     } else {
1159       m2 = mitr->second.Mom1 / mitr->second.Mom0 ;
1160       e2 = ((mitr->second.Mom2/mitr->second.Mom0)-m2*m2)/mitr->second.Mom0;
1161     }
1162     mean -= m2;
1163     error = std::sqrt(error*error+e2);
1164   }
1165   return std::pair<double,double>(mean,error);
1166 }
1167
1168 std::map<unsigned int,RecJet::MyInfo> RecJet::LoadNoise() {
1169
1170   std::map<unsigned int,RecJet::MyInfo> m_;
1171   std::string det[4] = {"HB", "HE", "HO", "HF"};
1172   for (int subd=0; subd<4; ++subd) {
1173     char fname[80];
1174     sprintf (fname,"%s_Noise.txt",det[subd].c_str());
1175     ifstream myfile(fname);
1176     if (!myfile.is_open()) {
1177       std::cerr << "** ERROR: Can't open '" << fname << "' for the noise file"
1178         << std::endl;
1179     } else {
1180       unsigned int nrec(0), ndets(0);
1181       while(1) {
1182     unsigned int detId;
1183     double       mom0, mom1, mom2, mom3, mom4;
1184     myfile >> detId >> mom0 >> mom1 >> mom2 >> mom3 >> mom4;
1185     if (!myfile.good()) break;
1186     std::map<unsigned int,RecJet::MyInfo>::iterator mitr = m_.find(detId);
1187     if (mitr == m_.end()) {
1188       RecJet::MyInfo info;
1189       m_[detId] = info ;
1190       mitr = m_.find(detId) ;
1191       ndets++;
1192     }
1193     mitr->second.Mom0 += mom0_MB;
1194     mitr->second.Mom1 += mom1_MB;
1195     mitr->second.Mom2 += mom2_MB;
1196     mitr->second.Mom3 += mom3_MB;
1197     mitr->second.Mom4 += mom4_MB;
1198     nrec++;
1199       }
1200       myfile.close();
1201       std::cout << "Reads " << nrec << " records for " << ndets << " elements"
1202         << " in file " << fname << std::endl;
1203     }
1204   }
1205   return m_;
1206 }
1207
1208 void RecJet::MeanVariance(std::map<unsigned int,RecJet::MyInfo>::iterator &mitr,
1209               std::pair<double,double>& mean,
1210               std::pair<double,double>& variance) {
1211   mean.first     = mitr->second.Mom1 / mitr->second.Mom0;
1212   double mean2   = mitr->second.Mom2 / mitr->second.Mom0;
1213   double mean3   = mitr->second.Mom3 / mitr->second.Mom0;
1214   double mean4   = mitr->second.Mom4 / mitr->second.Mom0;
1215   variance.first = (mean2 - mean.first*mean.first);
1216   mean.second    = std::sqrt(variance.first / mitr->second.Mom0);
1217   variance.second= std::sqrt((mean4+8*mean.first*mean.first*mean2-
1218                   4*mean.first*mean3-mean2*mean2-
1219                   4*mean.first*mean.first*mean.first*mean.first)/
1220                  mitr->second.Mom0);
1221 }
1222
1223 /* Calculation of correction factor following staggered geometry which is not requird*/
1224
1225
1226  std::vector<std::pair<double,double> > RecJet::Staggered_CorrFactor(int subd, std::map<unsigned int,RecJet::Hists>::iterator &hitr, bool varmethod) {
1227   std::vector<std::pair<double,double> > corrfac;
1228   double       mean1(0.), mean2(0.), error1(0.), error2(0.);
1229   unsigned int i_in = (subd == 0) ? 0 : 1;
1230   unsigned int i_end = (subd == 0) ? 71 : 72 ;
1231   if (varmethod) {
1232     for (unsigned int i=i_in; i<i_end; i=i+4) {
1233       mean1 += (i_in == 0) ? hitr->second.h2->GetBinContent(1)   + hitr->second.h2->GetBinContent(72) :
1234                      hitr->second.h2->GetBinContent(i)   + hitr->second.h2->GetBinContent(i+1) ;
1235       mean2 +=               hitr->second.h2->GetBinContent(i+2) + hitr->second.h2->GetBinContent(i+3) ;
1236
1237       error1 += (i_in == 0) ? pow((hitr->second.h2->GetBinError(1)), 2) + pow((hitr->second.h2->GetBinError(72)),2) :
1238                   pow((hitr->second.h2->GetBinError(i)), 2) + pow((hitr->second.h2->GetBinError(i+1)),2) ;
1239       error2 +=               pow((hitr->second.h2->GetBinError(i+2)),2)+ pow((hitr->second.h2->GetBinError(i+3)),2) ;
1240     }
1241     mean1  /= (hitr->second.h2->GetEntries()/2);
1242     mean2  /= (hitr->second.h2->GetEntries()/2);
1243     error1 /= (hitr->second.h2->GetEntries()/2);
1244     error2 /= (hitr->second.h2->GetEntries()/2);
1245   } else {
1246     for (unsigned int i=i_in; i<i_end; i=i+4) {
1247         mean1 += (i_in == 0) ? hitr->second.h1->GetBinContent(1)   + hitr->second.h1->GetBinContent(72) :
1248                    hitr->second.h1->GetBinContent(i)   + hitr->second.h1->GetBinContent(i+1) ;
1249         mean2 +=               hitr->second.h1->GetBinContent(i+2) + hitr->second.h1->GetBinContent(i+3) ;
1250
1251        error1 += (i_in == 0) ? pow((hitr->second.h1->GetBinError(1)), 2) + pow((hitr->second.h1->GetBinError(72)),2) :
1252                    pow((hitr->second.h1->GetBinError(i)), 2) + pow((hitr->second.h1->GetBinError(i+1)),2) ;
1253        error2 +=               pow((hitr->second.h1->GetBinError(i+2)),2)+ pow((hitr->second.h1->GetBinError(i+3)),2) ;
1254
1255     }
1256     mean1  /= (hitr->second.h1->GetEntries()/2);
1257     mean2  /= (hitr->second.h1->GetEntries()/2);
1258     error1 /= (hitr->second.h1->GetEntries()/2);
1259     error2 /= (hitr->second.h1->GetEntries()/2);
1260   }
1261
1262   double n     = (varmethod) ? hitr->second.h2->GetEntries() : hitr->second.h1->GetEntries();
1263   for (int i=1; i<=hitr->second.h1->GetNbinsX(); ++i) {
1264     double  mean =    (subd == 0) ? (((i-1)%4 == 0) || (i%4 == 0) ? mean1 : mean2) :
1265               (subd == 1) ? (((i+1)%4 == 0) || (i%4 == 0) ? mean2 : mean1) : (mean1 + mean2)/2 ;
1266
1267    double error =     (subd == 0) ? (((i-1)%4 == 0) || (i%4 == 0) ? error1 : error2) :
1268               (subd == 1) ? (((i+1)%4 == 0) || (i%4 == 0) ? error2 : error1) : (error2 + error1)/2;
1269    if(varmethod) {
1270      hitr->second.h8->SetBinContent(i,mean);
1271      hitr->second.h8->SetBinError(i,error);
1272      //     std::cout << error << std::endl;
1273    } else {
1274      hitr->second.h9->SetBinContent(i,mean);
1275      hitr->second.h9->SetBinError(i,error);
1276      std::cout << error << std::endl;
1277    }
1278
1279    double vcont = (varmethod) ? hitr->second.h2->GetBinContent(i) : hitr->second.h1->GetBinContent(i);
1280    double econt = (varmethod) ? hitr->second.h2->GetBinError(i)   : hitr->second.h1->GetBinError(i);
1281    double vcfac = (varmethod) ? ((vcont > 0) ? (std::sqrt(mean/vcont)) : 1) : ((vcont > 0) ? (mean/vcont) : 1);
1282    double ecfac = (varmethod) ? ((vcont>0 && mean>0) ? (0.5*std::sqrt(((2*error)/(vcont*mean*n)) + (econt*econt)/(vcont*vcont) *((mean/vcont) - (4.0/n)))) : 0) : ((vcont>0 && mean>0) ? (std::sqrt(((error*2.)/(n*vcont*vcont)) + (econt/vcont)*(econt/vcont)*((mean/vcont)*(mean/vcont)-(4.0/n)*(mean/vcont)))) : 0);
1283     corrfac.push_back(std::pair<double,double>(vcfac,ecfac));
1284   }
1285   return corrfac;
1286 }
1287
1288
1289 /* Calculation of correction factor where staggered geometry is not followed. This should be followed not previous one*/
1290
1291 std::vector<std::pair<double,double> > RecJet::CorrFactor(int , std::map<unsigned int,RecJet::Hists>::iterator &hitr, bool varmethod) {
1292   std::vector<std::pair<double,double> > corrfac;
1293   double       mean(0.), error(0.);
1294   if (varmethod) {
1295     for (unsigned int i=1; i<=72; i++) {
1296       mean  +=  hitr->second.h2->GetBinContent(i);
1297       error +=  pow((hitr->second.h2->GetBinError(i)), 2);
1298     }
1299     mean  /= hitr->second.h2->GetEntries();
1300     error /= hitr->second.h2->GetEntries();
1301   } else {
1302     for (unsigned int i=1; i<=72; i++) {
1303       mean  += hitr->second.h1->GetBinContent(i);
1304       error += pow((hitr->second.h1->GetBinError(i)), 2) ;
1305     }
1306     mean  /= hitr->second.h1->GetEntries();
1307     error /= hitr->second.h1->GetEntries();
1308   }
1309
1310   for (int i=1; i<=hitr->second.h1->GetNbinsX(); ++i) {
1311     if(varmethod) {
1312       hitr->second.h8->SetBinContent(i,mean);
1313       hitr->second.h8->SetBinError(i,std::sqrt(error));
1314     } else {
1315       hitr->second.h9->SetBinContent(i,mean);
1316       hitr->second.h9->SetBinError(i,std::sqrt(error));
1317     }
1318     double n     = (varmethod) ? hitr->second.h2->GetEntries() : hitr->second.h1->GetEntries();
1319     double vcont = (varmethod) ? hitr->second.h2->GetBinContent(i) : hitr->second.h1->GetBinContent(i);
1320     double econt = (varmethod) ? hitr->second.h2->GetBinError(i)   : hitr->second.h1->GetBinError(i);
1321     double vcfac = (varmethod) ? ((vcont > 0) ? (std::sqrt(mean/vcont)) : 1) : ((vcont > 0) ? (mean/vcont) : 1);
1322     double ecfac = (varmethod) ? ((vcont>0 && mean>0) ? (0.5*std::sqrt(((error)/(vcont*mean*n)) + (econt*econt)/(vcont*vcont) *((mean/vcont)-(2.0/n)))) : 0) : ((vcont>0 && mean>0) ? (std::sqrt(((error)/(n*vcont*vcont)) + (econt/vcont)*(econt/vcont)*((mean/vcont)*(mean/vcont)-(2.0/n)*(mean/vcont)))) : 0);
1323     corrfac.push_back(std::pair<double,double>(vcfac,ecfac));
1324   }
1325   return corrfac;
1326 }
1327
1328 void RecJet::StoreCorrFactor(const std::map<unsigned int,RecJet::Hists>& hh_,
1329                  std::map<unsigned int,RecJet::CFactors>& cfacs,
1330                  bool clear) {
1331
1332   if (clear) cfacs.clear();
1333   for (std::map<unsigned int,RecJet::Hists>::const_iterator hitr = hh_.begin();
1334        hitr != hh_.end(); ++hitr) {
1335     unsigned int detId2  = (hitr->first);
1336     if (hitr->second.h3 != 0) {
1337       for (int i=1; i<=hitr->second.h3->GetNbinsX(); ++i) {
1338     if (hitr->second.h3->GetBinContent(i) > 0) {
1339       unsigned int detId = (detId2 | i);
1340       std::map<unsigned int,RecJet::CFactors>::iterator itr = cfacs.find(detId);
1341       if (itr == cfacs.end()) {
1342         RecJet::CFactors fac;
1343         cfacs[detId] = fac;
1344         itr          = cfacs.find(detId);
1345       }
1346       itr->second.cfac1 = (hitr->second.h3->GetBinContent(i));
1347       itr->second.efac1 = (hitr->second.h3->GetBinError(i));
1348     }
1349       }
1350     }
1351     if (hitr->second.h4 != 0) {
1352       for (int i=1; i<=hitr->second.h4->GetNbinsX(); ++i) {
1353     if (hitr->second.h4->GetBinContent(i) > 0) {
1354       unsigned int detId = (detId2 | i);
1355       std::map<unsigned int,RecJet::CFactors>::iterator itr = cfacs.find(detId);
1356       if (itr == cfacs.end()) {
1357         RecJet::CFactors fac;
1358         cfacs[detId] = fac;
1359         itr          = cfacs.find(detId);
1360       }
1361       itr->second.cfac2 = (hitr->second.h4->GetBinContent(i));
1362       itr->second.efac2 = (hitr->second.h4->GetBinError(i));
1363     }
1364       }
1365     }
1366   }
1367 }
1368
1369 std::pair<double,double> RecJet::MeanDeviation(int sdet, std::map<unsigned int,RecJet::CFactors>& cfacs) {
1370   double dev1(0), dev2(0);
1371   int    kount(0);
1372   for (std::map<unsigned int,RecJet::CFactors>::iterator itr=cfacs.begin();
1373        itr != cfacs.end(); ++itr) {
1374     int sdet0 = ((itr->first)>>20)&0x7;
1375     if (sdet0 == sdet) {
1376       kount++;
1377       dev1 +=  (((itr->second.cfac1) > 1.0) ? (1.0-(1.0/(itr->second.cfac1))) : (1.0-(itr->second.cfac1)));
1378       dev2 +=  (((itr->second.cfac2) > 1.0) ? (1.0-(1.0/(itr->second.cfac2))) : (1.0-(itr->second.cfac2)));
1379     }
1380   }
1381   if (kount > 0) {
1382     dev1 /= kount;
1383     dev2 /= kount;
1384   }
1385   std::cout << kount << " dev " << dev1 << ":" << dev2 << std::endl;
1386   return std::pair<double,double>(dev1,dev2);
1387 }
1388
1389 void RecJet::WriteCorrFactor(std::string& outFile, std::string& inFile, bool var) {
1390
1391   // Read (if exists) correction factors from an old file
1392   std::map<unsigned int, std::pair<double,double> > corrFactor;
1393   ifstream myfile;
1394   myfile.open(inFile.c_str());
1395   if (!myfile.is_open()) {
1396     std::cout << "** ERROR: Can't open file '" << inFile << "' for i/p"
1397           << std::endl;
1398   } else {
1399     while(1) {
1400       unsigned int detId;
1401       double       cfac, efac;
1402       myfile >> detId >> cfac >> efac;
1403       if (!myfile.good()) break;
1404       std::map<unsigned int,std::pair<double,double> >::iterator itr = corrFactor.find(detId);
1405       if (itr == corrFactor.end()) {
1406     corrFactor[detId].first = cfac;
1407     corrFactor[detId].second= efac;
1408       } else {
1409     itr->second.first  = cfac;
1410     itr->second.second = efac;
1411       }
1412     }
1413     myfile.close();
1414   }
1415
1416   // See if it is the same as the ones calculated
1417   double diff(0);
1418   int    kount(0);
1419   const  unsigned int det(4);
1420   for (std::map<unsigned int,RecJet::CFactors>::iterator itr = corrFactor_.begin();
1421        itr != corrFactor_.end(); ++itr) {
1422     unsigned int sdet    = ((itr->first)>>20)&0x7;
1423     unsigned int detId   = ((det&0xf)<<28)|((sdet&0x7)<<25);
1424     detId               |= ((itr->first)&0xfffff);
1425     std::map<unsigned int,std::pair<double,double> >::iterator citr = corrFactor.find(detId);
1426     double cfac = (var) ? (itr->second.cfac2) : (itr->second.cfac1);
1427 //    double efac = (var) ? (itr->second.efac2) : (itr->second.efac1);
1428     if (citr != corrFactor.end() && cfac != 0) {
1429       kount++;
1430       //std::cout << cfac << "\t" << citr->second << "\t" << kount << std::endl;
1431       diff += fabs(cfac - citr->second.first);
1432     }
1433   }
1434   bool same = (kount > 0) ? (diff/kount < 0.0001) : false;
1435   std::cout << "Results from comparison: " << diff << " in " << kount << " --> "
1436         << same << std::endl;
1437
1438   // Now write it out
1439   ofstream myfile1;
1440   myfile1.open(outFile.c_str());
1441   if (!myfile1.is_open()) {
1442     std::cout << "** ERROR: Can't open '" << outFile << std::endl;
1443   } else {
1444     myfile1 << std::setprecision(4) << std::setw(10) << "detId"
1445         << std::setw(10) << "sub_det" << std::setw(10) << "ieta"
1446         << std::setw(10) << "iphi" << std::setw(10) << "depth"
1447         << std::setw(10) << "corrfact" << std::setw(15) << "error"
1448         << std::endl;
1449     for (std::map<unsigned int,RecJet::CFactors>::iterator itr = corrFactor_.begin();
1450      itr != corrFactor_.end(); ++itr) {
1451       double cfac = (var) ? (itr->second.cfac2) : (itr->second.cfac1);
1452       double efac = (var) ? (itr->second.efac2) : (itr->second.efac1);
1453       if (cfac != 0) {
1454     unsigned int sdet    = ((itr->first)>>20)&0x7;
1455     unsigned int detId   = ((det&0xf)<<28)|((sdet&0x7)<<25);
1456     detId               |= ((itr->first)&0xfffff);
1457     std::map<unsigned int,std::pair<double,double> >::iterator citr = corrFactor.find(detId);
1458     if (citr != corrFactor.end() && (!same)) {
1459       cfac *= (citr->second.first);
1460       efac *= (citr->second.second);
1461     }
1462     int keta  = ((detId&0x2000) ? ((detId>>7)&0x3f) : -((detId>>7)&0x3f));
1463         int depthx= ((detId>>14)&0x1f);
1464         int kphi  = (detId&0x7f);
1465         myfile1 << std::setw(10) << std::hex << detId
1466         << std::setw(10) << std::dec << sdet
1467         << std::setw(10) << keta
1468                 << std::setw(10) << kphi
1469         << std::setw(10) << depthx
1470         << std::setw(10) << cfac
1471         << std::setw(15) << efac
1472             << std::endl;
1473       }
1474     }
1475     myfile1.close();
1476   }
1477 }
1478
1479 void RecJet::Error_study(std::string fname, std::string det, int depth){
1480
1481   int detector = (det == "HB") ? 1 : ( (det == "HE") ? 2 : 4);
1482   char name[100], title[100];
1483
1484   sprintf(name,"%s_%d_error", det.c_str(), depth);
1485   sprintf(title,"%s_depth%d_error", det.c_str(), depth);
1486   TH1F* h_err = new TH1F(name,title,100, 0., 0.005);
1487   h_err->GetXaxis()->SetTitle("Error");  h_err->GetXaxis()->CenterTitle();
1488
1489   sprintf(name,"%s_%d_corr", det.c_str(), depth);
1490   sprintf(title,"%s_depth%d_corr", det.c_str(), depth);
1491   TH1F* h_corr = new TH1F(name,title,100, 0., 2);
1492   h_corr->GetXaxis()->SetTitle("Corr_Factor"); h_corr->GetXaxis()->CenterTitle();
1493
1494   sprintf (name, "%s_depth%d", det.c_str(), depth);
1495   TCanvas* c  = new TCanvas(name, name, 500, 500);
1496   c->SetRightMargin(0.10);
1497   c->SetTopMargin(0.10);
1498   gStyle->SetOptStat(111111);
1499   c->Divide(2,1);
1500   c->cd(1);
1501
1502   ifstream myfile;
1503   sprintf(name, "%s", fname.c_str());
1504   std::cout << name << std::endl;
1505   myfile.open(name);
1506   if (!myfile.is_open()) {
1507     std::cout << "** ERROR: Can't open file '" << std::endl;
1508   } else {
1509     while(1) {
1510       unsigned int detId;
1511       int subdet, ieta, iphi, idepth;
1512       double       cfac, efac;
1513       myfile >> std::hex >> detId >> std::dec >> subdet >> ieta >> iphi >> idepth >> cfac >> efac;
1514       //      std::cout << std::hex << detId << std::dec << subdet << ieta << iphi << idepth << cfac << efac << std::endl;
1515       if (!myfile.good()) break;
1516       if(subdet != detector || idepth != depth) continue;
1517       h_err->Fill(efac);
1518       h_corr->Fill(cfac);
1519     }
1520     myfile.close();
1521   }
1522   h_corr->Draw();
1523   c->cd(2);
1524   h_err->Draw();
1525 }
1526
1527 void RecJet::LoopMapTrig() {
1528   //   In a ROOT session, you can do:
1529   //      Root > .L RecJet.C
1530   //      Root > RecJet t
1531   //      Root > t.GetEntry(12); // Fill t data members with entry number 12
1532   //      Root > t.Show();       // Show values of entry 12
1533   //      Root > t.Show(16);     // Read and show values of entry 16
1534   //      Root > t.Loop();       // Loop on all entries
1535   //
1536   //     This is the loop skeleton where:
1537   //    jentry is the global entry number in the chain
1538   //    ientry is the entry number in the current Tree
1539   //  Note that the argument to GetEntry must be:
1540   //    jentry for TChain::GetEntry
1541   //    ientry for TTree::GetEntry and TBranch::GetEntry
1542   //
1543   //       To read only selected branches, Insert statements like:
1544   // METHOD1:
1545   //    fChain->SetBranchStatus("*",0);  // disable all branches
1546   //    fChain->SetBranchStatus("branchname",1);  // activate branchname
1547   // METHOD2: replace line
1548   //    fChain->GetEntry(jentry);       //read all branches
1549   //by  b_branchname->GetEntry(ientry); //read only this branch
1550
1551   mTrig_.clear();
1552   if (fChain != 0) {
1553     Long64_t nentries = fChain->GetEntriesFast();
1554     Long64_t nbytes = 0, nb = 0;
1555     for (Long64_t jentry=0; jentry<nentries;jentry++) {
1556       Long64_t ientry = LoadTree(jentry);
1557       if (ientry < 0) break;
1558       nb = fChain->GetEntry(jentry);   nbytes += nb;
1559       unsigned int detId = ((mysubd<<20) | ((depth&0x1f)<<14) |
1560                 ((ieta>0)?(0x2000|(ieta<<7)):((-ieta)<<7)) |
1561                 (iphi&0x7f));
1562       std::pair<unsigned int,int>key(detId,trigbit);
1563       std::map<std::pair<unsigned int,int>,RecJet::MyInfo>::iterator mitr = mTrig_.find(key);
1564       if (mitr == mTrig_.end()) {
1565     RecJet::MyInfo info;
1566     mTrig_[key] = info ;
1567     mitr = mTrig_.find(key) ;
1568       }
1569       mitr->second.Mom0 += mom0_MB;
1570       mitr->second.Mom1 += mom1_MB;
1571       mitr->second.Mom2 += mom2_MB;
1572       mitr->second.Mom3 += mom3_MB;
1573       mitr->second.Mom4 += mom4_MB;
1574     }
1575     loadTrig_ = true;
1576   }
1577 }
1578
1579 void RecJet::LoopTrig(int sdet, int eta, int phi, int dep, std::string indx) {
1580
1581   std::string det[4] = {"HB", "HE", "HO", "HF"};
1582   int subd = ((sdet>=1 && sdet<=4) ? sdet : 1);
1583
1584   if (!loadTrig_) LoopMapTrig();
1585   unsigned int detId = ((subd<<20) | ((dep&0x1f)<<14) |
1586             ((eta>0)?(0x2000|(eta<<7)):((-eta)<<7)) | (phi&0x7f));
1587   char fname[80],name[80],title[100];
1588   sprintf (fname, "dTrigbit%s_%s%d%d%d.root",indx.c_str(),det[subd-1].c_str(),eta,phi,dep);
1589   TFile f(fname, "recreate");
1590   sprintf(name, "mom0%s%d%d%d",det[subd-1].c_str(),eta,phi,dep);
1591   sprintf(title, "Entries for %s i#eta=%d i#phi=%d depth=%d",det[subd-1].c_str(),eta,phi,dep);
1592   TH1I* h1 = new TH1I(name, title, 128, 0, 128);
1593   sprintf(name, "edm0%s%d%d%d",det[subd-1].c_str(),eta,phi,dep);
1594   sprintf(title, "Error/Mean for %s i#eta=%d i#phi=%d depth=%d",det[subd-1].c_str(),eta,phi,dep);
1595   TH1D* h2 = new TH1D(name, title, 128, 0, 128);
1596   std::string bitn[113] = {"ZeroBias", "AlwaysTrue",
1597                "Mu16er_TauJet40erORCenJet72er", "Mu12_EG10",
1598                "Mu5_EG15", "DoubleTauJet40er",
1599                "IsoEG20er_TauJet20er_NotWdEta0",
1600                "Mu16er_TauJet20er", "QuadJetC84",
1601                "DoubleJetC72", "DoubleJetC120",
1602                "SingleJet52", "SingleJet68", "SingleJet92",
1603                "SingleJet128", "SingleJet176", "SingleJet200",
1604                "SingleJet36", "DoubleTauJet36er",
1605                "DoubleTauJet44er", "DoubleMu0",
1606                "SingleMu30er", "Mu3_JetC92_WdEtaPhi2",
1607                "Mu3_JetC16_WdEtaPhi2", "Mu3_JetC52_WdEtaPhi2",
1608                "EG25er_HTT125", "SingleEG35er", "SingleIsoEG25er",
1609                "SingleIsoEG25", "SingleIsoEG28er",
1610                "SingleIsoEG30er", "SingleEG10", "SingleIsoEG22er",
1611                "SingleJet240", "DoubleJetC52", "DoubleJetC84",
1612                "SingleMu14_Eta2p1", "DoubleJetC112",
1613                "DoubleMu_10_Open", "DoubleMu_10_3p5", "QuadJetC40",
1614                "SingleEG5", "SingleEG25", "SingleEG40",
1615                "SingleIsoEG18", "SingleIsoEG20er", "SingleEG20",
1616                "SingleEG30", "SingleEG35", "SingleMuOpen",
1617                "SingleMu16", "SingleMu5", "SingleMu20er",
1618                "SingleMu12", "SingleMu20", "SingleMu25er",
1619                "SingleMu25", "SingleMu30", "ETM30", "ETM50",
1620                "ETM70", "ETM100", "HTT125", "HTT150", "HTT175",
1621                "HTT200", "Mu20_EG10", "Mu5_EG20", "Mu5_IsoEG18",
1622                "Mu6_DoubleEG", "SingleJetC32_NotBptxOR",
1623                "ETM40", "HTT250", "Mu20_EG8", "Mu6_HTT150",
1624                "Mu10er_ETM50", "Mu14er_ETM30", "DoubleMu7_EG7",
1625                "SingleMu16_Eta2p1", "Mu8_HTT125", "Mu4_EG18",
1626                "SingleMu6_NotBptxOR", "DoubleMu6_EG6",
1627                "Mu5_DoubleEG5", "DoubleEG6_HTT150",
1628                "QuadJetC36_TauJet52", "TripleMu0", "TripleMu_5_5_3",
1629                "TripleEG_14_10_8", "DoubleEG_15_10",
1630                "DoubleEG_22_10", "DoubleEG_20_10_1LegIso",
1631                "Mu0er_ETM55", "DoubleJetC60_ETM60",
1632                "DoubleJetC32_WdPhi7_HTT125",
1633                "Jet32_DoubleMu_Open_10_MuMuNotWdPhi23_JetMuWdPhi1",
1634                "Jet32_MuOpen_EG10_MuEGNotWdPhi3_JetMuWdPhi1",
1635                "ETM60", "DoubleJetC100", "QuadJetC60",
1636                "SingleJetC20_NotBptxOR", "DoubleJetC56_ETM60",
1637                "DoubleMu0_Eta1p6_WdEta18", "ETM60_NotJet52WdPhi2",
1638                "ETM70_NotJet52WdPhi2", "TripleEG5",
1639                "TripleJet_92_76_64", "QuadMu0", "SingleMu18er",
1640                "DoubleMu0_Eta1p6_WdEta18_OS", "DoubleMu_12_5",
1641                "DoubleMu_10_0_WdEta18", "SingleMuBeamHalo"};
1642   int ibit[113] = {0,     1,   8,   9,  10,  11,  12,  13,  14,  15,  16,
1643            17,   18,  19,  20,  21,  22,  23,  24,  25,  26,  27,
1644            28,   29,  30,  31,  32,  33,  34,  35,  36,  37,  38,
1645            39,   40,  41,  42,  43,  44,  45,  46,  47,  48,  49,
1646            50,   51,  52,  53,  54,  55,  56,  58,  60,  61,  62,
1647            63,   64,  65,  66,  67,  68,  69,  70,  71,  72,  73,
1648            74,   75,  76,  77,  78,  79,  80,  81,  82,  83,  84,
1649            85,   86,  87,  88,  89,  90,  91,  92,  93,  97,  98,
1650            100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110,
1651            111, 112, 113, 114, 115, 116, 117, 119, 121, 122, 123,
1652            124, 126, 127};
1653   for (unsigned int k=0; k<113; ++k) {
1654     if (k%5 == 0) {
1655       h1->GetXaxis()->SetBinLabel(ibit[k]+1, bitn[k].c_str());
1656       h2->GetXaxis()->SetBinLabel(ibit[k]+1, bitn[k].c_str());
1657     }
1658   }
1659
1660   for (std::map<std::pair<unsigned int,int>,RecJet::MyInfo>::iterator mitr = mTrig_.begin();
1661        mitr !=mTrig_.end(); mitr++) {
1662     if (mitr->first.first == detId) {
1663       int bin     = mitr->first.second;
1664       double mean = mitr->second.Mom1/mitr->second.Mom0;
1665       double error= std::sqrt(((mitr->second.Mom2/mitr->second.Mom0)-mean*mean)/mitr->second.Mom0);
1666       std::pair<double,double> nsub = SubtractNoise(detId,mean,error,false);
1667       double edm = nsub.second/nsub.first;
1668       std::cout << "Before noise " << mean << " +- " << error << " after " << nsub.first << " +- " << nsub.second << " edm " << edm << std::endl;
1669       if (bin > 0 && bin <= 128) {
1670     h1->SetBinContent(bin+1,mitr->second.Mom0);
1671     h2->SetBinContent(bin+1,edm);
1672       }
1673     }
1674   }
1675   h1->Write();
1676   h2->Write();
1677   f.Close();
1678 }
1679
1680 std::map<unsigned int,RecJet::MyInfo> RecJet::BuildMap(int sdet, bool cfirst, double emin, double emax) {
1681
1682   std::map<unsigned int,RecJet::MyInfo> m_;
1683   int subd = ((sdet>=1 && sdet<=4) ? (sdet) : 0);
1684   std::string det[4] = {"HB", "HE", "HO", "HF"};
1685   int etamin[4] = {-16, -29, -15, -41};
1686   int etamax[4] = { 16,  29,  15,  41};
1687   int depmax[4] = {  2,   3,   4,   2};
1688   for (int eta=etamin[subd]; eta<=etamax[subd]; ++eta) {
1689     for (int phi=1; phi<=72; ++phi) {
1690       for (int depthx=1; depthx<=depmax[subd]; ++depthx) {
1691     char name[20];
1692     sprintf (name, "%seta%dphi%ddep%d", det[subd].c_str(), eta, phi, depthx);
1693     TH1D* hist = (TH1D*)file->FindObjectAny(name);
1694     if (hist != 0) {
1695       unsigned int detId = ((subd<<20) | ((depthx&0x1f)<<14) |
1696                 ((eta>0)?(0x2000|(eta<<7)):((-eta)<<7)) |
1697                 (phi&0x7f));
1698       double cfac = (cfirst) ? 1.0 : corrFactor_[detId].cfac1;
1699       std::map<unsigned int,RecJet::MyInfo>::iterator mitr = m_.find(detId);
1700       if (mitr == m_.end()) {
1701         RecJet::MyInfo info;
1702         m_[detId] = info ;
1703         mitr = m_.find(detId) ;
1704       }
1705       for (int i=1; i<=hist->GetNbinsX(); ++i) {
1706         double e = cfac*(hist->GetBinLowEdge(i)+0.5*hist->GetBinWidth(i));
1707         if (e > emin && e < emax) {
1708           double cont = hist->GetBinContent(i);
1709           mitr->second.Mom0 += cont;
1710           mitr->second.Mom1 += (cont*e);
1711           mitr->second.Mom2 += (cont*e*e);
1712           mitr->second.Mom3 += (cont*e*e*e);
1713           mitr->second.Mom4 += (cont*e*e*e*e);
1714         }
1715       }
1716     }
1717       }
1718     }
1719   }
1720   return m_;
1721 }