Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:51 

0001 //////////////////////////////////////////////////////////////////////////////
0002 //
0003 // Usage:
0004 // .L CalibRho.C+g
0005 //
0006 //  EHcalVsRho c1(inpFilName, dupFileName);
0007 //  c1.LoopFill(maxEta, outFile1, logFile);       Fills E vs rho plots
0008 //  FitEvsRho(logFile, parameterFile, outFile2);  Fits the area vs iEta
0009 //  c1.LoopTest(maxEta, parameterFile, rootFile); Makes the corrected histos
0010 //  FitEovPwithRho(maxEta, rootFile, outFile3);   Fits E/p plots
0011 //  PlotEvsRho(inFile, eta, type, save);          Make the plots
0012 //
0013 //  inpFileName   (const char*)  File name of the input ROOT tree
0014 //                               or name of the file containing a list of
0015 //                               file names of input ROOT trees
0016 //  dupFileName   (const char*)  Name of the file containing list of entries
0017 //                               of duplicate events
0018 //  maxEta        (int)          Maximum value of |iEta|
0019 //  outFile1      (const char*)  Output ROOT file name which will contain the
0020 //                               scatter plots and profile histograms which
0021 //                               will provide the estimate of "effective area"
0022 //  logFile       (const char*)  Name of the text file which will contain the
0023 //                               effective area for each ieta value
0024 //  parameterFile (const char*)  Name of the text file with values of the
0025 //                               fitted parameter set
0026 //  outFile2      (const char*)  Name of the ROOT file with the results of
0027 //                               the fit
0028 //  rootFile      (const char*)  Name of the ROOT file with corrected and
0029 //                               uncorrected histograms of E/p
0030 //  outFile3      (const char*)  Name of the ROOT file containing the E/p
0031 //                               histograms and the fit results
0032 //  eta           (int)          Plot the histograms for the given iEta value
0033 //                               (if it is 0; plots for all iEta's)
0034 //  type          (int)          Plot E vs Rho (type=0) or E/p fits (type=1)
0035 //  save          (bool)         Save the plots as pdf file
0036 //////////////////////////////////////////////////////////////////////////////
0037 
0038 #include <TCanvas.h>
0039 #include <TChain.h>
0040 #include <TFile.h>
0041 #include <TFitResult.h>
0042 #include <TFitResultPtr.h>
0043 #include <TGraph.h>
0044 #include <TGraphErrors.h>
0045 #include <TGraphAsymmErrors.h>
0046 #include <TH1.h>
0047 #include <TH2.h>
0048 #include <TMultiGraph.h>
0049 #include <TPaveStats.h>
0050 #include <TPaveText.h>
0051 #include <TProfile.h>
0052 #include <TROOT.h>
0053 #include <TStyle.h>
0054 
0055 #include <algorithm>
0056 #include <iostream>
0057 #include <fstream>
0058 #include <vector>
0059 
0060 #include "CalibCorr.C"
0061 
0062 class EHcalVsRho {
0063 public:
0064   TChain *fChain;  //!pointer to the analyzed TTree or TChain
0065   Int_t fCurrent;  //!current Tree number in a TChain
0066 
0067   EHcalVsRho(const char *inFile, const char *dupFile);
0068   virtual ~EHcalVsRho();
0069   virtual Int_t Cut(Long64_t entry);
0070   virtual Int_t GetEntry(Long64_t entry);
0071   virtual Long64_t LoadTree(Long64_t entry);
0072   virtual void Init(TChain *tree, const char *dupFile);
0073   virtual Bool_t Notify();
0074   virtual void Show(Long64_t entry = -1);
0075   void LoopFill(int maxEta, const char *outFile, const char *logFile);
0076   void LoopTest(int maxEta, const char *inFile, const char *outFile);
0077   double EffCalc(TH1D *, double, double &, double &);
0078   double getEA(const int ieta, const double *par);
0079 
0080 private:
0081   // Declaration of leaf types
0082   Int_t t_Run;
0083   Int_t t_Event;
0084   Int_t t_DataType;
0085   Int_t t_ieta;
0086   Int_t t_iphi;
0087   Double_t t_EventWeight;
0088   Int_t t_nVtx;
0089   Int_t t_nTrk;
0090   Int_t t_goodPV;
0091   Double_t t_l1pt;
0092   Double_t t_l1eta;
0093   Double_t t_l1phi;
0094   Double_t t_l3pt;
0095   Double_t t_l3eta;
0096   Double_t t_l3phi;
0097   Double_t t_p;
0098   Double_t t_pt;
0099   Double_t t_phi;
0100   Double_t t_mindR1;
0101   Double_t t_mindR2;
0102   Double_t t_eMipDR;
0103   Double_t t_eHcal;
0104   Double_t t_eHcal10;
0105   Double_t t_eHcal30;
0106   Double_t t_hmaxNearP;
0107   Double_t t_rhoh;
0108   Bool_t t_selectTk;
0109   Bool_t t_qltyFlag;
0110   Bool_t t_qltyMissFlag;
0111   Bool_t t_qltyPVFlag;
0112   Double_t t_gentrackP;
0113   std::vector<unsigned int> *t_DetIds;
0114   std::vector<double> *t_HitEnergies;
0115   std::vector<bool> *t_trgbits;
0116   std::vector<unsigned int> *t_DetIds1;
0117   std::vector<unsigned int> *t_DetIds3;
0118   std::vector<double> *t_HitEnergies1;
0119   std::vector<double> *t_HitEnergies3;
0120 
0121   // List of branches
0122   TBranch *b_t_Run;           //!
0123   TBranch *b_t_Event;         //!
0124   TBranch *b_t_DataType;      //!
0125   TBranch *b_t_ieta;          //!
0126   TBranch *b_t_iphi;          //!
0127   TBranch *b_t_EventWeight;   //!
0128   TBranch *b_t_nVtx;          //!
0129   TBranch *b_t_nTrk;          //!
0130   TBranch *b_t_goodPV;        //!
0131   TBranch *b_t_l1pt;          //!
0132   TBranch *b_t_l1eta;         //!
0133   TBranch *b_t_l1phi;         //!
0134   TBranch *b_t_l3pt;          //!
0135   TBranch *b_t_l3eta;         //!
0136   TBranch *b_t_l3phi;         //!
0137   TBranch *b_t_p;             //!
0138   TBranch *b_t_pt;            //!
0139   TBranch *b_t_phi;           //!
0140   TBranch *b_t_mindR1;        //!
0141   TBranch *b_t_mindR2;        //!
0142   TBranch *b_t_eMipDR;        //!
0143   TBranch *b_t_eHcal;         //!
0144   TBranch *b_t_eHcal10;       //!
0145   TBranch *b_t_eHcal30;       //!
0146   TBranch *b_t_hmaxNearP;     //!
0147   TBranch *b_t_rhoh;          //!
0148   TBranch *b_t_selectTk;      //!
0149   TBranch *b_t_qltyFlag;      //!
0150   TBranch *b_t_qltyMissFlag;  //!
0151   TBranch *b_t_qltyPVFlag;    //!
0152   TBranch *b_t_gentrackP;     //!
0153   TBranch *b_t_DetIds;        //!
0154   TBranch *b_t_HitEnergies;   //!
0155   TBranch *b_t_trgbits;       //!
0156   TBranch *b_t_DetIds1;       //!
0157   TBranch *b_t_DetIds3;       //!
0158   TBranch *b_t_HitEnergies1;  //!
0159   TBranch *b_t_HitEnergies3;  //!
0160 
0161   std::vector<Long64_t> entries_;
0162 };
0163 
0164 EHcalVsRho::EHcalVsRho(const char *inFile, const char *dupFile) : fChain(0) {
0165   char treeName[400];
0166   sprintf(treeName, "HcalIsoTrkAnalyzer/CalibTree");
0167   TChain *chain = new TChain(treeName);
0168   std::cout << "Create a chain for " << treeName << " from " << inFile << std::endl;
0169   if (!fillChain(chain, inFile)) {
0170     std::cout << "*****No valid tree chain can be obtained*****" << std::endl;
0171   } else {
0172     std::cout << "Proceed with a tree chain with " << chain->GetEntries() << " entries" << std::endl;
0173     Init(chain, dupFile);
0174   }
0175 }
0176 
0177 EHcalVsRho::~EHcalVsRho() {
0178   if (!fChain)
0179     return;
0180   delete fChain->GetCurrentFile();
0181 }
0182 
0183 Int_t EHcalVsRho::GetEntry(Long64_t entry) {
0184   // Read contents of entry.
0185   if (!fChain)
0186     return 0;
0187   return fChain->GetEntry(entry);
0188 }
0189 
0190 Long64_t EHcalVsRho::LoadTree(Long64_t entry) {
0191   // Set the environment to read one entry
0192   if (!fChain)
0193     return -5;
0194   Long64_t centry = fChain->LoadTree(entry);
0195   if (centry < 0)
0196     return centry;
0197   if (fChain->GetTreeNumber() != fCurrent) {
0198     fCurrent = fChain->GetTreeNumber();
0199     Notify();
0200   }
0201   return centry;
0202 }
0203 
0204 void EHcalVsRho::Init(TChain *tree, const char *dupFile) {
0205   // The Init() function is called when the selector needs to initialize
0206   // a new tree or chain. Typically here the branch addresses and branch
0207   // pointers of the tree will be set.
0208   // It is normally not necessary to make changes to the generated
0209   // code, but the routine can be extended by the user if needed.
0210   // Init() will be called many times when running on PROOF
0211   // (once per file to be processed).
0212 
0213   // Set object pointer
0214   t_DetIds = 0;
0215   t_HitEnergies = 0;
0216   t_trgbits = 0;
0217   t_DetIds1 = 0;
0218   t_DetIds3 = 0;
0219   t_HitEnergies1 = 0;
0220   t_HitEnergies3 = 0;
0221   // Set branch addresses and branch pointers
0222   if (!tree)
0223     return;
0224   fChain = tree;
0225   fCurrent = -1;
0226   fChain->SetMakeClass(1);
0227 
0228   fChain->SetBranchAddress("t_Run", &t_Run, &b_t_Run);
0229   fChain->SetBranchAddress("t_Event", &t_Event, &b_t_Event);
0230   fChain->SetBranchAddress("t_DataType", &t_DataType, &b_t_DataType);
0231   fChain->SetBranchAddress("t_ieta", &t_ieta, &b_t_ieta);
0232   fChain->SetBranchAddress("t_iphi", &t_iphi, &b_t_iphi);
0233   fChain->SetBranchAddress("t_EventWeight", &t_EventWeight, &b_t_EventWeight);
0234   fChain->SetBranchAddress("t_nVtx", &t_nVtx, &b_t_nVtx);
0235   fChain->SetBranchAddress("t_nTrk", &t_nTrk, &b_t_nTrk);
0236   fChain->SetBranchAddress("t_goodPV", &t_goodPV, &b_t_goodPV);
0237   fChain->SetBranchAddress("t_l1pt", &t_l1pt, &b_t_l1pt);
0238   fChain->SetBranchAddress("t_l1eta", &t_l1eta, &b_t_l1eta);
0239   fChain->SetBranchAddress("t_l1phi", &t_l1phi, &b_t_l1phi);
0240   fChain->SetBranchAddress("t_l3pt", &t_l3pt, &b_t_l3pt);
0241   fChain->SetBranchAddress("t_l3eta", &t_l3eta, &b_t_l3eta);
0242   fChain->SetBranchAddress("t_l3phi", &t_l3phi, &b_t_l3phi);
0243   fChain->SetBranchAddress("t_p", &t_p, &b_t_p);
0244   fChain->SetBranchAddress("t_pt", &t_pt, &b_t_pt);
0245   fChain->SetBranchAddress("t_phi", &t_phi, &b_t_phi);
0246   fChain->SetBranchAddress("t_mindR1", &t_mindR1, &b_t_mindR1);
0247   fChain->SetBranchAddress("t_mindR2", &t_mindR2, &b_t_mindR2);
0248   fChain->SetBranchAddress("t_eMipDR", &t_eMipDR, &b_t_eMipDR);
0249   fChain->SetBranchAddress("t_eHcal", &t_eHcal, &b_t_eHcal);
0250   fChain->SetBranchAddress("t_eHcal10", &t_eHcal10, &b_t_eHcal10);
0251   fChain->SetBranchAddress("t_eHcal30", &t_eHcal30, &b_t_eHcal30);
0252   fChain->SetBranchAddress("t_hmaxNearP", &t_hmaxNearP, &b_t_hmaxNearP);
0253   fChain->SetBranchAddress("t_rhoh", &t_rhoh, &b_t_rhoh);
0254   fChain->SetBranchAddress("t_selectTk", &t_selectTk, &b_t_selectTk);
0255   fChain->SetBranchAddress("t_qltyFlag", &t_qltyFlag, &b_t_qltyFlag);
0256   fChain->SetBranchAddress("t_qltyMissFlag", &t_qltyMissFlag, &b_t_qltyMissFlag);
0257   fChain->SetBranchAddress("t_qltyPVFlag", &t_qltyPVFlag, &b_t_qltyPVFlag);
0258   fChain->SetBranchAddress("t_gentrackP", &t_gentrackP, &b_t_gentrackP);
0259   fChain->SetBranchAddress("t_DetIds", &t_DetIds, &b_t_DetIds);
0260   fChain->SetBranchAddress("t_HitEnergies", &t_HitEnergies, &b_t_HitEnergies);
0261   fChain->SetBranchAddress("t_trgbits", &t_trgbits, &b_t_trgbits);
0262   fChain->SetBranchAddress("t_DetIds1", &t_DetIds1, &b_t_DetIds1);
0263   fChain->SetBranchAddress("t_DetIds3", &t_DetIds3, &b_t_DetIds3);
0264   fChain->SetBranchAddress("t_HitEnergies1", &t_HitEnergies1, &b_t_HitEnergies1);
0265   fChain->SetBranchAddress("t_HitEnergies3", &t_HitEnergies3, &b_t_HitEnergies3);
0266   Notify();
0267 
0268   if (std::string(dupFile) != "") {
0269     ifstream infile(dupFile);
0270     if (!infile.is_open()) {
0271       std::cout << "Cannot open " << dupFile << std::endl;
0272     } else {
0273       while (1) {
0274         Long64_t jentry;
0275         infile >> jentry;
0276         if (!infile.good())
0277           break;
0278         entries_.push_back(jentry);
0279       }
0280       infile.close();
0281       std::cout << "Reads a list of " << entries_.size() << " events from " << dupFile << std::endl;
0282     }
0283   }
0284 }
0285 
0286 Bool_t EHcalVsRho::Notify() {
0287   // The Notify() function is called when a new file is opened. This
0288   // can be either for a new TTree in a TChain or when when a new TTree
0289   // is started when using PROOF. It is normally not necessary to make changes
0290   // to the generated code, but the routine can be extended by the
0291   // user if needed. The return value is currently not used.
0292 
0293   return kTRUE;
0294 }
0295 
0296 void EHcalVsRho::Show(Long64_t entry) {
0297   // Print contents of entry.
0298   // If entry is not specified, print current entry
0299   if (!fChain)
0300     return;
0301   fChain->Show(entry);
0302 }
0303 
0304 Int_t EHcalVsRho::Cut(Long64_t) {
0305   // This function may be called from Loop.
0306   // returns  1 if entry is accepted.
0307   // returns -1 otherwise.
0308   return 1;
0309 }
0310 
0311 void EHcalVsRho::LoopFill(int maxEta, const char *outFile, const char *logFile) {
0312   if (fChain == 0)
0313     return;
0314   TFile *f1 = new TFile(outFile, "RECREATE");
0315   char name[100], Title[100], graph[100], proji[100];
0316 
0317   std::vector<TH2D *> VIsoRho;
0318   std::vector<TProfile *> Hcal_corr;
0319   for (int ieta = 0; ieta < maxEta; ieta++) {
0320     sprintf(name, "IsoRho2d%d", ieta + 1);
0321     sprintf(Title, "Iso vs Rho %d", ieta + 1);
0322     VIsoRho.push_back(new TH2D(name, Title, 30, 0, 30, 25000, 0, 250));
0323     sprintf(name, "IsoRhoProfile%d", ieta + 1);
0324     Hcal_corr.push_back(new TProfile(name, Title, 30, 0, 30));
0325   }
0326 
0327   Long64_t nentries = fChain->GetEntriesFast();
0328   std::cout << "Total # of entries: " << nentries << std::endl;
0329   Long64_t nbytes = 0, nb = 0;
0330   Long64_t kount(0), duplicate(0), good(0);
0331   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0332     Long64_t ientry = LoadTree(jentry);
0333     if (ientry < 0)
0334       break;
0335     nb = fChain->GetEntry(jentry);
0336     nbytes += nb;
0337 
0338     ++kount;
0339     if (kount % 100000 == 0)
0340       std::cout << "Processing Entry " << kount << std::endl;
0341     bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
0342     if (!select) {
0343       ++duplicate;
0344       continue;
0345     }
0346 
0347     int absIeta = abs(t_ieta);
0348     if ((absIeta <= maxEta) && (t_p >= 40) && (t_p <= 60)) {
0349       VIsoRho[absIeta - 1]->Fill(t_rhoh, t_eHcal);
0350       Hcal_corr[absIeta - 1]->Fill(t_rhoh, t_eHcal);
0351       ++good;
0352     }
0353   }
0354   std::cout << "Uses " << good << " events out of " << kount << " excluding " << duplicate << " duplicate events"
0355             << std::endl;
0356 
0357   gStyle->SetCanvasBorderMode(0);
0358   gStyle->SetCanvasColor(kWhite);
0359   gStyle->SetPadColor(kWhite);
0360   gStyle->SetFillColor(kWhite);
0361   gStyle->SetOptTitle(0);
0362   gStyle->SetOptFit(1);
0363   std::ofstream myfile;
0364   myfile.open(logFile);
0365   for (int ieta = 0; ieta < maxEta; ieta++) {
0366     VIsoRho[ieta]->Write();
0367     Hcal_corr[ieta]->Write();
0368     TH2D *his_i = dynamic_cast<TH2D *>(VIsoRho[ieta]->Clone());
0369     his_i->GetEntries();
0370     int dim = his_i->GetXaxis()->GetNbins();
0371     double *xcut, *binc, *errXL, *errXH, *errYL, *errYH;
0372     double *errX, *errY;
0373     double errX1, errX2, xmax(0);
0374 
0375     xcut = new double[dim];
0376     binc = new double[dim];
0377     errX = new double[dim];
0378     errY = new double[dim];
0379     errXL = new double[dim];
0380     errXH = new double[dim];
0381     errYL = new double[dim];
0382     errYH = new double[dim];
0383 
0384     for (int j = 0; j < dim; j++) {
0385       sprintf(proji, "proj%d-%d", ieta + 1, j);
0386       TH1D *h_proj = dynamic_cast<TH1D *>(his_i->ProjectionY(proji, j, j + 1, " "));
0387       binc[j] = his_i->GetXaxis()->GetBinCenter(j + 1);
0388       xcut[j] = EffCalc(h_proj, 0.90, errX1, errX2);
0389 
0390       errXL[j] = 0.0;
0391       errXH[j] = 0.0;
0392       errYL[j] = errX1;
0393       errYH[j] = errX2;
0394 
0395       errX[j] = 0.0;
0396       errY[j] = 0.0;
0397       h_proj->Write();
0398       if (xcut[j] > xmax)
0399         xmax = xcut[j];
0400     }
0401 
0402     TGraphAsymmErrors *Isovsrho = new TGraphAsymmErrors(dim, binc, xcut, errXL, errXH, errYL, errYH);
0403     sprintf(graph, "IsovsRho%d", ieta + 1);
0404     sprintf(name, "EvsRho%d", ieta + 1);
0405 
0406     TF1 *fnc = new TF1("fnc", "[1]*x + [0]", 4, 13);
0407     TFitResultPtr fitI = Isovsrho->Fit("fnc", "+QSR");
0408     double ic = fnc->GetParameter(1);
0409     double err = fitI->FitResult::Error(1);
0410     myfile << ieta + 1 << " " << ic << " " << err << std::endl;
0411     std::cout << "Fit " << ieta + 1 << " " << fnc->GetParameter(0) << " " << fitI->FitResult::Error(0) << " " << ic
0412               << " " << err << "\n";
0413     gStyle->SetOptFit(1);
0414     TCanvas *pad = new TCanvas(graph, name, 0, 10, 1200, 400);
0415     pad->SetRightMargin(0.10);
0416     pad->SetTopMargin(0.10);
0417     Isovsrho->SetMarkerStyle(24);
0418     Isovsrho->SetMarkerSize(0.4);
0419     Isovsrho->GetXaxis()->SetRangeUser(0, 15);
0420     Isovsrho->GetXaxis()->SetTitle("#rho");
0421     Isovsrho->GetXaxis()->SetLabelSize(0.04);
0422     Isovsrho->GetXaxis()->SetTitleSize(0.06);
0423     Isovsrho->GetXaxis()->SetTitleOffset(0.8);
0424     Isovsrho->GetYaxis()->SetRangeUser(0, 1.25 * xmax);
0425     Isovsrho->GetYaxis()->SetTitle("Energy (GeV)");
0426     Isovsrho->GetYaxis()->SetLabelSize(0.04);
0427     Isovsrho->GetYaxis()->SetTitleSize(0.06);
0428     Isovsrho->GetYaxis()->SetTitleOffset(0.6);
0429     Isovsrho->Draw("AP");
0430     pad->Update();
0431     TPaveStats *st1 = (TPaveStats *)Isovsrho->GetListOfFunctions()->FindObject("stats");
0432     if (st1 != nullptr) {
0433       st1->SetY1NDC(0.78);
0434       st1->SetY2NDC(0.90);
0435       st1->SetX1NDC(0.65);
0436       st1->SetX2NDC(0.90);
0437     }
0438     pad->Write();
0439   }
0440 
0441   myfile.close();
0442   f1->Close();
0443 }
0444 
0445 double EHcalVsRho::EffCalc(TH1D *h, double perc, double &errXL, double &errXH) {
0446   double eff, eff_err = 0.0, xCut = 0.0;
0447   int tot = h->GetEntries();
0448   int integ = 0;
0449   errXL = 0.0;
0450   errXH = 0.0;
0451   for (int i = 0; i < (h->GetXaxis()->GetNbins() + 1); i++) {
0452     xCut = h->GetXaxis()->GetBinLowEdge(i);
0453     integ += h->GetBinContent(i);
0454 
0455     if (integ != 0 && tot != 0) {
0456       eff = (integ * 1.0 / tot);
0457       eff_err = sqrt(eff * (1 - eff) / tot);
0458     } else {
0459       eff = 0.0;
0460     }
0461     if (eff > perc)
0462       break;
0463   }
0464   if (eff == 0.0)
0465     xCut = 0.0;
0466   errXL = eff_err;
0467   errXH = eff_err;
0468   return xCut;
0469 }
0470 
0471 void EHcalVsRho::LoopTest(int maxEta, const char *inFile, const char *outFile) {
0472   if (fChain == 0)
0473     return;
0474 
0475   TFile *f1 = new TFile(outFile, "RECREATE");
0476   std::map<int, TH1D *> histo, histo_uncorr;
0477   char name[100], title[100];
0478   for (int ieta = -maxEta; ieta <= maxEta; ieta++) {
0479     sprintf(name, "MPV%d", ieta);
0480     sprintf(title, "Corrected Response (i#eta = %d)", ieta - 30);
0481     histo[ieta] = new TH1D(name, title, 100, 0, 2);
0482     sprintf(name, "MPVUn%d", ieta);
0483     sprintf(title, "Uncorrected Response (i#eta = %d)", ieta - 30);
0484     histo_uncorr[ieta] = new TH1D(name, title, 100, 0, 2);
0485   }
0486   std::cout << "Initialized histograms from " << -maxEta << ":" << maxEta << "\n";
0487 
0488   double par[10];
0489   ifstream myReadFile;
0490   myReadFile.open(inFile);
0491   int npar = 0;
0492   if (myReadFile.is_open()) {
0493     while (!myReadFile.eof()) {
0494       myReadFile >> par[npar];
0495       ++npar;
0496     }
0497   }
0498   myReadFile.close();
0499   std::cout << "Reads " << npar << " parameters:";
0500   for (int k = 0; k < npar; ++k)
0501     std::cout << " [" << k << "] = " << par[k];
0502   std::cout << std::endl;
0503 
0504   Long64_t nentries = fChain->GetEntriesFast();
0505   Long64_t nbytes = 0, nb = 0;
0506   Long64_t kount(0), duplicate(0), good(0);
0507   for (Long64_t jentry = 0; jentry < nentries; jentry++) {
0508     Long64_t ientry = LoadTree(jentry);
0509     if (ientry < 0)
0510       break;
0511     nb = fChain->GetEntry(jentry);
0512     nbytes += nb;
0513     ++kount;
0514     if (kount % 100000 == 0)
0515       std::cout << "Processing Entry " << kount << std::endl;
0516     bool select = (std::find(entries_.begin(), entries_.end(), jentry) == entries_.end());
0517     if (!select) {
0518       ++duplicate;
0519       continue;
0520     }
0521 
0522     select = ((t_qltyFlag) && (t_selectTk) && (t_hmaxNearP < 10.0) && (t_eMipDR < 1.0) && (t_p > 40) && (t_p < 60.0));
0523     if (select) {
0524       double corr_eHcal = 0.0;
0525       int absIeta = abs(t_ieta);
0526       ++good;
0527       if (absIeta <= maxEta) {
0528         corr_eHcal = t_eHcal - t_rhoh * getEA(absIeta, par);
0529         double myEovP = corr_eHcal / (t_p - t_eMipDR);
0530         double myEovP_uncorr = t_eHcal / (t_p - t_eMipDR);
0531         histo[t_ieta]->Fill(myEovP);
0532         histo_uncorr[t_ieta]->Fill(myEovP_uncorr);
0533       }
0534     }
0535   }
0536 
0537   for (std::map<int, TH1D *>::iterator itr = histo.begin(); itr != histo.end(); ++itr)
0538     itr->second->Write();
0539   for (std::map<int, TH1D *>::iterator itr = histo_uncorr.begin(); itr != histo_uncorr.end(); ++itr)
0540     itr->second->Write();
0541   f1->Close();
0542   std::cout << "Processes " << good << " out of " << kount << " events with " << duplicate << " duplicate entries"
0543             << std::endl;
0544 }
0545 
0546 double EHcalVsRho::getEA(const int eta, const double *par) {
0547   double eA;
0548   if (eta < 20)
0549     eA = par[0];
0550   else
0551     eA = (((par[5] * eta + par[4]) * eta + par[3]) * eta + par[2]) * eta + par[1];
0552   return eA;
0553 }
0554 
0555 void FitEvsRho(const char *inFile, const char *outFile, const char *rootFile) {
0556   const int ndim = 30;
0557   double EA[ndim] = {0.0};
0558   double errEA[ndim] = {0.0};
0559   double ietaEA[ndim] = {0.0};
0560   ifstream myReadFile;
0561   myReadFile.open(inFile);
0562 
0563   int ii = 0;
0564   if (myReadFile.is_open()) {
0565     while (!myReadFile.eof()) {
0566       myReadFile >> ietaEA[ii] >> EA[ii] >> errEA[ii];
0567       if (EA[ii] < 0)
0568         EA[ii] = 0;
0569       ii++;
0570     }
0571   }
0572   myReadFile.close();
0573   std::cout << "Reads " << ii << " points from " << inFile << std::endl;
0574 
0575   gStyle->SetCanvasBorderMode(0);
0576   gStyle->SetCanvasColor(kWhite);
0577   gStyle->SetPadColor(kWhite);
0578   gStyle->SetFillColor(kWhite);
0579   gStyle->SetOptTitle(0);
0580   gStyle->SetOptStat(0);
0581   gStyle->SetOptFit(0);
0582   TFile *f1 = new TFile(rootFile, "RECREATE");
0583   TGraphErrors *eA = new TGraphErrors(ii, ietaEA, EA, errEA, errEA);
0584   eA->SetMarkerStyle(20);
0585   eA->SetMarkerColor(4);
0586   eA->SetLineColor(2);
0587   eA->GetXaxis()->SetTitle("i#eta");
0588   eA->GetXaxis()->SetTitleOffset(0.6);
0589   eA->GetXaxis()->SetTitleSize(0.06);
0590   eA->GetYaxis()->SetTitle("Effective Area");
0591   eA->GetYaxis()->SetTitleOffset(0.6);
0592   eA->GetYaxis()->SetTitleSize(0.06);
0593 
0594   Double_t par[6];
0595   const int nmid = 19;
0596   TF1 *g1 = new TF1("g1", "pol0", 1, nmid);
0597   TF1 *g2 = new TF1("g2", "pol4", nmid, ii);
0598 
0599   eA->Fit(g1, "R");
0600   eA->Fit(g2, "R+");
0601   g1->GetParameters(&par[0]);
0602   g2->GetParameters(&par[1]);
0603 
0604   TCanvas *c2 = new TCanvas("EA vs #eta", "EA vs ieta", 0, 10, 1200, 400);
0605   eA->Draw("AP");
0606   c2->Write();
0607   f1->Close();
0608 
0609   ofstream params;
0610   params.open(outFile);
0611   for (int i = 0; i < 6; i++) {
0612     params << par[i] << std::endl;
0613     std::cout << "Parameter[" << i << "] = " << par[i] << std::endl;
0614   }
0615   params.close();
0616 }
0617 
0618 void FitEovPwithRho(int maxEta, const char *inFile, const char *outFile) {
0619   TFile *file = new TFile(inFile);
0620   std::map<int, TH1D *> histo, histo_uncorr;
0621   char name[100];
0622   for (int ieta = -maxEta; ieta <= maxEta; ieta++) {
0623     sprintf(name, "MPV%d", ieta);
0624     TH1D *h0 = (TH1D *)file->FindObjectAny(name);
0625     histo[ieta] = (h0 != 0) ? (TH1D *)(h0->Clone()) : 0;
0626     sprintf(name, "MPVUn%d", ieta);
0627     TH1D *h1 = (TH1D *)file->FindObjectAny(name);
0628     histo_uncorr[ieta] = (h1 != 0) ? (TH1D *)(h1->Clone()) : 0;
0629   }
0630 
0631   //TFile *f1 =
0632   new TFile(outFile, "RECREATE");
0633   double xlim = maxEta + 0.5;
0634   TH1D *EovPvsieta = new TH1D("Corrected", "Corrected", 2 * maxEta + 1, -xlim, xlim);
0635   TH1D *EovPvsieta_uncorr = new TH1D("Uncorrect", "Uncorrect", 2 * maxEta + 1, -xlim, xlim);
0636 
0637   TF1 *fnc = new TF1("fnc", "gaus");
0638   unsigned int k1(0), k2(0);
0639   for (int ieta = -maxEta; ieta <= maxEta; ieta++) {
0640     if (ieta == 0)
0641       continue;
0642     if (histo[ieta] != 0) {
0643       double mean = histo[ieta]->GetMean();
0644       double rms = histo[ieta]->GetRMS();
0645       TFitResultPtr FitG = histo[ieta]->Fit("fnc", "QRWLS", "", mean - rms, mean + rms);
0646       double a = fnc->GetParameter(1);
0647       double err = FitG->FitResult::Error(1);
0648       histo[ieta]->Write();
0649       ++k1;
0650       int ibin = ieta + maxEta + 1;
0651       EovPvsieta->SetBinContent(ibin, a);
0652       EovPvsieta->SetBinError(ibin, err);
0653       std::cout << "Correct[" << k1 << "] " << ieta << " a " << a << " +- " << err << std::endl;
0654     }
0655 
0656     if (histo_uncorr[ieta] != 0) {
0657       double mean = histo_uncorr[ieta]->GetMean();
0658       double rms = histo_uncorr[ieta]->GetRMS();
0659       TFitResultPtr FitG = histo_uncorr[ieta]->Fit("fnc", "QRWLS", "", mean - rms, mean + rms);
0660       double a = fnc->GetParameter(1);
0661       double err = FitG->FitResult::Error(1);
0662       histo_uncorr[ieta]->Write();
0663       ++k2;
0664       int ibin = ieta + maxEta + 1;
0665       EovPvsieta_uncorr->SetBinContent(ibin, a);
0666       EovPvsieta_uncorr->SetBinError(ibin, err);
0667       std::cout << "Correct[" << k2 << "] " << ieta << " a " << a << " +- " << err << std::endl;
0668     }
0669   }
0670 
0671   gStyle->SetCanvasBorderMode(0);
0672   gStyle->SetCanvasColor(kWhite);
0673   gStyle->SetPadColor(kWhite);
0674   gStyle->SetFillColor(kWhite);
0675   gStyle->SetOptTitle(0);
0676   gStyle->SetOptStat(10);
0677   gStyle->SetOptFit(1);
0678   TCanvas *c3 = new TCanvas("E/P vs ieta", "E/P vs ieta", 0, 10, 1200, 400);
0679   EovPvsieta->GetXaxis()->SetTitle("i#eta");
0680   EovPvsieta->GetYaxis()->SetTitle("MPV[E_{Hcal}/(p_{Track}-E_{Ecal})]");
0681   EovPvsieta->SetMarkerStyle(20);
0682   EovPvsieta->SetMarkerColor(2);
0683   EovPvsieta->SetMarkerSize(1.0);
0684   EovPvsieta->Fit("pol0", "+QRWLS", "", -maxEta, maxEta);
0685 
0686   EovPvsieta_uncorr->SetMarkerStyle(24);
0687   EovPvsieta_uncorr->SetMarkerColor(4);
0688   EovPvsieta_uncorr->SetMarkerSize(1.0);
0689 
0690   EovPvsieta->GetYaxis()->SetRangeUser(0.5, 2.0);
0691   EovPvsieta->Draw();
0692   c3->Update();
0693   TPaveStats *st1 = (TPaveStats *)EovPvsieta->GetListOfFunctions()->FindObject("stats");
0694   if (st1 != nullptr) {
0695     st1->SetY1NDC(0.81);
0696     st1->SetY2NDC(0.90);
0697     st1->SetX1NDC(0.65);
0698     st1->SetX2NDC(0.90);
0699   }
0700   EovPvsieta_uncorr->Draw("sames");
0701   c3->Update();
0702   st1 = (TPaveStats *)EovPvsieta_uncorr->GetListOfFunctions()->FindObject("stats");
0703   std::cout << st1 << std::endl;
0704   if (st1 != nullptr) {
0705     st1->SetY1NDC(0.78);
0706     st1->SetY2NDC(0.81);
0707     st1->SetX1NDC(0.65);
0708     st1->SetX2NDC(0.90);
0709   }
0710   c3->Modified();
0711   c3->Update();
0712   EovPvsieta->Write();
0713   EovPvsieta_uncorr->Write();
0714 }
0715 
0716 void PlotEvsRho(const char *inFile, int eta = 0, int type = 0, bool save = false) {
0717   gStyle->SetCanvasBorderMode(0);
0718   gStyle->SetCanvasColor(kWhite);
0719   gStyle->SetPadColor(kWhite);
0720   gStyle->SetFillColor(kWhite);
0721   gStyle->SetOptTitle(0);
0722   gStyle->SetOptStat(1110);
0723   gStyle->SetOptFit(10);
0724 
0725   TFile *file = new TFile(inFile);
0726   int etamin = (eta == 0) ? ((type == 0) ? 1 : 25) : eta;
0727   int etamax = (eta == 0) ? 25 : eta;
0728   for (int it = etamin; it <= etamax; ++it) {
0729     char name[50];
0730     sprintf(name, "IsovsRho%d", it);
0731     TCanvas *pad;
0732     if (type == 0) {
0733       pad = (TCanvas *)(file->FindObjectAny(name));
0734       pad->Draw();
0735     } else {
0736       sprintf(name, "MPV%d", it);
0737       pad = new TCanvas(name, name, 0, 10, 800, 500);
0738       pad->SetRightMargin(0.10);
0739       pad->SetTopMargin(0.10);
0740       TH1D *h1 = (TH1D *)(file->FindObjectAny(name));
0741       sprintf(name, "MPVUn%d", it);
0742       TH1D *h2 = (TH1D *)(file->FindObjectAny(name));
0743       double ymx1 = h1->GetMaximum();
0744       double ymx2 = h2->GetMaximum();
0745       double ymax = (ymx1 > ymx2) ? ymx1 : ymx2;
0746       h1->GetXaxis()->SetRangeUser(0.5, 2.0);
0747       h2->GetXaxis()->SetRangeUser(0.5, 2.0);
0748       h1->GetYaxis()->SetRangeUser(0, 1.25 * ymax);
0749       h2->GetYaxis()->SetRangeUser(0, 1.25 * ymax);
0750       h1->GetXaxis()->SetTitleSize(0.048);
0751       h1->GetXaxis()->SetTitleOffset(0.8);
0752       h1->GetXaxis()->SetTitle("E_{Hcal}/(p-E_{Ecal})");
0753       h1->GetYaxis()->SetTitleSize(0.048);
0754       h1->GetYaxis()->SetTitleOffset(0.8);
0755       h1->GetYaxis()->SetTitle("Tracks");
0756       h1->SetLineColor(2);
0757       h1->Draw();
0758       pad->Modified();
0759       pad->Update();
0760       TPaveStats *st1 = (TPaveStats *)h1->GetListOfFunctions()->FindObject("stats");
0761       if (st1 != nullptr) {
0762         st1->SetLineColor(2);
0763         st1->SetTextColor(2);
0764         st1->SetY1NDC(0.75);
0765         st1->SetY2NDC(0.90);
0766         st1->SetX1NDC(0.65);
0767         st1->SetX2NDC(0.90);
0768       }
0769       h2->SetLineColor(4);
0770       h2->Draw("sames");
0771       pad->Modified();
0772       pad->Update();
0773       TPaveStats *st2 = (TPaveStats *)h2->GetListOfFunctions()->FindObject("stats");
0774       if (st2 != nullptr) {
0775         st2->SetLineColor(4);
0776         st2->SetTextColor(4);
0777         st2->SetY1NDC(0.60);
0778         st2->SetY2NDC(0.75);
0779         st2->SetX1NDC(0.65);
0780         st2->SetX2NDC(0.90);
0781       }
0782       pad->Modified();
0783       pad->Update();
0784       TF1 *f1 = (TF1 *)h1->GetListOfFunctions()->FindObject("fnc");
0785       if (f1 != nullptr)
0786         f1->SetLineColor(2);
0787       TF1 *f2 = (TF1 *)h2->GetListOfFunctions()->FindObject("fnc");
0788       if (f2 != nullptr)
0789         f2->SetLineColor(4);
0790       pad->Modified();
0791       pad->Update();
0792     }
0793     if (save) {
0794       sprintf(name, "%s.pdf", pad->GetName());
0795       pad->Print(name);
0796     }
0797   }
0798 }
[ 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!