Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:12 

0001 //////////////////////////////////////////////////////////
0002 // L3 iterative procedure 
0003 // for IsoTrack calibration
0004 //
0005 // CalibTree class contains ROOT-tree
0006 // generated with IsoTrackCalibration plugin
0007 //
0008 // version 5.0   September 2016
0009 // modified detID for CMSSW_8_X !!!!!!!!!!!!!!!!!!!!!
0010 // M. Chadeeva
0011 //////////////////////////////////////////////////////////
0012 
0013 #include <TSystem.h>
0014 #include <TStyle.h>
0015 #include <TCanvas.h>
0016 #include <TROOT.h>
0017 #include <TChain.h>
0018 #include <TFile.h>
0019 #include <TTree.h>
0020 #include <TH1.h>
0021 #include <TGraph.h>
0022 #include <TGraphErrors.h>
0023 //#include <TProfile.h>
0024 #include <TLegend.h>
0025 #include <TString.h>
0026 #include <TF1.h>
0027 
0028 #include <vector>
0029 #include <string>
0030 #include <iostream>
0031 #include <fstream>
0032 #include <map>
0033 #include <utility>
0034 
0035 //**********************************************************
0036 // Constants
0037 //**********************************************************
0038 const unsigned int MAXNUM_SUBDET = 100;
0039 const int FIRST_IETA_TR = 15;
0040 const int FIRST_IETA_HE = 18;
0041 const int FIRST_IETA_FWD_1 = 25;
0042 const int FIRST_IETA_FWD_2 = 26;
0043 const unsigned int N_DEPTHS = 3;
0044 /*
0045 // old detID format for CMSSW_7_X 
0046 const unsigned int PHI_MASK       = 0x7F;
0047 const unsigned int ETA_OFFSET     = 7;
0048 const unsigned int ETA_MASK       = 0x3F;
0049 const unsigned int ZSIDE_MASK     = 0x2000;
0050 const unsigned int DEPTH_OFFSET   = 14;
0051 const unsigned int DEPTH_MASK     = 0x1F;
0052 const unsigned int DEPTH_SET      = 0x1C000;
0053 unsigned int MASK(0xFF80); //merge phi 
0054 */
0055 // new detID format for CMSSW_8_X
0056 const unsigned int PHI_MASK       = 0x3FF;
0057 const unsigned int ETA_OFFSET     = 10;
0058 const unsigned int ETA_MASK       = 0x1FF;
0059 const unsigned int ZSIDE_MASK     = 0x80000;
0060 const unsigned int DEPTH_OFFSET   = 20;
0061 const unsigned int DEPTH_MASK     = 0xF;
0062 const unsigned int DEPTH_SET      = 0xF00000;
0063 
0064 //--------------------------------------------
0065 // merging depths
0066 
0067 const unsigned int MERGE_PHI_AND_DEPTHS = 1;
0068 const unsigned int MASK(0xFFC00); //merge phi and depth
0069 /*
0070 const unsigned int MERGE_PHI_AND_DEPTHS = 0;
0071 const unsigned int MASK(0xFFFC00); //merge phi
0072 */
0073 //-------------------------------------------
0074 
0075 // individual ieta rings
0076 const unsigned int MASK2(0); // no second mask
0077 const int N_ETA_RINGS_PER_BIN = 1;
0078 /*
0079 // twin (even+odd) ieta rings
0080 const unsigned int MASK2(0x80);
0081 const int N_ETA_RINGS_PER_BIN = 2;
0082 
0083 // 4-fold ieta rings
0084 const unsigned int MASK2(0x180);
0085 const int N_ETA_RINGS_PER_BIN = 4;
0086 */
0087 const int MAX_ONESIDE_ETA_RINGS = 30;
0088 const int HALF_NUM_ETA_BINS =
0089   (MAX_ONESIDE_ETA_RINGS + 1*(N_ETA_RINGS_PER_BIN>1))/N_ETA_RINGS_PER_BIN;
0090 const int NUM_ETA_BINS = 2*HALF_NUM_ETA_BINS + 1;
0091 
0092 const int MIN_N_TRACKS_PER_CELL = 50;
0093 const int MIN_N_ENTRIES_FOR_FIT = 150;
0094 const double MAX_REL_UNC_FACTOR = 0.2;
0095 
0096 const bool APPLY_CORRECTION_FOR_PU = true;
0097 const bool SINGLE_REFERENCE_RESPONSE = false;
0098 /*
0099 //--- base correction for PU as of 2015 MC studies
0100 const char *l3prefix1 = "base";
0101 const double DELTA_CUT = 0.0; 
0102 const double LINEAR_COR_COEF[5] = { -0.375, -0.375, -0.375, -0.375, -0.375 };
0103 const double SQUARE_COR_COEF[5] = { -0.450, -0.450, -0.450, -0.450, -0.450 };
0104 //const double UPPER_LIMIT_DELTA_PU_COR = 2.0;
0105 */
0106 //--- optimized correction for PU as of 2016 MC studies
0107 // tuned to get MPV after correction
0108 // at the level of that of single pion response w/o PU and w/o correction
0109 const char *l3prefix1 = "opt";
0110 const double DELTA_CUT = 0.02; 
0111 const double LINEAR_COR_COEF[5] = { -0.35, -0.35, -0.35, -0.35, -0.45 };
0112 const double SQUARE_COR_COEF[5] = { -0.65, -0.65, -0.65, -0.30, -0.10 };
0113 //const double UPPER_LIMIT_DELTA_PU_COR = 1.5;
0114 
0115 const double UPPER_LIMIT_RESPONSE_BEFORE_COR = 3.0;
0116 const double FLEX_SEL_FIRST_CONST = 20.0;  // 20.(for exp); or 16*2
0117 const double FLEX_SEL_SECOND_CONST = 18.0; // 18.(for exp);
0118 
0119 const double MIN_RESPONSE_HIST = 0.0;
0120 const double MAX_RESPONSE_HIST = UPPER_LIMIT_RESPONSE_BEFORE_COR;
0121 const int NBIN_RESPONSE_HIST = 480;
0122 const int NBIN_RESPONSE_HIST_IND = 120;
0123 const double FIT_RMS_INTERVAL = 1.5;
0124 const double RESOLUTION_HCAL = 0.3;
0125 const double LOW_RESPONSE = 0.5;
0126 const double HIGH_RESPONSE = 1.5;
0127 
0128 //std::cout.precision(3);
0129 
0130 //**********************************************************
0131 // Header with CalibTree class definition
0132 //**********************************************************
0133 
0134 //////////////////////////////////////////////////////////
0135 // Header with CalibTree class
0136 // for L3 iterative procedure
0137 // implemented in L3_IsoTrackCalibration.C
0138 // 
0139 // version 2.0 January 2016
0140 // M. Chadeeva
0141 //////////////////////////////////////////////////////////
0142 
0143 #include <TSystem.h>
0144 #include <TStyle.h>
0145 #include <TCanvas.h>
0146 #include <TROOT.h>
0147 #include <TChain.h>
0148 #include <TFile.h>
0149 #include <TTree.h>
0150 #include <TH1.h>
0151 #include <TGraph.h>
0152 #include <TGraphErrors.h>
0153 #include <TProfile.h>
0154 #include <TLegend.h>
0155 #include <TString.h>
0156 #include <TF1.h>
0157 
0158 #include <vector>
0159 #include <string>
0160 #include <iostream>
0161 #include <fstream>
0162 #include <map>
0163 #include <utility>
0164 
0165 
0166 //**********************************************************
0167 // Class with TTree containing parameters of selected events
0168 //**********************************************************
0169 class CalibTree {
0170 public :
0171   TChain          *fChain;   //!pointer to the analyzed TTree
0172   //TChain          *inChain;   //!pointer to the analyzed TChain
0173   Int_t           fCurrent; //!current Tree number in a TChain
0174 
0175   // Declaration of leaf types
0176   Int_t           t_Run;
0177   Int_t           t_Event;
0178   Int_t           t_nVtx;
0179   Int_t           t_nTrk;
0180   Double_t        t_EventWeight;
0181   Double_t        t_p;
0182   Double_t        t_pt;
0183   Int_t           t_ieta;
0184   Double_t        t_phi;
0185   Double_t        t_eMipDR;
0186   Double_t        t_eHcal;
0187   Double_t        t_eHcal10;
0188   Double_t        t_eHcal30;
0189   Double_t        t_hmaxNearP;
0190   Bool_t          t_selectTk;
0191   Bool_t          t_qltyMissFlag;
0192   Bool_t          t_qltyPVFlag;
0193 /*
0194   Double_t        t_l1pt;
0195   Double_t        t_l1eta;
0196   Double_t        t_l1phi;
0197   Double_t        t_l3pt;
0198   Double_t        t_l3eta;
0199   Double_t        t_l3phi;
0200 */
0201   Double_t        t_mindR1;
0202   Double_t        t_mindR2;
0203   std::vector<unsigned int> *t_DetIds;
0204   //std::vector<unsigned int> *t_DetIds1;
0205   //std::vector<unsigned int> *t_DetIds3;
0206   std::vector<double>  *t_HitEnergies;
0207   //std::vector<double>  *t_HitEnergies1;
0208   //std::vector<double>  *t_HitEnergies3;
0209   
0210   // List of branches
0211   TBranch        *b_t_Run;   //!
0212   TBranch        *b_t_Event;   //!
0213   TBranch        *b_t_nVtx;
0214   TBranch        *b_t_nTrk;
0215   TBranch        *b_t_EventWeight;   //!
0216   TBranch        *b_t_p;   //!
0217   TBranch        *b_t_pt;   //!
0218   TBranch        *b_t_ieta;   //!
0219   TBranch        *b_t_phi;   //!
0220   TBranch        *b_t_eMipDR;   //!
0221   TBranch        *b_t_eHcal;   //!
0222   TBranch        *b_t_eHcal10;   //!
0223   TBranch        *b_t_eHcal30;   //!
0224   TBranch        *b_t_hmaxNearP;   //!
0225   TBranch        *b_t_selectTk;   //!
0226   TBranch        *b_t_qltyMissFlag;   //!
0227   TBranch        *b_t_qltyPVFlag;   //!
0228 /*
0229   TBranch        *b_t_l1pt;   //!
0230   TBranch        *b_t_l1eta;   //!
0231   TBranch        *b_t_l1phi;   //!
0232   TBranch        *b_t_l3pt;   //!
0233   TBranch        *b_t_l3eta;   //!
0234   TBranch        *b_t_l3phi;   //!
0235 */
0236   TBranch        *b_t_mindR1;   //!
0237   TBranch        *b_t_mindR2;   //!
0238   TBranch        *b_t_DetIds;   //!
0239   //TBranch        *b_t_DetIds1;   //!
0240   //TBranch        *b_t_DetIds3;   //!
0241   TBranch        *b_t_HitEnergies;   //!
0242   //TBranch        *b_t_HitEnergies1;   //!
0243   //TBranch        *b_t_HitEnergies3;   //!
0244 
0245   //--- constructor & destructor
0246   //CalibTree(TTree *tree=0);
0247   CalibTree(TChain *tree,
0248         double min_enrHcal, double min_pt,
0249         double lim_mipEcal, double lim_charIso,
0250         double min_trackMom, double max_trackMom);
0251   virtual ~CalibTree();
0252   
0253   //--- functions
0254   virtual Int_t      GetEntry(Long64_t entry);
0255   virtual Long64_t   LoadTree(Long64_t entry);
0256   //virtual void     Init(TTree *tree);
0257   virtual void       Init(TChain *tree);
0258   virtual Bool_t     Notify();
0259   virtual Int_t      firstLoop(unsigned, bool, unsigned);
0260   virtual Double_t   loopForIteration(unsigned, unsigned, unsigned);
0261   virtual Double_t   lastLoop(unsigned, unsigned, bool, unsigned);
0262   Bool_t             goodTrack(int);
0263   Bool_t             getFactorsFromFile(std::string, unsigned);
0264   unsigned int       saveFactorsInFile(std::string);
0265   Bool_t             openOutputRootFile(std::string);
0266 
0267   //--- variables for iterations
0268   Double_t referenceResponse;
0269   Double_t referenceResponseHB;
0270   Double_t referenceResponseTR;
0271   Double_t referenceResponseHE;
0272   Double_t maxZtestFromWeights;
0273   Double_t maxSys2StatRatio;
0274   int maxNumOfTracksForIeta;
0275   std::map<unsigned int, double> factors;
0276   std::map<unsigned int, double> uncFromWeights;
0277   std::map<unsigned int, double> uncFromDeviation;
0278   std::map<unsigned int, int> subDetector_trk;
0279   std::map<unsigned int, int> subDetector_final;
0280   std::map<unsigned int, int> nTrks;
0281   std::map<unsigned int, int> nSubdetInEvent;
0282   std::map<unsigned int, int> nPhiMergedInEvent;
0283   std::map<unsigned int, double> sumOfResponse;
0284   std::map<unsigned int, double> sumOfResponseSquared;
0285 
0286   //--- variables for selection
0287   double minEnrHcal;
0288   double minTrackPt;
0289   double minTrackMom;
0290   double maxTrackMom;
0291   double limMipEcal;
0292   double limCharIso;
0293   double constForFlexSel;
0294   
0295   //--- variables for plotting
0296   TFile *foutRootFile;
0297   TH1F* e2p_init;
0298   TH1F* e2pHB_init;
0299   TH1F* e2pTR_init;
0300   TH1F* e2pHE_init;
0301   TH1F* e2p_last;
0302   TH1F* e2pHB_last;
0303   TH1F* e2pTR_last;
0304   TH1F* e2pHE_last;
0305   
0306   TH1F* ieta_lefttail;
0307   TH1F* ieta_righttail;
0308   /*
0309   TProfile* deltaVSieta;
0310   TProfile* deltaNorm;
0311   TProfile* eHcalDeltaHB;
0312   TProfile* eHcalDeltaTR;
0313   TProfile* eHcalDeltaHE;
0314   TProfile* eHcalDeltaHEfwd;
0315   TProfile* frespDeltaHB;
0316   TProfile* frespDeltaTR;
0317   TProfile* frespDeltaHE;
0318   TProfile* frespDeltaHEfwd;
0319   */
0320 };
0321 
0322 
0323 //**********************************************************
0324 // Description of function to run iteration
0325 //**********************************************************
0326 unsigned int runIterations(const char *inFileDir = ".", 
0327                const char *inFileNamePrefix = "outputFromAnalyzer",
0328                const int firstInputFileEnum = 0,
0329                const int lastInputFileEnum = 1,
0330                const unsigned maxNumberOfIterations = 10,
0331                //const char *l3prefix0 = "cf",
0332                const double minHcalEnergy = 10.0,
0333                const double minPt = 7.0,
0334                const double limitForChargeIsolation = 2.0,
0335                const double minTrackMomentum = 40.0,
0336                const double maxTrackMomentum = 60.0,
0337                const double limitForMipInEcal = 1.0,
0338                const bool shiftResponse = 1,
0339                const unsigned int subSample = 2,
0340                const bool isCrosscheck = false,
0341                const char *inTxtFilePrefix = "test",
0342                const char *treeDirName = "IsoTrackCalibration", 
0343                const char *treeName = "CalibTree",
0344                unsigned int Debug = 0)
0345 {
0346   // Debug:  0-no debugging; 1-short debug; >1 - number of events to be shown in detail
0347   // subSample: extract factors from odd (0), even(1) or all(2) events
0348   // limitForChargeIsolation:  <0 - flex. sel.
0349   // and corr. for PU
0350   
0351   if ( isCrosscheck )
0352     std::cout << "Test with previously extracted factors..." << std::endl;
0353   else
0354     std::cout << "Extracting factors using L3 algorithm and isolated tracks..." << std::endl;
0355 
0356   char l3prefix0[10] = "ref1";
0357   if ( !SINGLE_REFERENCE_RESPONSE ) sprintf(l3prefix0,"ref3");
0358   char l3prefix2[20] = "_noCor";
0359   if ( APPLY_CORRECTION_FOR_PU ) sprintf(l3prefix2,"_%s%02d",
0360                      l3prefix1, int(100*DELTA_CUT)
0361                      );             
0362   char l3prefix3[10] = "_mean";
0363   if ( shiftResponse ) sprintf(l3prefix3,"_mpv");
0364   char l3prefix4[8] = "_3dep";
0365   if ( MERGE_PHI_AND_DEPTHS ) sprintf(l3prefix4,"_merged");
0366   char l3prefix[40];
0367   sprintf(l3prefix,"%s%s%s%s", l3prefix0, l3prefix2, l3prefix3, l3prefix4);
0368 
0369   char isoPrefix[10] = "hard";
0370   if ( limitForChargeIsolation < 0 ) sprintf(isoPrefix, "flex");
0371     
0372   char fnameInput[120];
0373   char fnameOutRoot[120];
0374   char fnameOutTxt[120] = "dummy";
0375   char fnameInTxt[120]  = "dummy";
0376   char tname[100];
0377   
0378   TGraph *g_converge1 = new TGraph(maxNumberOfIterations);
0379   TGraph *g_converge2 = new TGraph(maxNumberOfIterations);
0380   TGraph *g_converge3 = new TGraph(maxNumberOfIterations);
0381 
0382   sprintf(tname, "%s/%s", treeDirName, treeName );
0383   TChain tree(tname);
0384 
0385   //--- combine tree from several enumerated files with the same prefix
0386   //    or one file w/o number (firstInputFileEnum = lastInputFileEnum < 0 )
0387 
0388   for ( int ik = firstInputFileEnum; ik <= lastInputFileEnum; ik++ ) {
0389     if ( ik < 0 ) 
0390       sprintf(fnameInput, "%s/%s.root", inFileDir, inFileNamePrefix);
0391     else if (ik < 10 )
0392       sprintf(fnameInput, "%s/%s_%1d.root", inFileDir, inFileNamePrefix, ik);
0393     else if (ik < 100 )
0394       sprintf(fnameInput, "%s/%s_%2d.root", inFileDir, inFileNamePrefix, ik);
0395     else if (ik < 1000 )
0396       sprintf(fnameInput, "%s/%s_%3d.root", inFileDir, inFileNamePrefix, ik);
0397     else
0398       sprintf(fnameInput, "%s/%s_%4d.root", inFileDir, inFileNamePrefix, ik);
0399 
0400     if ( !gSystem->Which("./", fnameInput ) ) { // check file availability
0401       std::cout << "File " << fnameInput << " doesn't exist." << std::endl;
0402     }
0403     else {
0404       tree.Add(fnameInput);
0405       std::cout << "Add tree from " << fnameInput 
0406             << "   total number of entries (tracks): "
0407         << tree.GetEntries() << std::endl;
0408     }
0409   }
0410   if ( tree.GetEntries() == 0 ) {
0411     std:: cout << "Tree is empty." << std::endl;
0412     return -2;
0413   }
0414 
0415   //--- Initialize tree
0416   CalibTree t(&tree,
0417           minHcalEnergy, minPt,
0418           limitForMipInEcal, limitForChargeIsolation,
0419           minTrackMomentum, maxTrackMomentum);
0420   
0421   //--- Define files
0422   if ( isCrosscheck ) {
0423     sprintf(fnameInTxt, "%s_%s%02d-%02d-%02d_p%02d-%02d_pt%02d_eh%02d_ee%1d_step%1d.txt",
0424         inTxtFilePrefix,
0425         isoPrefix, int(t.limCharIso),
0426         int(abs(FLEX_SEL_FIRST_CONST)),
0427         int(abs(FLEX_SEL_SECOND_CONST)),
0428         int(minTrackMomentum), int(maxTrackMomentum), int(minPt),
0429         int(minHcalEnergy), int(limitForMipInEcal),
0430         N_ETA_RINGS_PER_BIN
0431         );
0432     sprintf(fnameOutRoot,
0433         "test_%1d_%s_by_%s_%s%02d-%02d-%02d_p%02d-%02d_pt%02d_eh%02d_ee%1d_step%1d.root",
0434         subSample,
0435         inFileNamePrefix,
0436         inTxtFilePrefix,
0437         isoPrefix, int(t.limCharIso),
0438         int(abs(FLEX_SEL_FIRST_CONST)),
0439         int(abs(FLEX_SEL_SECOND_CONST)),
0440         int(minTrackMomentum), int(maxTrackMomentum), int(minPt),
0441         int(minHcalEnergy), int(limitForMipInEcal),
0442         N_ETA_RINGS_PER_BIN
0443         );
0444   }
0445   else {    
0446     sprintf(fnameOutTxt,
0447         "%s_%1d_%s_i%02d_%s%02d-%02d-%02d_p%02d-%02d_pt%02d_eh%02d_ee%1d_step%1d.txt",
0448         l3prefix,
0449         subSample,
0450         inFileNamePrefix,
0451         maxNumberOfIterations,
0452         isoPrefix, int(t.limCharIso),
0453         int(abs(FLEX_SEL_FIRST_CONST)),
0454         int(abs(FLEX_SEL_SECOND_CONST)),
0455         int(minTrackMomentum), int(maxTrackMomentum), int(minPt),
0456         int(minHcalEnergy), int(limitForMipInEcal),
0457         N_ETA_RINGS_PER_BIN
0458         );
0459     sprintf(fnameOutRoot,
0460         "%s_%1d_%s_i%02d_%s%02d-%02d-%02d_p%02d-%02d_pt%02d_eh%02d_ee%1d_step%1d.root",
0461         l3prefix,
0462         subSample,
0463         inFileNamePrefix,
0464         maxNumberOfIterations,
0465         isoPrefix, int(t.limCharIso),
0466         int(abs(FLEX_SEL_FIRST_CONST)),
0467         int(abs(FLEX_SEL_SECOND_CONST)),
0468         int(minTrackMomentum), int(maxTrackMomentum), int(minPt),
0469         int(minHcalEnergy), int(limitForMipInEcal),
0470         N_ETA_RINGS_PER_BIN
0471         );
0472   }  
0473   if ( !t.openOutputRootFile(fnameOutRoot) ) {
0474     std::cout << "Problems with booking output file " << fnameOutRoot << std::endl;
0475     return -1;
0476   }
0477   std::cout << "Correction for PU: ";
0478   if ( APPLY_CORRECTION_FOR_PU ) {
0479     std::cout << " applied for delta > " << DELTA_CUT << std::endl;
0480     std::cout << " for HB: " << LINEAR_COR_COEF[0]
0481           << " ; " << SQUARE_COR_COEF[0]
0482           << std::endl;
0483     std::cout << " for TR: " << LINEAR_COR_COEF[1]
0484           << " ; " << SQUARE_COR_COEF[1]
0485           << std::endl;
0486     std::cout << " for HE(<=24): " << LINEAR_COR_COEF[2]
0487           << " ; " << SQUARE_COR_COEF[2]
0488           << std::endl;
0489     std::cout << " for ieta=25: " << LINEAR_COR_COEF[3]
0490           << " ; " << SQUARE_COR_COEF[3]
0491           << std::endl;
0492     std::cout << " for ieta=26: " << LINEAR_COR_COEF[4]
0493           << " ; " << SQUARE_COR_COEF[4]
0494           << std::endl;
0495   }
0496   else
0497     std::cout << " no " << std::endl;
0498     
0499   /*
0500   std::cout << "Constant coefficient from charge isolation: "
0501         << t.constForFlexSel << std::endl; 
0502   */
0503 
0504   unsigned int numOfSavedFactors(0);
0505   int nEventsWithGoodTrack(0);
0506   double MPVfromLastFit(0);
0507   
0508   if ( isCrosscheck ) {
0509     // open txt file and fill map with factors
0510     if ( t.getFactorsFromFile(fnameInTxt, Debug) ) {
0511       nEventsWithGoodTrack = t.firstLoop(subSample, false, Debug);
0512       std::cout << "Number of events with good track = "
0513         << nEventsWithGoodTrack << std::endl;
0514       MPVfromLastFit = t.lastLoop(subSample, maxNumberOfIterations, true, Debug);
0515       std::cout << "Finish testing " << t.factors.size() << " factors from file "
0516         << fnameInTxt << std::endl;
0517       std::cout << "MPV from fit after last iteration = "
0518         << MPVfromLastFit << std::endl;
0519       std::cout << "Test plots saved in " << fnameOutRoot << std::endl;
0520     }
0521     else {
0522       std::cout << "File " << fnameInTxt << " doesn't exist." << std::endl;
0523     }
0524   }
0525   else {
0526     //--- Prepare initial histograms and count good track
0527     nEventsWithGoodTrack = t.firstLoop(subSample, shiftResponse, Debug);
0528     std::cout << "Number of events with good track = "
0529           << nEventsWithGoodTrack << std::endl;
0530     //--- Iterate
0531     for ( unsigned int k = 0; k < maxNumberOfIterations; ++k ) {
0532       g_converge1->SetPoint( k, k+1, t.loopForIteration(subSample, k+1, Debug) );
0533       g_converge2->SetPoint( k, k+1, t.maxZtestFromWeights );
0534       g_converge3->SetPoint( k, k+1, t.maxSys2StatRatio );
0535     }
0536     //--- Finish
0537     MPVfromLastFit = t.lastLoop(subSample, maxNumberOfIterations, false, Debug);
0538     numOfSavedFactors = t.saveFactorsInFile(fnameOutTxt);
0539 
0540     sprintf(tname,"Mean deviation for subdetectors with Ntrack>%d",
0541         MIN_N_TRACKS_PER_CELL);
0542     g_converge1->SetTitle(tname);
0543     g_converge1->GetXaxis()->SetTitle("iteration");
0544     t.foutRootFile->WriteTObject(g_converge1, "g_cvgD");
0545     sprintf(tname,"Max abs(Z-test) for factors");
0546     g_converge2->SetTitle(tname);
0547     g_converge2->GetXaxis()->SetTitle("iteration");
0548     t.foutRootFile->WriteTObject(g_converge2, "g_cvgW");
0549     sprintf(tname,"Max ratio of syst. to stat. uncertainty");
0550     g_converge3->SetTitle(tname);
0551     g_converge3->GetXaxis()->SetTitle("iteration");
0552     t.foutRootFile->WriteTObject(g_converge3, "g_cvgR");
0553     
0554     std::cout << "Finish adjusting factors after "
0555           << maxNumberOfIterations << " iterations" << std::endl;
0556     std::cout << "MPV from fit after last iteration = "
0557           << MPVfromLastFit << std::endl;
0558     std::cout << "Table with " << numOfSavedFactors << " factors"
0559           << " with more than " << MIN_N_TRACKS_PER_CELL << " tracks/subdetector"
0560           << " (from " << t.factors.size() << " available)"
0561           << " is written in file " << fnameOutTxt << std::endl;
0562     std::cout << "Plots saved in " << fnameOutRoot << std::endl;
0563   }
0564 
0565   return numOfSavedFactors;
0566 }
0567 
0568 //**********************************************************
0569 // CalibTree constructor
0570 //**********************************************************
0571 //CalibTree::CalibTree(TTree *tree) : fChain(0) {
0572 CalibTree::CalibTree(TChain *tree,
0573              double min_enrHcal,
0574              double min_pt,
0575              double lim_mipEcal,
0576              double lim_charIso,
0577              double min_trackMom,
0578              double max_trackMom )
0579 { //: fChain(0) {
0580   // if parameter tree is not specified (or zero), connect the file
0581   // used to generate this class and read the Tree.
0582   if (tree == 0) {
0583     TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("output.root");
0584     if (!f || !f->IsOpen()) {
0585       f = new TFile("output.root");
0586     }
0587     TDirectory * dir = (TDirectory*)f->Get("IsoTrackCalibration");
0588     dir->GetObject("CalibTree",tree);
0589   }
0590   Init(tree);
0591 
0592   referenceResponse = 1;
0593   maxNumOfTracksForIeta = 0;
0594   // initialization of maps
0595   factors.clear();
0596   uncFromWeights.clear();
0597   uncFromDeviation.clear();
0598   subDetector_trk.clear();
0599   subDetector_final.clear();
0600   nTrks.clear();
0601   nSubdetInEvent.clear();
0602   nPhiMergedInEvent.clear();
0603   sumOfResponse.clear();
0604   sumOfResponseSquared.clear();
0605   
0606   // selection
0607   minEnrHcal = min_enrHcal;
0608   minTrackPt = min_pt;
0609   minTrackMom = min_trackMom;
0610   maxTrackMom = max_trackMom;
0611   limMipEcal = lim_mipEcal;
0612   limCharIso = abs(lim_charIso);
0613   if ( lim_charIso < 0 ) 
0614     constForFlexSel = log(FLEX_SEL_FIRST_CONST/limCharIso)/FLEX_SEL_SECOND_CONST;
0615   else constForFlexSel = 0;
0616 }
0617 
0618 //**********************************************************
0619 // CalibTree destructor
0620 //**********************************************************
0621 CalibTree::~CalibTree() {
0622 
0623     foutRootFile->cd();
0624     foutRootFile->Write();
0625     foutRootFile->Close();
0626 
0627   if (!fChain) return;
0628   delete fChain->GetCurrentFile();
0629 }
0630 
0631 //**********************************************************
0632 // Get entry function
0633 //**********************************************************
0634 Int_t CalibTree::GetEntry(Long64_t entry) {
0635   // Read contents of entry.
0636   if (!fChain) return 0;
0637   return fChain->GetEntry(entry);
0638 }
0639 
0640 //**********************************************************
0641 // Load tree function
0642 //**********************************************************
0643 Long64_t CalibTree::LoadTree(Long64_t entry) {
0644   // Set the environment to read one entry
0645   if (!fChain) return -5;
0646   Long64_t centry = fChain->LoadTree(entry);
0647   if (centry < 0) return centry;
0648   if (fChain->GetTreeNumber() != fCurrent) {
0649     fCurrent = fChain->GetTreeNumber();
0650     Notify();
0651   }
0652   return centry;
0653 }
0654 
0655 //**********************************************************
0656 // Initialisation of TTree
0657 //**********************************************************
0658 void CalibTree::Init(TChain *tree) {
0659   // Set object pointer
0660   t_DetIds = 0;
0661   //t_DetIds1 = 0;
0662   //t_DetIds3 = 0;
0663   t_HitEnergies = 0;
0664   //t_HitEnergies1 = 0;
0665   //t_HitEnergies3 = 0;
0666   // Set branch addresses and branch pointers
0667   if (!tree) return;
0668   fChain = tree;
0669   fCurrent = -1;
0670   fChain->SetMakeClass(1);
0671   
0672   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0673   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0674   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0675   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0676   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0677   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0678   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0679   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0680   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0681   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0682   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0683   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0684   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0685   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0686   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0687   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0688   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0689 /*
0690   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0691   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0692   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0693   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0694   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0695   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0696 */
0697   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0698   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0699   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0700   //fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0701   //fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0702   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0703   //fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0704   //fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0705   Notify();
0706 }
0707 
0708 //**********************************************************
0709 // Notification when opening new file
0710 //**********************************************************
0711 Bool_t CalibTree::Notify() {
0712   // The Notify() function is called when a new file is opened. This
0713   // can be either for a new TTree in a TChain or when when a new TTree
0714   // is started when using PROOF. It is normally not necessary to make changes
0715   // to the generated code, but the routine can be extended by the
0716   // user if needed. The return value is currently not used.
0717   
0718   return kTRUE;
0719 }
0720 //**********************************************************
0721 // Open file and book histograms
0722 //**********************************************************
0723 bool CalibTree::openOutputRootFile(std::string fname)
0724 {
0725   bool decision = false;
0726   
0727   foutRootFile = new TFile(fname.c_str(), "RECREATE");
0728   if ( foutRootFile != NULL ) decision = true;  
0729   foutRootFile->cd();
0730 
0731   return decision;
0732 }
0733 
0734 //**********************************************************
0735 // Initial loop over events in the tree
0736 //**********************************************************
0737 Int_t CalibTree::firstLoop(unsigned int subsample,
0738                bool shiftResp,
0739                unsigned int debug)
0740 {
0741   char name[100];
0742   unsigned int ndebug(0);
0743   double maxRespForGoodTrack(0);
0744   double minRespForGoodTrack(1000);
0745   int nRespOverHistLimit(0);
0746   int ntrk_ieta[NUM_ETA_BINS];
0747   for ( int j = 0; j < NUM_ETA_BINS; j++ ) {
0748     ntrk_ieta[j] = 0;
0749   }
0750   
0751   char scorr[80] = "correction for PU";
0752   char sxlabel[80] ="(E^{cor}_{hcal} + E_{ecal})/p_{track}"; 
0753   if ( !APPLY_CORRECTION_FOR_PU ) {
0754     sprintf(scorr,"no correction for PU");
0755     sprintf(sxlabel,"(E_{hcal} + E_{ecal})/p_{track}");
0756   }
0757   
0758   TF1* f1 = new TF1("f1","gaus", MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
0759 
0760   //--------- initialize histograms for response -----------------------------------------
0761   sprintf(name,"Initial HB+HE: %s", scorr);
0762   e2p_init = new TH1F("e2p_init", name,
0763               NBIN_RESPONSE_HIST, MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
0764   e2p_init->Sumw2();
0765   e2p_init->GetXaxis()->SetTitle(sxlabel);
0766   
0767   sprintf(name,"Initial HB: %s", scorr);
0768   e2pHB_init = new TH1F("e2pHB_init", name,
0769             NBIN_RESPONSE_HIST/2, MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
0770   e2pHB_init->Sumw2();
0771   e2pHB_init->GetXaxis()->SetTitle(sxlabel);
0772 
0773   sprintf(name,"Initial TR: %s", scorr);
0774   e2pTR_init = new TH1F("e2pTR_init", name,
0775             NBIN_RESPONSE_HIST/10, MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
0776   e2pTR_init->Sumw2();
0777   e2pTR_init->GetXaxis()->SetTitle(sxlabel);
0778 
0779   sprintf(name,"Initial HE: %s", scorr);
0780   e2pHE_init = new TH1F("e2pHE_init", name,
0781             NBIN_RESPONSE_HIST/2, MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
0782   e2pHE_init->Sumw2();
0783   e2pHE_init->GetXaxis()->SetTitle(sxlabel);
0784 
0785   sprintf(name,"Response < %3.1f", LOW_RESPONSE);
0786   ieta_lefttail = new TH1F("ieta_lefttail", name,
0787                2*MAX_ONESIDE_ETA_RINGS,
0788                -MAX_ONESIDE_ETA_RINGS, MAX_ONESIDE_ETA_RINGS); 
0789   ieta_lefttail->GetXaxis()->SetTitle("i#eta");
0790 
0791   sprintf(name,"Response > %3.1f", HIGH_RESPONSE);
0792   ieta_righttail = new TH1F("ieta_righttail", name,
0793                2*MAX_ONESIDE_ETA_RINGS,
0794                -MAX_ONESIDE_ETA_RINGS, MAX_ONESIDE_ETA_RINGS);              
0795   ieta_righttail->GetXaxis()->SetTitle("i#eta");
0796   
0797 //--- initialize chain ----------------------------------------
0798   if (fChain == 0) return 0;
0799   Long64_t nentries = fChain->GetEntriesFast();
0800   Long64_t nb = 0;
0801   
0802   int nSelectedEvents(0);
0803 
0804   if ( debug > 0 ) { 
0805     std::cout << "---------- First loop -------------------------- " << std::endl;
0806   }
0807 // ----------------------- loop over events -------------------------------------  
0808   for (Long64_t jentry=0; jentry<nentries; jentry++) {
0809     Long64_t ientry = LoadTree(jentry);
0810     if ( ientry < 0 || ndebug > debug ) break;   
0811     nb = fChain->GetEntry(jentry);   //nbytes += nb;
0812     
0813     if ( (jentry%2 == subsample) ) continue;   // only odd or even events
0814  
0815 // --------------- selection of good track --------------------    
0816 
0817     if ( !goodTrack(t_ieta) ) continue;
0818     
0819     if ( debug > 1 ) {
0820       ndebug++;
0821       std::cout << "***Entry (Track) Number : " << ientry << "(" << jentry << ")"
0822         << " p/eHCal/eMipDR/nDets : " << t_p << "/" << t_eHcal 
0823         << "/" << t_eMipDR << "/" << (*t_DetIds).size() 
0824         << std::endl;
0825     }
0826     
0827     double eTotal(0.0);
0828     double eTotalWithEcal(0.0);
0829       
0830     // ---- loop over active subdetectors in the event for total energy ---
0831     unsigned int nDets = (*t_DetIds).size();
0832     for (unsigned int idet = 0; idet < nDets; idet++) { 
0833       eTotal += (*t_HitEnergies)[idet];
0834     }
0835     eTotalWithEcal = eTotal + t_eMipDR;    
0836 
0837 // --- Correction for PU  --------
0838     double eTotalCor(eTotal);
0839     double eTotalWithEcalCor(eTotalWithEcal);
0840     //double e10(0.0);
0841     //double e30(0.0);
0842     double correctionForPU(1.0);
0843     double de2p(0.0);
0844     int abs_t_ieta = abs(t_ieta);
0845     
0846     if ( APPLY_CORRECTION_FOR_PU ) { 
0847       /*
0848       for (unsigned int idet1 = 0; idet1 < (*t_DetIds1).size(); idet1++) { 
0849     e10 += (*t_HitEnergies1)[idet1];
0850       }
0851       for (unsigned int idet3 = 0; idet3 < (*t_DetIds3).size(); idet3++) { 
0852     e30 += (*t_HitEnergies3)[idet3];
0853       }
0854       
0855       de2p = (e30 - e10)/t_p;
0856       */
0857       de2p = (t_eHcal30 - t_eHcal10)/t_p;
0858 
0859       if ( de2p > DELTA_CUT ) {
0860     int icor = int(abs_t_ieta >= FIRST_IETA_TR) + int(abs_t_ieta >= FIRST_IETA_HE)
0861       + int(abs_t_ieta >= FIRST_IETA_FWD_1) + int(abs_t_ieta >= FIRST_IETA_FWD_2);
0862     correctionForPU = (1 + LINEAR_COR_COEF[icor]*(t_eHcal/t_p)*de2p
0863                *(1 + SQUARE_COR_COEF[icor]*de2p));
0864       }
0865     }    
0866 
0867     // check for possibility to correct for PU
0868     if ( correctionForPU <= 0 || correctionForPU > 1 ) continue;
0869     nSelectedEvents++;
0870 
0871     eTotalCor = eTotal*correctionForPU;
0872     eTotalWithEcalCor = eTotalCor + t_eMipDR;
0873 
0874     double response = eTotalWithEcalCor/t_p;
0875     std::map<unsigned int, bool> sameSubdet;
0876     sameSubdet.clear();
0877     double resp2 = response*response;
0878 
0879     for (unsigned int idet = 0; idet < nDets; idet++) { 
0880       unsigned int detId = ( (*t_DetIds)[idet] & MASK ) | MASK2 ;
0881 
0882       if ( debug > 1 ) {
0883     unsigned int detId0 = ( (*t_DetIds)[idet] & MASK ) ;
0884     std::cout << "jentry/idet/detId :: ieta/z/depth ::: "
0885           << std::dec
0886           << jentry << " / "
0887           << ((*t_DetIds)[idet]) << " / "
0888           << detId0 << "(" << detId << ")" << " :: "
0889           << ((detId0>>ETA_OFFSET) & ETA_MASK)
0890           << "(" << ((detId>>ETA_OFFSET) & ETA_MASK) << ")" << " / "
0891           << ((detId0&ZSIDE_MASK) ? 1 : -1)
0892           << "(" << ((detId&ZSIDE_MASK) ? 1 : -1) << ")" << " / "
0893           << ((detId0>>DEPTH_OFFSET)&DEPTH_MASK)
0894           << "(" << ((detId>>DEPTH_OFFSET)&DEPTH_MASK) << ")"
0895           << std::endl;
0896       }
0897       if (nPhiMergedInEvent.find(detId) != nPhiMergedInEvent.end()) 
0898     nPhiMergedInEvent[detId]++;
0899       else 
0900     nPhiMergedInEvent.insert(std::pair<unsigned int,int>(detId, 1));
0901         
0902       if (nTrks.find(detId) != nTrks.end()) {
0903     if ( sameSubdet.find(detId) == sameSubdet.end() ) {
0904       nTrks[detId]++;
0905       nSubdetInEvent[detId] += nDets;
0906       sumOfResponse[detId] += response;
0907       sumOfResponseSquared[detId] += resp2;
0908       sameSubdet.insert(std::pair<unsigned int,bool>(detId, true));
0909     }
0910       }
0911       else {
0912     nTrks.insert(std::pair<unsigned int,int>(detId, 1));
0913     nSubdetInEvent.insert(std::pair<unsigned int,int>(detId, nDets));
0914     sumOfResponse.insert(std::pair<unsigned int,double>(detId,response));
0915     sumOfResponseSquared.insert(std::pair<unsigned int,double>(detId,resp2));
0916     sameSubdet.insert(std::pair<unsigned int,bool>(detId, true));
0917     subDetector_trk.insert(std::pair<unsigned int,
0918                    int>( detId,((*t_DetIds)[idet] &0xe000000) / 0x2000000 ));
0919       }
0920       
0921     }
0922 
0923 // --- Fill initial histograms ---------------------------      
0924     e2p_init->Fill(response ,1.0);
0925 
0926     if ( abs_t_ieta < FIRST_IETA_TR )
0927       e2pHB_init->Fill(response ,1.0);
0928     else if ( abs_t_ieta < FIRST_IETA_HE )
0929       e2pTR_init->Fill(response ,1.0);
0930     else
0931       e2pHE_init->Fill(response ,1.0);
0932 
0933     if ( debug > 1 ) {
0934       std::cout << "***Entry : " << ientry
0935         << " ***ieta/p/Ecal/nDet : "
0936         << t_ieta << "/" << t_p
0937         << "/" << t_eMipDR << "/" << (*t_DetIds).size() 
0938         << " ***Etot/E10/E30/Ecor/cPU : " << t_eHcal
0939         << "/" << t_eHcal10 << "/" << t_eHcal30
0940         << "/" << eTotalCor << "/" << correctionForPU
0941         << "(" << de2p << ")"
0942         << std::endl;
0943     }
0944     if ( response > maxRespForGoodTrack  )
0945       maxRespForGoodTrack = response;
0946     if ( response < minRespForGoodTrack )
0947       minRespForGoodTrack = response;
0948     if ( response > MAX_RESPONSE_HIST )
0949       nRespOverHistLimit++;
0950 
0951     if ( response < LOW_RESPONSE ) ieta_lefttail->Fill(t_ieta);
0952     if ( response > HIGH_RESPONSE ) ieta_righttail->Fill(t_ieta);
0953 
0954     int jj = HALF_NUM_ETA_BINS + int(t_ieta/N_ETA_RINGS_PER_BIN) + (t_ieta>0);
0955     ntrk_ieta[jj]++;
0956     
0957   } // ------------------- end of loop over events -------------------------------------
0958 
0959   for ( int j = 0; j < NUM_ETA_BINS; j++ ) {
0960     if ( maxNumOfTracksForIeta < ntrk_ieta[j] ) maxNumOfTracksForIeta = ntrk_ieta[j];
0961   }
0962 
0963   double jeta[N_DEPTHS][MAXNUM_SUBDET];
0964   double nTrk[N_DEPTHS][MAXNUM_SUBDET];
0965   double nSub[N_DEPTHS][MAXNUM_SUBDET];
0966   double nPhi[N_DEPTHS][MAXNUM_SUBDET];
0967   double rms[N_DEPTHS][MAXNUM_SUBDET];
0968   unsigned int kdep[N_DEPTHS];
0969   for ( unsigned ik = 0; ik < N_DEPTHS; ik++ ) { kdep[ik] = 0; }
0970 
0971   // fill number of tracks
0972   std::map <unsigned int,int>::iterator nTrksItr = nTrks.begin();
0973   for (nTrksItr = nTrks.begin(); nTrksItr != nTrks.end(); nTrksItr++ ) {
0974     unsigned int detId = nTrksItr->first;
0975     int depth= ((detId>>DEPTH_OFFSET) & DEPTH_MASK) + int(MERGE_PHI_AND_DEPTHS);
0976     int zside= (detId&ZSIDE_MASK) ? 1 : -1;
0977     unsigned int kcur = kdep[depth-1];
0978     
0979     jeta[depth-1][kcur] = int((detId>>ETA_OFFSET) & ETA_MASK)*zside;
0980     nTrk[depth-1][kcur] = nTrksItr->second;
0981     nSub[depth-1][kcur] = double(nSubdetInEvent[detId])/double(nTrksItr->second);
0982     nPhi[depth-1][kcur] = double(nPhiMergedInEvent[detId])/double(nTrksItr->second);
0983     if ( nTrk[depth-1][kcur] > 1 ) 
0984       rms[depth-1][kcur] = sqrt((sumOfResponseSquared[detId] -
0985                  pow(sumOfResponse[detId],2)/nTrk[depth-1][kcur])
0986                 /(nTrk[depth-1][kcur] - 1));
0987     else rms[depth-1][kcur] = RESOLUTION_HCAL;
0988     kdep[depth-1]++;
0989   }
0990   for ( unsigned ik = 0; ik < N_DEPTHS; ik++ ) {
0991     double x[MAXNUM_SUBDET];
0992     double ytrk[MAXNUM_SUBDET], ysub[MAXNUM_SUBDET];
0993     double yphi[MAXNUM_SUBDET], yrms[MAXNUM_SUBDET];
0994     for ( unsigned im = 0; im < MAXNUM_SUBDET; im++ ) {
0995       x[im] = jeta[ik][im];
0996       ytrk[im] = nTrk[ik][im];
0997       ysub[im] = nSub[ik][im];
0998       yphi[im] = nPhi[ik][im];
0999       yrms[im] = rms[ik][im];
1000     }
1001     TGraph*  g_ntrk = new TGraph(kdep[ik], x, ytrk);
1002     sprintf(name, "Number of tracks for depth %1d", ik+1);
1003     g_ntrk->SetTitle(name);
1004     sprintf(name, "nTrk_depth%1d", ik+1);
1005     foutRootFile->WriteTObject(g_ntrk, name);
1006 
1007     TGraph*  g_nsub = new TGraph(kdep[ik], x, ysub);
1008     sprintf(name, "Mean number of active subdetectors, depth %1d", ik+1);
1009     g_nsub->SetTitle(name);
1010     sprintf(name, "nSub_depth%1d", ik+1);
1011     foutRootFile->WriteTObject(g_nsub, name);
1012 
1013     TGraph*  g_nphi = new TGraph(kdep[ik], x, yphi);
1014     sprintf(name, "Mean number of phi-merged subdetectors, depth %1d", ik+1);
1015     g_nphi->SetTitle(name);
1016     sprintf(name, "nPhi_depth%1d", ik+1);
1017     foutRootFile->WriteTObject(g_nphi, name);
1018 
1019     TGraph*  g_rms = new TGraph(kdep[ik], x, yrms);
1020     sprintf(name, "RMS of samples, depth %1d", ik+1);
1021     g_rms->SetTitle(name);
1022     sprintf(name, "rms_depth%1d", ik+1);
1023     foutRootFile->WriteTObject(g_rms, name);
1024   }
1025 
1026   //--- estimate ratio mean/MPV
1027   double xl = e2p_init->GetMean() - FIT_RMS_INTERVAL*e2p_init->GetRMS();
1028   double xr = e2p_init->GetMean() + FIT_RMS_INTERVAL*e2p_init->GetRMS();
1029   e2p_init->Fit("f1","QN", "R", xl, xr);
1030   xl = f1->GetParameter(1) - FIT_RMS_INTERVAL*f1->GetParameter(2);
1031   xr = f1->GetParameter(1) + FIT_RMS_INTERVAL*f1->GetParameter(2);
1032   e2p_init->Fit("f1","QN", "R", xl, xr);
1033 
1034   if ( shiftResp && (f1->GetParameter(1) != 0) ) {
1035     referenceResponse = e2p_init->GetMean()/f1->GetParameter(1);
1036     std::cout << "Use reference response=<mean from sample>/<mpv from fit>:"
1037           << e2p_init->GetMean() << "/" << f1->GetParameter(1)
1038           << " = " << referenceResponse //<< std::endl
1039           << " (chi2ndf = " << f1->GetChisquare()/f1->GetNDF() << ")"
1040           << std::endl;
1041   }
1042   else {
1043     referenceResponse = 1;
1044     std::cout << "Use reference response = 1" << std::endl
1045           << "<mean from sample>/<mpv from fit> = "
1046           << e2p_init->GetMean()/f1->GetParameter(1)
1047           << "  (chi2ndf = " << f1->GetChisquare()/f1->GetNDF() << ")"
1048           << std::endl;
1049   }
1050   //---- for HB
1051   xl = e2pHB_init->GetMean() - FIT_RMS_INTERVAL*e2pHB_init->GetRMS();
1052   xr = e2pHB_init->GetMean() + FIT_RMS_INTERVAL*e2pHB_init->GetRMS();
1053   e2pHB_init->Fit("f1","QN", "R", xl, xr);
1054   xl = f1->GetParameter(1) - FIT_RMS_INTERVAL*f1->GetParameter(2);
1055   xr = f1->GetParameter(1) + FIT_RMS_INTERVAL*f1->GetParameter(2);
1056   e2pHB_init->Fit("f1","QN", "R", xl, xr);
1057 
1058   if ( shiftResp && (f1->GetParameter(1) != 0) ) {
1059     referenceResponseHB = e2pHB_init->GetMean()/f1->GetParameter(1);
1060     std::cout << "In HB <mean from sample>/<mpv from fit> = "
1061           << e2pHB_init->GetMean() << "/" << f1->GetParameter(1)
1062           << " = " << referenceResponseHB //<< std::endl
1063           << " (chi2ndf = " << f1->GetChisquare()/f1->GetNDF() << ")"
1064           << std::endl;
1065   }
1066   else {
1067     referenceResponseHB = 1;
1068     std::cout << "Use reference response in HB = 1" << std::endl
1069           << "<mean from sample>/<mpv from fit> = "
1070           << e2pHB_init->GetMean()/f1->GetParameter(1)
1071           << "  (chi2ndf = " << f1->GetChisquare()/f1->GetNDF() << ")"
1072           << std::endl;
1073   }
1074   //---- for TR
1075   xl = e2pTR_init->GetMean() - FIT_RMS_INTERVAL*e2pTR_init->GetRMS();
1076   xr = e2pTR_init->GetMean() + FIT_RMS_INTERVAL*e2pTR_init->GetRMS();
1077   e2pTR_init->Fit("f1","QN", "R", xl, xr);
1078   xl = f1->GetParameter(1) - FIT_RMS_INTERVAL*f1->GetParameter(2);
1079   xr = f1->GetParameter(1) + FIT_RMS_INTERVAL*f1->GetParameter(2);
1080   e2pTR_init->Fit("f1","QN", "R", xl, xr);
1081 
1082   if ( shiftResp && (f1->GetParameter(1) != 0) ) {
1083     referenceResponseTR = e2pTR_init->GetMean()/f1->GetParameter(1);
1084     std::cout << "In TR <mean from sample>/<mpv from fit> = "
1085           << e2pTR_init->GetMean() << "/" << f1->GetParameter(1)
1086           << " = " << referenceResponseTR //<< std::endl
1087           << " (chi2ndf = " << f1->GetChisquare()/f1->GetNDF() << ")"
1088           << std::endl;
1089   }
1090   else {
1091     referenceResponseTR = 1;
1092     std::cout << "Use reference response in TR = 1" << std::endl
1093           << "<mean from sample>/<mpv from fit> = "
1094           << e2pTR_init->GetMean()/f1->GetParameter(1)
1095           << "  (chi2ndf = " << f1->GetChisquare()/f1->GetNDF() << ")"
1096           << std::endl;
1097   }
1098   //---- for HE
1099   xl = e2pHE_init->GetMean() - FIT_RMS_INTERVAL*e2pHE_init->GetRMS();
1100   xr = e2pHE_init->GetMean() + FIT_RMS_INTERVAL*e2pHE_init->GetRMS();
1101   e2pHE_init->Fit("f1","QN", "R", xl, xr);
1102   xl = f1->GetParameter(1) - FIT_RMS_INTERVAL*f1->GetParameter(2);
1103   xr = f1->GetParameter(1) + FIT_RMS_INTERVAL*f1->GetParameter(2);
1104   e2pHE_init->Fit("f1","QN", "R", xl, xr);
1105 
1106   if ( shiftResp && (f1->GetParameter(1) != 0) ) {
1107     referenceResponseHE = e2pHE_init->GetMean()/f1->GetParameter(1);
1108     std::cout << "In HE <mean from sample>/<mpv from fit> = "
1109           << e2pHE_init->GetMean() << "/" << f1->GetParameter(1)
1110           << " = " << referenceResponseHE //<< std::endl
1111           << " (chi2ndf = " << f1->GetChisquare()/f1->GetNDF() << ")"
1112           << std::endl;
1113   }
1114   else {
1115     referenceResponseHE = 1;
1116     std::cout << "Use reference response in HE = 1" << std::endl
1117           << "<mean from sample>/<mpv from fit> = "
1118           << e2pHE_init->GetMean()/f1->GetParameter(1)
1119           << "  (chi2ndf = " << f1->GetChisquare()/f1->GetNDF() << ")"
1120           << std::endl;
1121   }
1122 
1123   //----- print additional info
1124   std::cout << "Maximal response for good tracks = " 
1125         << maxRespForGoodTrack << std::endl
1126         << nRespOverHistLimit
1127         << " events with response > " << MAX_RESPONSE_HIST
1128         << "(hist limit for mean estimate)"
1129         << std::endl;
1130   std::cout << "Minimal response for good tracks = " 
1131         << minRespForGoodTrack
1132         << std::endl;
1133   std::cout << "Maximum number of selected tracks per ieta bin = " 
1134         << maxNumOfTracksForIeta
1135         << std::endl;
1136   std::cout << "Number of selected tracks in HB = "
1137         << e2pHB_init->GetEntries()
1138         << std::endl;
1139   std::cout << "Number of selected tracks in TR = "
1140         << e2pTR_init->GetEntries()
1141         << std::endl;
1142   std::cout << "Number of selected tracks in HE = "
1143         << e2pHE_init->GetEntries()
1144         << std::endl;
1145 
1146   /*
1147   xl = e2pHB_init->GetMean() - FIT_RMS_INTERVAL*e2pHB_init->GetRMS();
1148   xr = e2pHB_init->GetMean() + FIT_RMS_INTERVAL*e2pHB_init->GetRMS();
1149   e2pHB_init->Fit("f1","QN", "R", xl, xr);
1150 
1151   xl = e2pHE_init->GetMean() - FIT_RMS_INTERVAL*e2pHE_init->GetRMS();
1152   xr = e2pHE_init->GetMean() + FIT_RMS_INTERVAL*e2pHE_init->GetRMS();
1153   e2pHE_init->Fit("f1","QN", "R", xl, xr);
1154   */
1155   return nSelectedEvents;
1156 }
1157 
1158 //**********************************************************
1159 // Loop over events in the tree for current iteration
1160 //**********************************************************
1161 Double_t CalibTree::loopForIteration(unsigned int subsample,
1162                      unsigned int nIter,
1163                      unsigned int debug )
1164 {
1165   char name[500];
1166   double meanDeviation = 0;
1167   unsigned int ndebug(0);
1168     
1169   TF1* f1 = new TF1("f1","gaus", MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
1170   TH1F* e2p[NUM_ETA_BINS];
1171 
1172   int n_ieta_bins = 2*2.5*pow(maxNumOfTracksForIeta,1/3.0);
1173   for ( int i = 0; i < NUM_ETA_BINS; i++ ) {
1174     sprintf(name,"e2p[%02d]", i);
1175     e2p[i] = new TH1F(name, "",
1176               n_ieta_bins, //NBIN_RESPONSE_HIST_IND,
1177               MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
1178     e2p[i]->Sumw2();
1179   }
1180 
1181   std::map<unsigned int, std::pair<double,double> > sumsForFactorCorrection;
1182   std::map<unsigned int, double> sumOfWeightsSquared;
1183 
1184   if ( debug > 0 ) {
1185     std::cout.precision(3);
1186     std::cout << "-------------------------------------------- nIter = "
1187           << nIter << std::endl;
1188   }
1189 //--- initialize chain ----------------------------------------
1190   if (fChain == 0) return 0;
1191   Long64_t nentries = fChain->GetEntriesFast();
1192   Long64_t nb = 0;
1193   
1194 // ----------------------- loop over events -------------------------------------  
1195   for (Long64_t jentry=0; jentry<nentries; jentry++) {
1196     Long64_t ientry = LoadTree(jentry);
1197     if ( ientry < 0 || ndebug > debug ) break;   
1198     nb = fChain->GetEntry(jentry);   //nbytes += nb;
1199     
1200     if ( (jentry%2 == subsample) ) continue;   // only odd or even events
1201 
1202     // --------------- selection of good track --------------------
1203     
1204     if ( !goodTrack(t_ieta) ) continue;
1205 
1206     if ( debug > 1 ) {
1207       ndebug++;
1208       std::cout << "***Entry (Track) Number : " << ientry 
1209         << " p/eHCal/eMipDR/nDets : " << t_p << "/" << t_eHcal 
1210         << "/" << t_eMipDR << "/" << (*t_DetIds).size() 
1211         << std::endl;
1212     }
1213     
1214     double eTotal(0.0);
1215     double eTotalWithEcal(0.0);
1216       
1217     // ---- first loop over active subdetectors in the event for total energy ---
1218 
1219     for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) { 
1220       double hitEnergy(0);  
1221       unsigned int detId = ( (*t_DetIds)[idet] & MASK ) | MASK2 ;
1222     
1223       if (factors.find(detId) != factors.end()) 
1224     hitEnergy = factors[detId] * (*t_HitEnergies)[idet];
1225       else 
1226     hitEnergy = (*t_HitEnergies)[idet];
1227 
1228       eTotal += hitEnergy;
1229     }
1230 
1231     eTotalWithEcal = eTotal + t_eMipDR;    
1232 
1233 // --- Correction for PU   --------      
1234     double eTotalCor(eTotal);
1235     double eTotalWithEcalCor(eTotalWithEcal);
1236     //double e10(0.0);
1237     //double e30(0.0);
1238     double correctionForPU(1.0);
1239 
1240     if ( APPLY_CORRECTION_FOR_PU ) {
1241       /*
1242       for (unsigned int idet1 = 0; idet1 < (*t_DetIds1).size(); idet1++) { 
1243     double hitEnergy(0);
1244     unsigned int detId1 = ( (*t_DetIds1)[idet1] & MASK ) | MASK2;
1245     
1246     if (factors.find(detId1) != factors.end()) 
1247       hitEnergy = factors[detId1] * (*t_HitEnergies1)[idet1];
1248     else 
1249       hitEnergy = (*t_HitEnergies1)[idet1];
1250     
1251     e10 += hitEnergy;
1252       }
1253       for (unsigned int idet3 = 0; idet3 < (*t_DetIds3).size(); idet3++) { 
1254     double hitEnergy(0);
1255     unsigned int detId3 = ( (*t_DetIds3)[idet3] & MASK ) | MASK2;
1256     
1257     if (factors.find(detId3) != factors.end()) 
1258       hitEnergy = factors[detId3] * (*t_HitEnergies3)[idet3];
1259     else
1260       hitEnergy = (*t_HitEnergies3)[idet3];
1261     
1262     e30 += hitEnergy;
1263       }
1264       double de2p = (e30 - e10)/t_p;
1265       */
1266       double de2p = (t_eHcal30 - t_eHcal10)/t_p;
1267       if ( de2p > DELTA_CUT ) {
1268     int abs_t_ieta = abs(t_ieta);
1269     int icor = int(abs_t_ieta >= FIRST_IETA_TR) + int(abs_t_ieta >= FIRST_IETA_HE)
1270       + int(abs_t_ieta >= FIRST_IETA_FWD_1) + int(abs_t_ieta >= FIRST_IETA_FWD_2);
1271     correctionForPU = (1 + LINEAR_COR_COEF[icor]*(t_eHcal/t_p)*de2p
1272                *(1 + SQUARE_COR_COEF[icor]*de2p));
1273       }
1274     }    
1275 
1276     // check for possibility to correct for PU
1277     if ( correctionForPU <= 0 || correctionForPU > 1 ) continue;
1278 
1279     eTotalCor = eTotal*correctionForPU;
1280     eTotalWithEcalCor = eTotalCor + t_eMipDR;
1281      
1282     int jeta = HALF_NUM_ETA_BINS + int(t_ieta/N_ETA_RINGS_PER_BIN) + (t_ieta>0);
1283     e2p[jeta]->Fill(eTotalWithEcalCor/t_p ,1.0);
1284       
1285 // ---- second loop over active subdetectors in the event  -----------------
1286       
1287     double response = eTotalWithEcalCor/t_p; // - referenceResponse;
1288     
1289     for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) {
1290       double hitEnergy(0);
1291       unsigned int detId = ( (*t_DetIds)[idet] & MASK ) | MASK2 ;
1292          
1293       if (factors.find(detId) != factors.end())
1294     hitEnergy = factors[detId] * (*t_HitEnergies)[idet];
1295       else
1296     hitEnergy = (*t_HitEnergies)[idet];
1297 
1298       double cellWeight = hitEnergy/eTotal;   
1299       //double trackWeight = (cellWeight * t_p) / eTotalWithEcalCor; // old method
1300       double trackWeight = cellWeight*response;   // new method
1301       double cellweight2 = cellWeight*cellWeight;
1302       
1303       if( sumsForFactorCorrection.find(detId) != sumsForFactorCorrection.end() ) {
1304     cellWeight  += sumsForFactorCorrection[detId].first;
1305     trackWeight += sumsForFactorCorrection[detId].second;
1306     sumsForFactorCorrection[detId] = std::pair<double,double>(cellWeight,trackWeight);
1307     sumOfWeightsSquared[detId] += cellweight2;
1308       }
1309       else {
1310     sumsForFactorCorrection.insert(std::pair<unsigned int,
1311                        std::pair<double,double> >(detId,
1312                                   std::pair<double,double>(cellWeight,
1313                                                trackWeight)));
1314     sumOfWeightsSquared.insert(std::pair<unsigned int,double>(detId, cellweight2));
1315      }
1316     
1317       if ( debug > 1 ) { //|| hitEnergy < -0.5) {
1318     double f = 1;
1319     int zside= (detId&ZSIDE_MASK) ? 1 : -1;
1320     if (factors.find(detId) != factors.end()) f = factors[detId];
1321     std::cout << jentry << "::: "
1322           << " Ncells: " << (*t_DetIds).size()
1323           << " !! detId(ieta)/e/f : " 
1324         //    << std::hex << (*t_DetIds)[idet] << ":"
1325           << detId << "(" << int((detId>>ETA_OFFSET) & ETA_MASK)*zside << ")"
1326           << "/" << hitEnergy
1327           << "/" << f
1328           << " ||| cellW/trW : " << cellWeight << " / " << trackWeight
1329           << " ||| E/Ecor/p : " << eTotal
1330           << " / " << eTotalCor
1331           << " / " << t_p
1332           << " || e10/e30/cF : " << t_eHcal10
1333           << " / " << t_eHcal30
1334           << " / " << correctionForPU
1335           << std::endl;
1336       }
1337     }  // --------------- end of second loop over cells ----------
1338   } // ------------------- end of loop over events -------------------------------------
1339 
1340 //----- Graphs to be saved in root file ----------------
1341   if ( debug > 0 ) {
1342     std::cout << "Fit and calculate means..." << std::endl;
1343     std::cout << "Number of plots (ieta bins) = " << NUM_ETA_BINS << std::endl;
1344   }
1345   TGraph *g_chi = new TGraph(NUM_ETA_BINS);  
1346   TGraphErrors* g_e2pFit = new TGraphErrors(NUM_ETA_BINS);
1347   TGraphErrors* g_e2pMean = new TGraphErrors(NUM_ETA_BINS);
1348   TGraph *g_nhistentries = new TGraph(NUM_ETA_BINS);
1349   
1350   int ipoint(0);
1351   for ( int i = 0; i < NUM_ETA_BINS; i++ ) {
1352     int ieta = (i - HALF_NUM_ETA_BINS - (i>HALF_NUM_ETA_BINS))*N_ETA_RINGS_PER_BIN;
1353     if ( N_ETA_RINGS_PER_BIN > 1 ) {
1354       ieta = (i > HALF_NUM_ETA_BINS) ? ieta+1 : ieta-1;
1355     }
1356     int nhistentries = e2p[i]->GetEntries();
1357     /*
1358       if ( debug > 0 ) {
1359     std::cout << "i / entries / ieta :::"
1360           << i
1361           << " / " << nhistentries
1362           << " / " << ieta
1363           << std::endl;
1364       }
1365     */
1366      if ( nIter == 1 ) {
1367        g_nhistentries->SetPoint(i, ieta, nhistentries);
1368      }
1369 
1370     if ( nhistentries < 1 ) continue;
1371     else {
1372       g_e2pMean->SetPoint(ipoint, ieta, e2p[i]->GetMean());
1373       g_e2pMean->SetPointError(ipoint, 0, e2p[i]->GetMeanError());
1374 
1375       if ( nhistentries > MIN_N_ENTRIES_FOR_FIT ) {
1376     int nrebin = n_ieta_bins/(2*2.5*pow(nhistentries,1/3.0));
1377     if ( nrebin > 2 ) e2p[i]->Rebin(nrebin);
1378     
1379     double xl = e2p[i]->GetMean() - FIT_RMS_INTERVAL*e2p[i]->GetRMS();
1380     double xr = e2p[i]->GetMean() + FIT_RMS_INTERVAL*e2p[i]->GetRMS();
1381     e2p[i]->Fit("f1","QN", "R", xl, xr);
1382     xl = f1->GetParameter(1) - FIT_RMS_INTERVAL*f1->GetParameter(2);
1383     xr = f1->GetParameter(1) + FIT_RMS_INTERVAL*f1->GetParameter(2);
1384     e2p[i]->Fit("f1","QN", "R", xl, xr);
1385     g_e2pFit->SetPoint(ipoint, ieta, f1->GetParameter(1));
1386     g_e2pFit->SetPointError(ipoint, 0, f1->GetParError(1));
1387     g_chi->SetPoint(ipoint, ieta, f1->GetChisquare()/f1->GetNDF());
1388       }
1389       else {
1390     g_e2pFit->SetPoint(ipoint, ieta, e2p[i]->GetMean());
1391     g_e2pFit->SetPointError(ipoint, 0, e2p[i]->GetMeanError());
1392     g_chi->SetPoint(ipoint, ieta, 0);
1393       }
1394       ipoint++;
1395     }
1396   }
1397   // fill number of tracks per ieta
1398   if ( nIter == 1 ) {
1399       sprintf(name, "Number of selected tracks");
1400       g_nhistentries->SetTitle(name);
1401       g_nhistentries->GetXaxis()->SetTitle("i#eta");
1402       sprintf(name, "selTrks");
1403       foutRootFile->WriteTObject(g_nhistentries, name);   
1404   }
1405 
1406   for ( int k = ipoint; k < NUM_ETA_BINS; k++ ) {
1407     g_e2pFit->RemovePoint(ipoint);
1408     g_e2pMean->RemovePoint(ipoint);
1409   }
1410   sprintf(name, "Response from fit, iteration %2d", nIter);
1411   g_e2pFit->SetTitle(name);
1412   g_e2pFit->GetXaxis()->SetTitle("i#eta");
1413   sprintf(name, "respFit_%d", nIter);
1414   foutRootFile->WriteTObject(g_e2pFit, name);
1415 
1416   sprintf(name, "Mean response, iteration %2d", nIter);
1417   g_e2pMean->SetTitle(name);
1418   g_e2pMean->GetXaxis()->SetTitle("i#eta");
1419   sprintf(name, "respMean_%d", nIter);
1420   foutRootFile->WriteTObject(g_e2pMean, name);
1421 
1422   sprintf(name, "Chi2/NDF, iteration %2d", nIter);
1423   g_chi->SetTitle(name);
1424   g_chi->GetXaxis()->SetTitle("i#eta");
1425   sprintf(name, "chi2ndf_%d", nIter);
1426   foutRootFile->WriteTObject(g_chi, name);
1427 
1428 // --- convergence criteria and correction factors -----------------------------------
1429 
1430   double MeanConvergenceDelta(0),  MaxRelDeviationWeights(0), MaxRatioUncertainties(0);
1431   double dets[MAXNUM_SUBDET];
1432   double ztest[MAXNUM_SUBDET], sys2statRatio[MAXNUM_SUBDET];
1433 
1434   if ( debug > 0 ) std::cout << "Calculate correction factors..." << std::endl;
1435 
1436   unsigned int kount(0), mkount(0);
1437   unsigned int maxKountW(0), maxKountR(0);
1438 
1439   double fac[N_DEPTHS][MAXNUM_SUBDET];
1440   double dfac[N_DEPTHS][MAXNUM_SUBDET];
1441   double ieta[N_DEPTHS][MAXNUM_SUBDET];
1442   double dieta[N_DEPTHS][MAXNUM_SUBDET];
1443 
1444   unsigned int kdep[N_DEPTHS];
1445   for ( unsigned ik = 0; ik < N_DEPTHS; ik++ ) { kdep[ik] = 0; }
1446 
1447 //-------------- loop over all cells ---------------------------------
1448   std::map <unsigned int,
1449         std::pair<double,double> >::iterator sumsForFactorCorrectionItr
1450     = sumsForFactorCorrection.begin();
1451   for (; sumsForFactorCorrectionItr != sumsForFactorCorrection.end();
1452        sumsForFactorCorrectionItr++) {
1453 
1454     unsigned int detId = sumsForFactorCorrectionItr->first;
1455     double sumOfWeights = (sumsForFactorCorrectionItr->second).first;
1456     int nSubDetTracks(0);
1457     double subdetRMS(RESOLUTION_HCAL);
1458     if ( nTrks.find(detId) != nTrks.end() ) {
1459       nSubDetTracks = nTrks[detId];
1460       if ( nSubDetTracks > 1 ) 
1461     subdetRMS = sqrt((sumOfResponseSquared[detId] 
1462               - pow(sumOfResponse[detId],2)/double(nSubDetTracks))
1463              /double(nSubDetTracks - 1));
1464     }
1465     else {
1466       std::cout << "!!!!!!! No tracks for subdetector " << detId << std::endl;
1467       continue;
1468     }
1469     double NcellMean = double(nSubdetInEvent[detId])/double(nSubDetTracks);
1470 
1471     double ratioWeights(1);
1472     if ( abs(sumOfWeights) > 0 )
1473       ratioWeights = sqrt(sumOfWeightsSquared[detId])/sumOfWeights;
1474     double correctionRMS = subdetRMS*ratioWeights*sqrt(NcellMean);
1475     
1476     double absErrorW(0);
1477     double absErrorWprevious(0);
1478     double factorPrevious(1);
1479     double factorCorrection(1);
1480 
1481     int zside = (detId&ZSIDE_MASK) ? 1 : -1;
1482     int depth = ((detId>>DEPTH_OFFSET)&DEPTH_MASK) + int(MERGE_PHI_AND_DEPTHS);
1483     unsigned int kcur = kdep[depth-1];
1484     
1485     ieta[depth-1][kcur] = int((detId>>ETA_OFFSET) & ETA_MASK)*zside;
1486     dieta[depth-1][kcur] = 0;
1487     
1488     double refR = referenceResponse;
1489     if ( !SINGLE_REFERENCE_RESPONSE ) {
1490       if ( abs(ieta[depth-1][kcur]) < FIRST_IETA_TR )
1491     refR = referenceResponseHB;
1492       else if ( abs(ieta[depth-1][kcur]) < FIRST_IETA_HE )
1493     refR = referenceResponseTR;
1494       else
1495     refR = referenceResponseHE;
1496     }
1497     
1498     //--- old expression ---------------
1499     /*
1500     factorCorrection = (sumsForFactorCorrectionItr->second).second
1501                      / (sumsForFactorCorrectionItr->second).first;
1502     */
1503     //--- new expression --------------
1504     if ( abs(sumOfWeights) > 0 )  
1505       factorCorrection = 1 + refR
1506     - (sumsForFactorCorrectionItr->second).second / sumOfWeights;
1507     //---------------------------------
1508     if ( correctionRMS/factorCorrection > MAX_REL_UNC_FACTOR ||
1509      nSubDetTracks < MIN_N_TRACKS_PER_CELL ) {
1510       correctionRMS = sqrt(pow(correctionRMS,2) + pow((factorCorrection - 1),2));
1511       factorCorrection = 1;
1512     }
1513     
1514     if( nSubDetTracks > MIN_N_TRACKS_PER_CELL ) {
1515       if (factorCorrection > 1) MeanConvergenceDelta += (1 - 1/factorCorrection);
1516       else                      MeanConvergenceDelta += (1 - factorCorrection);
1517       mkount++;
1518     }
1519 
1520               
1521     if (factors.find(detId) != factors.end()) {
1522       factorPrevious = factors[detId];
1523       factors[detId] *= factorCorrection;
1524       absErrorWprevious = uncFromWeights[detId];
1525       absErrorW = factorPrevious*correctionRMS;
1526       uncFromWeights[detId] = absErrorW;
1527       uncFromDeviation[detId] = factorPrevious*abs(factorCorrection - 1);
1528     }
1529     else {
1530       factorPrevious = 1;
1531       factors.insert(std::pair<unsigned int, double>(detId, factorCorrection));
1532       subDetector_final.insert(std::pair<unsigned int, double>(detId,
1533                                    subDetector_trk[detId]));
1534       absErrorW = correctionRMS;
1535       absErrorWprevious = 0;
1536       uncFromWeights.insert(std::pair<unsigned int, double>(detId, absErrorW));
1537       uncFromDeviation.insert(std::pair<unsigned int, double>(detId,
1538                                   abs(factorCorrection - 1)));
1539     }
1540 
1541     if ( debug > 0 ) {
1542       //if ( ieta[depth-1][kcur] == 27 && depth == 2 ) {
1543       std::cout.precision(3);
1544       std::cout << detId // << " (" << mkount << ")"
1545             << " *** ieta/depth | rw | cw | tw | fCor | nTrk | Ncell | C |::: "
1546             << ieta[depth-1][kcur] << "/" << depth << " | "
1547         << ratioWeights << " | "
1548             << sumOfWeights << " | "
1549             << (sumsForFactorCorrectionItr->second).second << " | "
1550             << factorCorrection << " | "
1551         << nSubDetTracks << " | "
1552         << NcellMean << " | "
1553         << correctionRMS << " |"
1554             << std::endl;
1555     }
1556 
1557     dets[kount] = detId;
1558 
1559     fac[depth-1][kcur] = factors[detId];
1560     dfac[depth-1][kcur] =
1561       sqrt(pow(uncFromWeights[detId],2) + pow(uncFromDeviation[detId],2));
1562     
1563     sys2statRatio[kount] = abs(factorPrevious*(factorCorrection - 1))/absErrorW;
1564     if ( sys2statRatio[kount] > MaxRatioUncertainties ) {
1565       MaxRatioUncertainties = sys2statRatio[kount];
1566       maxKountR = kount;
1567     }
1568     ztest[kount] = factorPrevious*(factorCorrection - 1)
1569       /sqrt(pow(absErrorWprevious,2) + pow(absErrorW,2));
1570     if ( abs(ztest[kount]) > MaxRelDeviationWeights ) {
1571       MaxRelDeviationWeights = abs(ztest[kount]);
1572       maxKountW = kount;
1573     } 
1574     kount++;
1575     kdep[depth-1]++;
1576   }
1577 
1578 //---- write current plots -----------------------------
1579   if ( debug > 0 ) std::cout << "Write graphs..." << std::endl;
1580 
1581   for ( unsigned ik = 0; ik < N_DEPTHS; ik++ ) {
1582     double x[MAXNUM_SUBDET], dx[MAXNUM_SUBDET], y[MAXNUM_SUBDET], dy[MAXNUM_SUBDET];
1583     for ( unsigned im = 0; im < MAXNUM_SUBDET; im++ ) {
1584       x[im] = ieta[ik][im];
1585       dx[im] = dieta[ik][im];
1586       y[im] = fac[ik][im];
1587       dy[im] = dfac[ik][im];
1588     }
1589     TGraphErrors*  g_fac = new TGraphErrors(kdep[ik], x, y, dx, dy);
1590     sprintf(name, "Extracted correction factors, depth %1d", ik+1);
1591     g_fac->SetTitle(name);
1592     g_fac->GetXaxis()->SetTitle("i#eta");
1593     sprintf(name, "Cfacs_depth%1d_%02d", ik+1, nIter);
1594     foutRootFile->WriteTObject(g_fac, name);
1595   }
1596 
1597   TGraph  *g_ztest, *g_sys2stat;
1598 
1599   g_ztest = new TGraph(kount, dets, ztest); 
1600   sprintf(name, "Z-test (unc. from weights) vs detId for iter %d", nIter);
1601   g_ztest->SetTitle(name);
1602   sprintf(name, "Ztest_detId_%02d", nIter);
1603   foutRootFile->WriteTObject(g_ztest, name);
1604 
1605   g_sys2stat = new TGraph(kount, dets, sys2statRatio); 
1606   sprintf(name, "Ratio of syst. to stat. unc. vs detId for iter %d", nIter);
1607   g_sys2stat->SetTitle(name);
1608   sprintf(name, "Sys2stat_detId_%02d", nIter);
1609   foutRootFile->WriteTObject(g_sys2stat, name);
1610 
1611   std::cout << "----------Iteration " << nIter << "--------------------" << std::endl;
1612   maxZtestFromWeights = MaxRelDeviationWeights; 
1613   std::cout << "Max abs(Z-test) with stat errors from weights = "
1614         << maxZtestFromWeights << " for subdetector " << maxKountW << std::endl;
1615   maxSys2StatRatio = MaxRatioUncertainties; 
1616   std::cout << "Max ratio of syst.(f_cur - f_prev) to stat. uncertainty = "
1617         << maxSys2StatRatio << " for subdetector " << maxKountR << std::endl;
1618 
1619   meanDeviation = (mkount > 0) ? (MeanConvergenceDelta/mkount) : 0;
1620   std::cout << "Mean absolute deviation from previous iteration = " << meanDeviation
1621         << " for " << mkount
1622         << " from " << kount << " DetIds" << std::endl;
1623 
1624 //--- delete hists ---------------------------
1625   for ( int i = 0; i < NUM_ETA_BINS; i++ ) {
1626     delete e2p[i];
1627   }
1628 
1629   return meanDeviation;
1630 }
1631 //**********************************************************
1632 // Last loop over events in the tree
1633 //**********************************************************
1634 Double_t CalibTree::lastLoop(unsigned int subsample,
1635                  unsigned int maxIter,
1636                  bool isTest,
1637                  unsigned int debug)
1638 {
1639   char name[100];
1640   unsigned int ndebug(0);
1641   
1642   char stest[80] = "test";
1643   if ( !isTest )
1644     sprintf(stest,"after %2d iterations", maxIter);
1645   char scorr[80] = "correction for PU";
1646   char sxlabel[80] ="(E^{cor}_{hcal} + E_{ecal})/p_{track}"; 
1647   if ( !APPLY_CORRECTION_FOR_PU ) {
1648     sprintf(scorr,"no correction for PU");
1649     sprintf(sxlabel,"(E_{hcal} + E_{ecal})/p_{track}");
1650   } 
1651     
1652   TF1* f1 = new TF1("f1","gaus", MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
1653 
1654   sprintf(name,"HB+HE: %s, %s", stest, scorr);
1655   e2p_last = new TH1F("e2p_last", name,
1656               NBIN_RESPONSE_HIST, MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
1657   e2p_last->Sumw2();
1658   e2p_last->GetXaxis()->SetTitle(sxlabel);
1659   
1660   sprintf(name,"HB: %s, %s", stest, scorr);
1661   e2pHB_last = new TH1F("e2pHB_last", name,
1662             NBIN_RESPONSE_HIST/2, MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
1663   e2pHB_last->Sumw2();
1664   e2pHB_last->GetXaxis()->SetTitle(sxlabel);
1665 
1666   sprintf(name,"Initial TR: %s", scorr);
1667   e2pTR_last = new TH1F("e2pTR_last", name,
1668             NBIN_RESPONSE_HIST/10, MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
1669   e2pTR_last->Sumw2();
1670   e2pTR_last->GetXaxis()->SetTitle(sxlabel);
1671   
1672   sprintf(name,"HE: %s, %s", stest, scorr);
1673   e2pHE_last = new TH1F("e2pHE_last", name,
1674             NBIN_RESPONSE_HIST/2, MIN_RESPONSE_HIST, MAX_RESPONSE_HIST);
1675   e2pHE_last->Sumw2();
1676   e2pHE_last->GetXaxis()->SetTitle(sxlabel);
1677  
1678 //--- initialize chain ----------------------------------------
1679   if (fChain == 0) return 0;
1680   Long64_t nentries = fChain->GetEntriesFast();
1681   Long64_t nb = 0;
1682   
1683   int nSelectedEvents(0);
1684 
1685   if ( debug > 0 ) { 
1686     std::cout << "------------- Last loop after " << maxIter << " iterations"
1687           << std::endl;
1688   }
1689 // ----------------------- loop over events -------------------------------------  
1690   for (Long64_t jentry=0; jentry<nentries; jentry++) {
1691     Long64_t ientry = LoadTree(jentry);
1692     if ( ientry < 0 || ndebug > debug ) break;   
1693     nb = fChain->GetEntry(jentry);   //nbytes += nb;
1694     
1695     if ( (jentry%2 == subsample) ) continue;   // only odd or even events
1696 
1697     // --------------- selection of good track --------------------
1698     
1699     if ( !goodTrack(t_ieta) ) continue;
1700 
1701     nSelectedEvents++;
1702     
1703     if ( debug > 1 ) {
1704       ndebug++;
1705       std::cout << "***Entry (Track) Number : " << ientry 
1706         << " p/eHCal/eMipDR/nDets : " << t_p << "/" << t_eHcal 
1707         << "/" << t_eMipDR << "/" << (*t_DetIds).size() 
1708         << std::endl;
1709     }
1710     
1711     double eTotal(0.0);
1712     double eTotalWithEcal(0.0);
1713       
1714     // ---- loop over active cells in the event for total energy ---
1715 
1716     for (unsigned int idet = 0; idet < (*t_DetIds).size(); idet++) { 
1717       double hitEnergy(0);  
1718       unsigned int detId = ( (*t_DetIds)[idet] & MASK ) | MASK2 ;
1719     
1720       if (factors.find(detId) != factors.end()) 
1721     hitEnergy = factors[detId] * (*t_HitEnergies)[idet];
1722       else 
1723     hitEnergy = (*t_HitEnergies)[idet];
1724 
1725       eTotal += hitEnergy;
1726     }
1727 
1728     eTotalWithEcal = eTotal + t_eMipDR;    
1729 
1730 // --- Correction for PU   --------      
1731     double eTotalCor(eTotal);
1732     double eTotalWithEcalCor(eTotalWithEcal);
1733     //double e10(0.0);
1734     //double e30(0.0);
1735     double correctionForPU(1.0);
1736     int abs_t_ieta = abs(t_ieta);
1737 
1738     if ( APPLY_CORRECTION_FOR_PU ) { 
1739       /* 
1740       for (unsigned int idet1 = 0; idet1 < (*t_DetIds1).size(); idet1++) { 
1741     double hitEnergy(0);
1742     unsigned int detId1 = ( (*t_DetIds1)[idet1] & MASK ) | MASK2;
1743     
1744     if (factors.find(detId1) != factors.end()) 
1745       hitEnergy = factors[detId1] * (*t_HitEnergies1)[idet1];
1746     else 
1747       hitEnergy = (*t_HitEnergies1)[idet1];
1748     
1749     e10 += hitEnergy;
1750       }
1751       for (unsigned int idet3 = 0; idet3 < (*t_DetIds3).size(); idet3++) { 
1752     double hitEnergy(0);
1753     unsigned int detId3 = ( (*t_DetIds3)[idet3] & MASK ) | MASK2;
1754     
1755     if (factors.find(detId3) != factors.end()) 
1756       hitEnergy = factors[detId3] * (*t_HitEnergies3)[idet3];
1757     else
1758       hitEnergy = (*t_HitEnergies3)[idet3];
1759     
1760     e30 += hitEnergy;
1761       }
1762       double de2p = (e30 - e10)/t_p;
1763       */
1764       
1765       double de2p = (t_eHcal30 - t_eHcal10)/t_p;
1766       if ( de2p > DELTA_CUT ) {
1767     int icor = int(abs_t_ieta >= FIRST_IETA_TR) + int(abs_t_ieta >= FIRST_IETA_HE)
1768       + int(abs_t_ieta >= FIRST_IETA_FWD_1) + int(abs_t_ieta >= FIRST_IETA_FWD_2);
1769     correctionForPU = (1 + LINEAR_COR_COEF[icor]*(t_eHcal/t_p)*de2p
1770                *(1 + SQUARE_COR_COEF[icor]*de2p));
1771       }
1772     }      
1773     // check for possibility to correct for PU
1774     if ( correctionForPU <= 0 || correctionForPU > 1 ) continue;
1775 
1776     eTotalCor = eTotal*correctionForPU;
1777     eTotalWithEcalCor = eTotalCor + t_eMipDR;
1778       
1779     e2p_last->Fill(eTotalWithEcalCor/t_p ,1.0);
1780 
1781     if ( abs_t_ieta < FIRST_IETA_TR )
1782       e2pHB_last->Fill(eTotalWithEcalCor/t_p ,1.0);
1783     else if ( abs_t_ieta < FIRST_IETA_HE )
1784       e2pTR_last->Fill(eTotalWithEcalCor/t_p ,1.0);
1785     else
1786       e2pHE_last->Fill(eTotalWithEcalCor/t_p ,1.0);
1787       
1788   } // ------------------- end of loop over events -------------------------------------
1789 
1790   if ( isTest ) {
1791     double fac[N_DEPTHS][MAXNUM_SUBDET]; 
1792     double dfac[N_DEPTHS][MAXNUM_SUBDET];
1793     double ieta[N_DEPTHS][MAXNUM_SUBDET];
1794     double dieta[N_DEPTHS][MAXNUM_SUBDET];
1795     unsigned int kdep[N_DEPTHS];
1796     for ( unsigned ik = 0; ik < N_DEPTHS; ik++ ) { kdep[ik] = 0; }
1797   
1798     std::map<unsigned int, double>::iterator factorsItr = factors.begin();
1799     for (factorsItr=factors.begin(); factorsItr != factors.end(); factorsItr++){
1800 
1801       unsigned int detId = factorsItr->first;
1802       int zside = (detId&ZSIDE_MASK) ? 1 : -1;
1803       int depth = ((detId>>DEPTH_OFFSET)&DEPTH_MASK) + int(MERGE_PHI_AND_DEPTHS);
1804 
1805       unsigned int kcur = kdep[depth-1];
1806       ieta[depth-1][kcur] = int((detId>>ETA_OFFSET) & ETA_MASK)*zside;
1807       dieta[depth-1][kcur] = 0;
1808       fac[depth-1][kcur] = factorsItr->second;
1809       dfac[depth-1][kcur] = 0;
1810       kdep[depth-1]++;
1811     }
1812     for ( unsigned ik = 0; ik < N_DEPTHS; ik++ ) {
1813       double x[MAXNUM_SUBDET], dx[MAXNUM_SUBDET], y[MAXNUM_SUBDET], dy[MAXNUM_SUBDET];
1814       for ( unsigned im = 0; im < MAXNUM_SUBDET; im++ ) {
1815     x[im] = ieta[ik][im];
1816     dx[im] = dieta[ik][im];
1817     y[im] = fac[ik][im];
1818     dy[im] = dfac[ik][im];
1819       }
1820       TGraphErrors*  g_fac = new TGraphErrors(kdep[ik], x, y, dx, dy);
1821       sprintf(name, "Applied correction factors, depth %1d", ik+1);
1822       g_fac->SetTitle(name);
1823       g_fac->GetXaxis()->SetTitle("i#eta");
1824       sprintf(name, "Cfacs_depth%1d", ik+1);
1825       foutRootFile->WriteTObject(g_fac, name);
1826     }
1827   }
1828   //--- fit response distributions ---------------------------------
1829 
1830   double xl = e2p_last->GetMean() - FIT_RMS_INTERVAL*e2p_last->GetRMS();
1831   double xr = e2p_last->GetMean() + FIT_RMS_INTERVAL*e2p_last->GetRMS();
1832   e2p_last->Fit("f1","QN", "R", xl, xr);
1833   xl = f1->GetParameter(1) - FIT_RMS_INTERVAL*f1->GetParameter(2);
1834   xr = f1->GetParameter(1) + FIT_RMS_INTERVAL*f1->GetParameter(2);
1835   e2p_last->Fit("f1","QN", "R", xl, xr);
1836 
1837   double fitMPV = f1->GetParameter(1);
1838   /*
1839   xl = e2pHB_last->GetMean() - FIT_RMS_INTERVAL*e2pHB_last->GetRMS();
1840   xr = e2pHB_last->GetMean() + FIT_RMS_INTERVAL*e2pHB_last->GetRMS();
1841   e2pHB_last->Fit("f1","QN", "R", xl, xr);
1842 
1843   xl = e2pHE_last->GetMean() - FIT_RMS_INTERVAL*e2pHE_last->GetRMS();
1844   xr = e2pHE_last->GetMean() + FIT_RMS_INTERVAL*e2pHE_last->GetRMS();
1845   e2pHE_last->Fit("f1","QN", "R", xl, xr);
1846   */
1847   return fitMPV;
1848 }
1849 
1850 //**********************************************************
1851 // Isolated track selection
1852 //**********************************************************
1853 Bool_t CalibTree::goodTrack(int ieta) 
1854 {
1855 
1856   double maxCharIso = limCharIso*exp(abs(ieta)*constForFlexSel);
1857 
1858   bool ok = (    (t_selectTk)
1859           && (t_qltyMissFlag)
1860           && (t_hmaxNearP < maxCharIso)
1861           && (t_eMipDR < limMipEcal) 
1862           && (t_p > minTrackMom) && (t_p < maxTrackMom)
1863           && (t_pt >= minTrackPt)               // constraint on track pt
1864           && (t_eHcal >= minEnrHcal)            // constraint on Hcal energy
1865           && (t_eHcal/t_p < UPPER_LIMIT_RESPONSE_BEFORE_COR)
1866          // reject events with too big cluster energy
1867          //&& ((t_eHcal30 - t_eHcal10)/t_p < UPPER_LIMIT_DELTA_PU_COR)
1868          // reject events with too high PU in the ring around cluster
1869          );
1870   return ok;
1871 }
1872 //**********************************************************
1873 // Save txt file with calculated factors
1874 //**********************************************************
1875 unsigned int CalibTree::saveFactorsInFile(std::string txtFileName)
1876 {
1877   char sprnt[100];
1878   
1879   FILE* foutTxtFile = fopen(txtFileName.c_str(),"w+");
1880   fprintf(foutTxtFile,
1881       "%1s%16s%16s%16s%9s%11s\n","#", "eta", "depth", "det", "value", "DetId");
1882 
1883   std::cout << "New factors:" << std::endl;
1884   std::map<unsigned int, double>::iterator factorsItr = factors.begin();
1885   unsigned int indx(0);
1886   unsigned int isave(0);
1887   
1888   for (factorsItr=factors.begin(); factorsItr != factors.end(); factorsItr++, indx++){
1889     unsigned int detId = factorsItr->first;
1890     int ieta = (detId>>ETA_OFFSET) & ETA_MASK;
1891     int zside= (detId&ZSIDE_MASK) ? 1 : -1;
1892     int depth= ((detId>>DEPTH_OFFSET)&DEPTH_MASK) + int(MERGE_PHI_AND_DEPTHS);
1893 
1894     double erWeight = 100*uncFromWeights[detId]/factorsItr->second;
1895     double erDev = 100*uncFromDeviation[detId]/factorsItr->second;
1896     double erTotal = 100*sqrt(pow(uncFromWeights[detId],2)
1897               + pow(uncFromDeviation[detId],2))/factorsItr->second;
1898     
1899     if ( N_ETA_RINGS_PER_BIN < 2 ) { 
1900       sprintf(sprnt,
1901           "DetId[%3d] %x (%3d,%1d)  %6.4f  : %6d  [%8.3f%% + %8.3f%% = %8.3f%%]",
1902           indx, detId, ieta*zside, depth,
1903           factorsItr->second, nTrks[detId],
1904           erWeight, erDev, erTotal);
1905       std::cout << sprnt << std::endl;
1906     }
1907     else {
1908       int ieta_min = ieta - (N_ETA_RINGS_PER_BIN - 1);
1909       sprintf(sprnt,
1910           "DetId[%3d] %x (%3d:%3d,%1d)  %6.4f  : %6d  [%8.3f%% + %8.3f%% = %8.3f%%]",
1911           indx, detId, ieta_min*zside, ieta*zside, depth,
1912           factorsItr->second, nTrks[detId],
1913           erWeight, erDev, erTotal);
1914       std::cout << sprnt << std::endl;
1915     }
1916         /*
1917     std::cout << "DetId[" << indx << "] " << std::hex  << (detId) << std::dec 
1918           << "(" << ieta*zside << "," << depth << ") ( nTrks:" 
1919           << nTrks[detId] << ") : " << factorsItr->second
1920           << ""
1921           << std::endl;
1922         */
1923     
1924     const char* subDetector[2] = {"HB","HE"};
1925     if ( nTrks[detId] < MIN_N_TRACKS_PER_CELL ) continue;
1926     isave++;
1927     fprintf(foutTxtFile, "%17i%16i%16s%9.5f%11X\n", 
1928         ieta*zside, depth, subDetector[subDetector_final[detId]-1],
1929         factorsItr->second, detId);
1930   }
1931   fclose(foutTxtFile);
1932   foutTxtFile = NULL;
1933   return isave;
1934 }
1935 //**********************************************************
1936 // Get factors from txt file
1937 //**********************************************************
1938 Bool_t CalibTree::getFactorsFromFile(std::string txtFileName,
1939                      unsigned int dbg)
1940 {
1941 
1942   if ( !gSystem->Which("./", txtFileName.c_str() ) ) return false;
1943 
1944   FILE* finTxtFile = fopen(txtFileName.c_str(),"r");
1945   int flag;
1946 
1947   char header[80]; 
1948   for ( unsigned int i = 0; i < 6; i++ ) { 
1949     flag = fscanf(finTxtFile, "%7s", header);
1950   }
1951 
1952   int eta;
1953   int depth;
1954   char det[2]; 
1955   double cellFactor;
1956   unsigned int detId;
1957   unsigned int nReadFactors(0);
1958   
1959   while ( fscanf(finTxtFile, "%3d", &eta) != EOF )
1960     {
1961       flag = fscanf(finTxtFile, "%2d", &depth);
1962       flag = fscanf(finTxtFile, "%10s", det);
1963       flag = fscanf(finTxtFile, "%lf", &cellFactor);
1964       flag = fscanf(finTxtFile, "%x", &detId);
1965       factors.insert( std::pair<unsigned int, double>(detId, cellFactor) );
1966       nReadFactors++;
1967       if ( dbg > 0 ) 
1968     std::cout << "  " << std::dec << cellFactor
1969           << "  " << std::hex << detId << std::endl; 
1970     }
1971 
1972   std::cout << std::dec << nReadFactors << " factors read from file "
1973         << txtFileName
1974         << std::endl;
1975   
1976   return true;
1977 }
1978 //**********************************************************
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!