Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-08-12 22:39:59

0001 // STL headers
0002 #include <iostream>
0003 #include <iomanip>
0004 #include <sstream>
0005 
0006 // ROOT headers
0007 #include <TStyle.h>
0008 #include <TCanvas.h>
0009 #include <TTree.h>
0010 #include <TString.h>
0011 #include <TAxis.h>
0012 #include <TProfile.h>
0013 #include <TF1.h>
0014 #include <TH1.h>
0015 #include <TH2F.h>
0016 #include <TGraphErrors.h>
0017 #include <TROOT.h>
0018 #include <TDirectory.h>
0019 #include <TFile.h>
0020 #include <TDirectoryFile.h>
0021 #include <TLegend.h>
0022 #include <TChain.h>
0023 #include <TMath.h>
0024 #include <TLatex.h>
0025 #include <TVirtualFitter.h>
0026 #include <TLorentzVector.h>
0027 #include <TMatrixD.h>
0028 #include <vector>
0029 
0030 // CMSSW headers
0031 #include "Alignment/OfflineValidation/interface/EopElecVariables.h"
0032 
0033 // New enumeration for kind of plots
0034 enum ModeType { TRACK_ETA, TRACK_PHI };
0035 
0036 // New structure for storing information relative to one file
0037 struct HistoType {
0038   TString label;
0039   TFile* file;
0040   TTree* tree;
0041 
0042   std::vector<Double_t> etaRange;
0043   std::vector<Double_t> zRange;
0044   std::vector<Double_t> eRange;
0045 
0046   // Basic histograms
0047   // dimension = numberOfBins x numberOfSteps
0048   UInt_t** selectedTracks;
0049   TH1F*** xAxisBin;
0050   TH1F*** negative;
0051   TH1F*** positive;
0052   TH2F*** Enegative;
0053   TH2F*** Epositive;
0054 
0055   // Intermediate histograms
0056   TH1F** fit;
0057   TH1F** combinedxAxisBin;
0058 
0059   // Final histograms
0060   TGraphErrors* overallGraph;
0061   TH1F* overallhisto;
0062 
0063   // Final fit
0064   TF1* f2;
0065   TString misalignmentfit;
0066 
0067   HistoType() {
0068     label = "";
0069     file = 0;
0070     tree = 0;
0071     overallGraph = 0;
0072     overallhisto = 0;
0073     f2 = 0;
0074     misalignmentfit = "";
0075   }
0076 
0077   ~HistoType() {}
0078 };
0079 
0080 // Prototype of functions
0081 void initializeHistograms(ModeType mode, HistoType& histos);
0082 Bool_t checkArguments(TString variable,  //input
0083                       TString path,
0084                       TString alignmentWithLabel,
0085                       TString outputType,
0086                       Double_t radius,
0087                       Bool_t verbose,
0088                       Double_t givenMin,
0089                       Double_t givenMax,
0090                       ModeType& mode,  // ouput
0091                       std::vector<TFile*>& files,
0092                       std::vector<TString>& labels);
0093 void configureROOTstyle(Bool_t verbose);
0094 Bool_t initializeTree(std::vector<TFile*>& files, std::vector<TTree*>& trees, EopElecVariables* track);
0095 void readTree(ModeType mode, TString variable, HistoType& histo, EopElecVariables* track, Double_t radius);
0096 void fillIntermediateHisto(
0097     ModeType mode, Bool_t verbose, HistoType& histo, TCanvas* ccontrol, TString outputType, Double_t radius);
0098 void fillFinalHisto(ModeType mode, Bool_t verbose, HistoType& histo);
0099 void layoutFinalPlot(ModeType mode,
0100                      std::vector<HistoType>& histos,
0101                      TString outputType,
0102                      TString variable,
0103                      TCanvas* c,
0104                      Double_t givenMin,
0105                      Double_t givenMax);
0106 
0107 std::vector<Double_t> extractgausparams(TH1F* histo, TString option1, TString option2);
0108 
0109 // -----------------------------------------------------------------------------
0110 //
0111 //    Main function : momentumBiasValidation
0112 //
0113 // -----------------------------------------------------------------------------
0114 void momentumBiasValidation(TString variable,
0115                             TString path,
0116                             TString alignmentWithLabel,
0117                             TString outputType,
0118                             Double_t radius = 1.,
0119                             Bool_t verbose = false,
0120                             Double_t givenMin = 0.,
0121                             Double_t givenMax = 0.) {
0122   // Displaying first message
0123   std::cout << "!!! Welcome to momentumBiasValidation !!!" << std::endl;
0124   std::cout << std::endl;
0125   time_t start = time(0);
0126 
0127   // Checking and decoding the arguments
0128   std::cout << "Checking arguments ..." << std::endl;
0129   ModeType mode;                // variable to study
0130   std::vector<TFile*> files;    // list of input files
0131   std::vector<TString> labels;  // list of input labels
0132   if (!checkArguments(
0133           variable, path, alignmentWithLabel, outputType, radius, verbose, givenMin, givenMax, mode, files, labels))
0134     return;
0135   else {
0136     std::cout << "-> Number of files: " << files.size() << std::endl;
0137   }
0138 
0139   // Initializing the TTree
0140   std::cout << "Initializing the TTree ..." << std::endl;
0141   std::vector<TTree*> trees(files.size(), 0);
0142   EopElecVariables* track = new EopElecVariables();
0143   if (!initializeTree(files, trees, track))
0144     delete track;
0145   return;
0146 
0147   // Configuring ROOT style
0148   std::cout << "Configuring the ROOT style ..." << std::endl;
0149   configureROOTstyle(verbose);
0150 
0151   TCanvas* c = new TCanvas("c", "Canvas", 0, 0, 1150, 800);
0152   TCanvas* ccontrol = new TCanvas("ccontrol", "controlCanvas", 0, 0, 600, 300);
0153 
0154   // Creating histos
0155   std::vector<HistoType> histos(files.size());
0156 
0157   // Loop over files
0158 
0159   // To save the table cut
0160   TFile* histo1;
0161   std::string fileName;
0162 
0163   for (unsigned int ifile = 0; ifile < histos.size(); ifile++) {
0164     HistoType& histo = histos[ifile];
0165 
0166     // setting labels
0167     histo.label = labels[ifile];
0168     histo.file = files[ifile];
0169     histo.tree = trees[ifile];
0170 
0171     fileName = histo.file->GetName();
0172     fileName = fileName.substr(path.Sizeof() - 1);
0173     fileName = "Plots" + fileName;
0174 
0175     std::cout << "NAME :" << fileName << std::endl;
0176     histo1 = new TFile(fileName.c_str(), "RECREATE");
0177 
0178     // Getting and saving cut flow values from tree producer step
0179     TH1D* Nevents;
0180     TH1D* NeventsTriggered;
0181     TH1D* Nevents2elec;
0182     TH1D* Ntracks;
0183     TH1D* NtracksFiltered;
0184     TH1D* NtracksFirstPtcut;
0185     TH1D* NtracksOneSC;
0186 
0187     Nevents = dynamic_cast<TH1D*>(histo.file->Get("energyOverMomentumTree/nEvents"));
0188     NeventsTriggered = dynamic_cast<TH1D*>(histo.file->Get("energyOverMomentumTree/nEventsTriggered"));
0189     Nevents2elec = dynamic_cast<TH1D*>(histo.file->Get("energyOverMomentumTree/nEvents2Elec"));
0190     Ntracks = dynamic_cast<TH1D*>(histo.file->Get("energyOverMomentumTree/nTracks"));
0191     NtracksFiltered = dynamic_cast<TH1D*>(histo.file->Get("energyOverMomentumTree/nTracksFiltered"));
0192     NtracksFirstPtcut = dynamic_cast<TH1D*>(histo.file->Get("energyOverMomentumTree/cut_Ptmin"));
0193     NtracksOneSC = dynamic_cast<TH1D*>(histo.file->Get("energyOverMomentumTree/cut_OneSCmatch"));
0194 
0195     Nevents->Write();
0196     NeventsTriggered->Write();
0197     Nevents2elec->Write();
0198     Ntracks->Write();
0199     NtracksFiltered->Write();
0200     NtracksFirstPtcut->Write();
0201     NtracksOneSC->Write();
0202 
0203     // Initializing the range in Energy
0204     histo.eRange.clear();
0205     histo.eRange.push_back(30);
0206     histo.eRange.push_back(40);
0207     histo.eRange.push_back(60);
0208     //histo.eRange.push_back(100);
0209 
0210     // Dump at screen range Energy
0211     if (ifile == 0) {
0212       std::cout << "Range used for energy : ";
0213       for (unsigned int i = 0; i < (histo.eRange.size() - 1); i++)
0214         std::cout << "[" << histo.eRange[i] << "]-[" << histo.eRange[i + 1] << "] ";
0215       std::cout << std::endl;
0216     }
0217 
0218     // Initializing the range in Phi
0219     if (mode == TRACK_ETA) {
0220       Double_t begin = -0.9;  //-1.40;//-1.65;
0221       Double_t end = 0.9;     //1.40;//1.65;
0222       //      UInt_t   nsteps  = 18;
0223       UInt_t nsteps = 11;  //8;
0224       Double_t binsize = (end - begin) / static_cast<Float_t>(nsteps);
0225 
0226       histo.etaRange.resize(nsteps + 1, 0.);
0227       histo.zRange.resize(nsteps + 1, 0.);
0228       for (UInt_t i = 0; i < histo.etaRange.size(); i++) {
0229         histo.etaRange[i] = begin + i * binsize;
0230         histo.zRange[i] = radius * 100. / TMath::Tan(2 * TMath::ATan(TMath::Exp(-(begin + i * binsize))));
0231       }
0232     }
0233 
0234     // Initializing the range in Eta
0235     else if (mode == TRACK_PHI) {
0236       Double_t begin = -TMath::Pi();
0237       Double_t end = +TMath::Pi();
0238       UInt_t nsteps = 8;  //
0239       Double_t binsize = (end - begin) / static_cast<Float_t>(nsteps);
0240 
0241       histo.etaRange.resize(nsteps + 1, 0.);
0242       for (UInt_t i = 0; i < histo.etaRange.size(); i++) {
0243         histo.etaRange[i] = begin + i * binsize;
0244       }
0245     }
0246 
0247     // Dump at screen the range in Eta (or Phi)
0248     if (ifile == 0) {
0249       std::cout << "Range used for ";
0250       if (mode == TRACK_ETA)
0251         std::cout << "eta";
0252       else
0253         std::cout << "phi";
0254       std::cout << " : " << std::endl;
0255       ;
0256       for (UInt_t i = 0; i < (histo.etaRange.size() - 1); i++)
0257         std::cout << "[" << histo.etaRange[i] << "]-[" << histo.etaRange[i + 1] << "] ";
0258       std::cout << std::endl;
0259     }
0260 
0261     // Initializing histos
0262     std::cout << "Initialzing histograms ..." << std::endl;
0263     initializeHistograms(mode, histo);
0264 
0265     // Filling with events
0266     std::cout << "Reading the TFile ..." << std::endl;
0267     readTree(mode, variable, histo, track, radius);
0268 
0269     //Filling histograms
0270     std::cout << "Filling histograms ..." << std::endl;
0271     fillIntermediateHisto(mode, verbose, histo, ccontrol, outputType, radius);
0272     fillFinalHisto(mode, verbose, histo);
0273   }
0274 
0275   // Achieving the final plot
0276   std::cout << "Final plot ..." << std::endl;
0277   layoutFinalPlot(mode, histos, outputType, variable, c, givenMin, givenMax);
0278 
0279   // Displaying final message
0280   time_t end = time(0);
0281   std::cout << "Done in " << static_cast<int>(difftime(end, start)) / 60 << " min and "
0282             << static_cast<int>(difftime(end, start)) % 60 << " sec." << std::endl;
0283 
0284   delete track;
0285   delete histo1;
0286 }
0287 
0288 // -----------------------------------------------------------------------------
0289 //
0290 //    Auxiliary function : fillIntermediateHisto
0291 //
0292 // -----------------------------------------------------------------------------
0293 void fillIntermediateHisto(
0294     ModeType mode, Bool_t verbose, HistoType& histo, TCanvas* ccontrol, TString outputType, Double_t radius) {
0295   // Loop over eta or phi range
0296   for (UInt_t j = 0; j < (histo.etaRange.size() - 1); j++) {
0297     // Display eta (or phi) range
0298     if (verbose) {
0299       std::cout << histo.etaRange[j] << " < ";
0300       if (mode == TRACK_ETA)
0301         std::cout << "eta";
0302       else
0303         std::cout << "phi";
0304       std::cout << " < " << histo.etaRange[j + 1] << std::endl;
0305     }
0306 
0307     // Loop over energy range
0308     for (UInt_t i = 0; i < (histo.eRange.size() - 1); i++) {
0309       // Verbose mode : saving histo positive and negative
0310       TString controlName;
0311       if (verbose) {
0312         // filename for control plots
0313         controlName = "controlEOP/control_eop_";
0314         controlName += histo.label;
0315         controlName += "_bin";
0316         controlName += j;
0317         controlName += "_energy";
0318         controlName += i;
0319         controlName += ".";
0320         controlName += outputType;
0321         ccontrol->cd();
0322 
0323         Double_t posMax = histo.positive[i][j]->GetMaximum();
0324         Double_t negMax = histo.negative[i][j]->GetMaximum();
0325         if (posMax > negMax)
0326           histo.negative[i][j]->SetMaximum(1.1 * posMax);
0327 
0328         histo.negative[i][j]->DrawClone();
0329         histo.positive[i][j]->DrawClone("same");
0330       }
0331 
0332       // Fitting by a gaussian and extracting fit parameters
0333       std::vector<Double_t> curvNeg = extractgausparams(histo.negative[i][j], "", "");
0334       std::vector<Double_t> curvPos = extractgausparams(histo.positive[i][j], "", "same");
0335 
0336       // Verbose mode : saving gaussian fit plots
0337       if (verbose) {
0338         ccontrol->Print(controlName);
0339       }
0340 
0341       // Initial misalignment value
0342       Double_t misalignment = 0.;
0343       Double_t misaliUncert = 1000.;
0344 
0345       // Compute misalignment only if there are selected tracks
0346       if (histo.selectedTracks[i][j] != 0) {
0347         // Setting gaussian range
0348         histo.Enegative[i][j]->GetYaxis()->SetRangeUser(curvNeg[2], curvNeg[3]);
0349         histo.Epositive[i][j]->GetYaxis()->SetRangeUser(curvPos[2], curvPos[3]);
0350 
0351         // Getting mean values from histogram
0352         Double_t meanEnergyNeg = histo.Enegative[i][j]->GetMean();
0353         Double_t meanEnergyPos = histo.Epositive[i][j]->GetMean();
0354         Double_t meanEnergy = (meanEnergyNeg + meanEnergyPos) /
0355                               2.;  // use mean of positive and negative tracks to reduce energy dependence
0356 
0357         // Verbose mode : displaying difference between positive and negative means
0358         if (verbose) {
0359           std::cout << "difference in energy between positive and negative tracks: " << meanEnergyNeg - meanEnergyPos
0360                     << std::endl;
0361         }
0362 
0363         // Compute misalignment
0364         if (mode == TRACK_ETA) {
0365           misalignment = 1000000. * 0.5 *
0366                          (-TMath::ASin((0.57 * radius / meanEnergy) * curvNeg[0]) +
0367                           TMath::ASin((0.57 * radius / meanEnergy) * curvPos[0]));
0368           misaliUncert =
0369               1000000. * 0.5 *
0370               (TMath::Sqrt((0.57 * 0.57 * radius * radius * curvNeg[1] * curvNeg[1]) /
0371                                (meanEnergy * meanEnergy - 0.57 * 0.57 * radius * radius * curvNeg[0] * curvNeg[0]) +
0372                            (0.57 * 0.57 * radius * radius * curvPos[1] * curvPos[1]) /
0373                                (meanEnergy * meanEnergy - 0.57 * 0.57 * radius * radius * curvPos[0] * curvPos[0])));
0374         } else if (mode == TRACK_PHI) {
0375           misalignment = 1000. * (curvPos[0] - curvNeg[0]) / (curvPos[0] + curvNeg[0]);
0376           misaliUncert = 1000. * 2 / ((curvPos[0] + curvNeg[0]) * (curvPos[0] + curvNeg[0])) *
0377                          TMath::Sqrt((curvPos[0] * curvPos[0] * curvPos[1] * curvPos[1]) +
0378                                      (curvNeg[0] * curvNeg[0] * curvNeg[1] * curvNeg[1]));
0379         }
0380       }
0381 
0382       // Verbose mode : displaying computed misalignment
0383       if (verbose)
0384         std::cout << "misalignment: " << misalignment << "+-" << misaliUncert << std::endl << std::endl;
0385 
0386       // Fill intermediate histogram : histo.fit
0387       histo.fit[j]->SetBinContent(i + 1, misalignment);
0388       histo.fit[j]->SetBinError(i + 1, misaliUncert);
0389 
0390       // Fill intermediate histogram : histo.combinedxAxisBin
0391       Double_t xBinCentre = histo.xAxisBin[i][j]->GetMean();
0392       Double_t xBinCenUnc = histo.xAxisBin[i][j]->GetMeanError();
0393       histo.combinedxAxisBin[j]->SetBinContent(i + 1, xBinCentre);
0394       histo.combinedxAxisBin[j]->SetBinError(i + 1, xBinCenUnc);
0395     }
0396   }
0397 }
0398 
0399 // -----------------------------------------------------------------------------
0400 //
0401 //    Auxiliary function : fillFinalHisto
0402 //
0403 // -----------------------------------------------------------------------------
0404 void fillFinalHisto(ModeType mode, Bool_t verbose, HistoType& histo) {
0405   TString fitOption;
0406   if (verbose)
0407     fitOption = "";
0408   else
0409     fitOption = "Q";
0410   for (UInt_t i = 0; i < (histo.etaRange.size() - 1); i++) {
0411     Double_t overallmisalignment;
0412     Double_t overallmisaliUncert;
0413     Double_t overallxBin;
0414     Double_t overallxBinUncert;
0415 
0416     // calculate mean of different energy bins
0417     TF1* fit = new TF1("fit", "pol0", 0, histo.eRange.size() - 1);
0418     TF1* fit2 = new TF1("fit2", "pol0", 0, histo.eRange.size() - 1);
0419     histo.fit[i]->Fit("fit", fitOption + "0");
0420     histo.combinedxAxisBin[i]->Fit("fit2", "Q0");
0421     overallmisalignment = fit->GetParameter(0);
0422     overallmisaliUncert = fit->GetParError(0);
0423     overallxBin = fit2->GetParameter(0);
0424     overallxBinUncert = fit2->GetParError(0);
0425     fit->Delete();
0426     fit2->Delete();
0427 
0428     // Fill final histograms
0429     histo.overallhisto->SetBinContent(i + 1, overallmisalignment);
0430     histo.overallhisto->SetBinError(i + 1, overallmisaliUncert);
0431     histo.overallGraph->SetPoint(i, overallxBin, overallmisalignment);
0432     histo.overallGraph->SetPointError(i, overallxBinUncert, overallmisaliUncert);
0433   }
0434 
0435   // Fit to final histogram
0436   TString func = "func";
0437   func += histo.label;
0438   if (mode == TRACK_ETA)
0439     histo.f2 = new TF1(func, "[0]+[1]*x/100.", -500, 500);  //Divide by 100. cm->m
0440   if (mode == TRACK_PHI)
0441     histo.f2 = new TF1(func, "[0]+[1]*TMath::Cos(x+[2])", -500, 500);
0442 
0443   // Fitting final histogram
0444   histo.overallGraph->Fit(func, fitOption + "mR0+");
0445 
0446   // Verbose mode : displaying covariance from fit
0447   if (verbose) {
0448     std::cout << "Covariance Matrix:" << std::endl;
0449     TVirtualFitter* fitter = TVirtualFitter::GetFitter();
0450     TMatrixD matrix(2, 2, fitter->GetCovarianceMatrix());
0451     Double_t oneOne = fitter->GetCovarianceMatrixElement(0, 0);
0452     Double_t oneTwo = fitter->GetCovarianceMatrixElement(0, 1);
0453     Double_t twoOne = fitter->GetCovarianceMatrixElement(1, 0);
0454     Double_t twoTwo = fitter->GetCovarianceMatrixElement(1, 1);
0455 
0456     std::cout << "( " << oneOne << ", " << twoOne << ")" << std::endl;
0457     std::cout << "( " << oneTwo << ", " << twoTwo << ")" << std::endl;
0458   }
0459 
0460   // Displaying at screen  fit parameters
0461   if (mode == TRACK_ETA) {
0462     std::cout << "const: " << histo.f2->GetParameter(0) << "+-" << histo.f2->GetParError(0)
0463               << ", slope: " << histo.f2->GetParameter(1) << "+-" << histo.f2->GetParError(1) << std::endl;
0464   } else if (mode == TRACK_PHI) {
0465     std::cout << "const: " << histo.f2->GetParameter(0) << "+-" << histo.f2->GetParError(0)
0466               << ", amplitude: " << histo.f2->GetParameter(1) << "+-" << histo.f2->GetParError(1)
0467               << ", shift: " << histo.f2->GetParameter(2) << "+-" << histo.f2->GetParError(2) << std::endl;
0468   }
0469   std::cout << "fit probability: " << histo.f2->GetProb() << std::endl;
0470   std::cout << "fit chi2       : " << histo.f2->GetChisquare() << std::endl;
0471   std::cout << "fit chi2/Ndof  : " << histo.f2->GetChisquare() / static_cast<Double_t>(histo.f2->GetNDF()) << std::endl;
0472 
0473   // Adding the fit function for the legend
0474   Char_t misalignmentfitchar[20];
0475   sprintf(misalignmentfitchar, "%1.f", histo.f2->GetParameter(0));
0476   histo.misalignmentfit += "(";
0477   histo.misalignmentfit += misalignmentfitchar;
0478   histo.misalignmentfit += "#pm";
0479   sprintf(misalignmentfitchar, "%1.f", histo.f2->GetParError(0));
0480   histo.misalignmentfit += misalignmentfitchar;
0481   if (mode == TRACK_ETA) {
0482     histo.misalignmentfit += ")#murad #upoint r[m] + (";
0483     sprintf(misalignmentfitchar, "%1.f", histo.f2->GetParameter(1));
0484     histo.misalignmentfit += misalignmentfitchar;
0485     histo.misalignmentfit += "#pm";
0486     sprintf(misalignmentfitchar, "%1.f", histo.f2->GetParError(1));
0487     histo.misalignmentfit += misalignmentfitchar;
0488     histo.misalignmentfit += ")#murad #upoint z[m]";
0489   } else if (mode == TRACK_PHI) {
0490     histo.misalignmentfit += ") + (";
0491     sprintf(misalignmentfitchar, "%1.f", histo.f2->GetParameter(1));
0492     histo.misalignmentfit += misalignmentfitchar;
0493     histo.misalignmentfit += "#pm";
0494     sprintf(misalignmentfitchar, "%1.f", histo.f2->GetParError(1));
0495     histo.misalignmentfit += misalignmentfitchar;
0496     histo.misalignmentfit += ") #upoint cos(#phi";
0497     if (histo.f2->GetParameter(2) > 0.)
0498       histo.misalignmentfit += "+";
0499     sprintf(misalignmentfitchar, "%1.1f", histo.f2->GetParameter(2));
0500     histo.misalignmentfit += misalignmentfitchar;
0501     histo.misalignmentfit += "#pm";
0502     sprintf(misalignmentfitchar, "%1.1f", histo.f2->GetParError(2));
0503     histo.misalignmentfit += misalignmentfitchar;
0504     histo.misalignmentfit += ")";
0505   }
0506 }
0507 
0508 // -----------------------------------------------------------------------------
0509 //
0510 //    Auxiliary function : LayoutFinalPlot
0511 //
0512 // -----------------------------------------------------------------------------
0513 void layoutFinalPlot(ModeType mode,
0514                      std::vector<HistoType>& histos,
0515                      TString outputType,
0516                      TString variable,
0517                      TCanvas* c,
0518                      Double_t givenMin,
0519                      Double_t givenMax) {
0520   // Create a legend
0521   TLegend* leg = new TLegend(0.13, 0.78, 0.98, 0.98);
0522   leg->SetFillColor(10);
0523   leg->SetTextFont(42);
0524   leg->SetTextSize(0.038);
0525   if (variable == "phi1")
0526     leg->SetHeader("low #eta tracks (#eta < -0.9)");
0527   if (variable == "phi2")
0528     leg->SetHeader("central tracks (|#eta| < 0.9)");
0529   if (variable == "phi3")
0530     leg->SetHeader("high #eta tracks (#eta > 0.9)");
0531 
0532   // Setting display options to final histos
0533   for (UInt_t i = 0; i < histos.size(); i++) {
0534     // Setting axis title to final histos
0535     if (mode == TRACK_ETA) {
0536       histos[i].overallhisto->GetXaxis()->SetTitle("z [cm]");
0537       histos[i].overallhisto->GetYaxis()->SetTitle("misalignment #Delta#phi[#murad]");
0538     }
0539     if (mode == TRACK_PHI) {
0540       histos[i].overallhisto->GetXaxis()->SetTitle("#phi");
0541       histos[i].overallhisto->GetYaxis()->SetTitle("misalignment #Delta#phi[a.u.]");
0542     }
0543 
0544     // Fit function
0545     histos[i].f2->SetLineColor(i + 1);
0546     histos[i].f2->SetLineStyle(i + 1);
0547     histos[i].f2->SetLineWidth(2);
0548 
0549     // Other settings
0550     histos[i].overallhisto->GetYaxis()->SetTitleOffset(1.05);
0551     histos[i].overallhisto->GetYaxis()->SetTitleSize(0.065);
0552     histos[i].overallhisto->GetYaxis()->SetLabelSize(0.065);
0553     histos[i].overallhisto->GetXaxis()->SetTitleOffset(0.8);
0554     histos[i].overallhisto->GetXaxis()->SetTitleSize(0.065);
0555     histos[i].overallhisto->GetXaxis()->SetLabelSize(0.065);
0556     histos[i].overallhisto->SetLineWidth(2);
0557     histos[i].overallhisto->SetLineColor(i + 1);
0558     histos[i].overallhisto->SetMarkerColor(i + 1);
0559     histos[i].overallGraph->SetLineWidth(2);
0560     histos[i].overallGraph->SetLineColor(i + 1);
0561     histos[i].overallGraph->SetMarkerColor(i + 1);
0562     histos[i].overallGraph->SetMarkerStyle(i + 20);
0563     histos[i].overallGraph->SetMarkerSize(2);
0564   }
0565 
0566   // set pad margins
0567   c->cd();
0568   gPad->SetTopMargin(0.02);
0569   gPad->SetBottomMargin(0.12);
0570   gPad->SetLeftMargin(0.13);
0571   gPad->SetRightMargin(0.02);
0572 
0573   // Determining common reasonable y-axis range
0574   Double_t overallmax = 0.;
0575   Double_t overallmin = 0.;
0576   for (UInt_t i = 0; i < histos.size(); i++) {
0577     // Getting maximum from overallmax
0578     overallmax = TMath::Max(overallmax,
0579                             histos[i].overallhisto->GetMaximum() +
0580                                 histos[i].overallhisto->GetBinError(histos[i].overallhisto->GetMaximumBin()) +
0581                                 0.55 * (histos[i].overallhisto->GetMaximum() +
0582                                         histos[i].overallhisto->GetBinError(histos[i].overallhisto->GetMaximumBin()) -
0583                                         histos[i].overallhisto->GetMinimum() +
0584                                         histos[i].overallhisto->GetBinError(histos[i].overallhisto->GetMinimumBin())));
0585 
0586     // Getting minimum from overallmin
0587     overallmin = TMath::Min(overallmin,
0588                             histos[i].overallhisto->GetMinimum() -
0589                                 fabs(histos[i].overallhisto->GetBinError(histos[i].overallhisto->GetMinimumBin())) -
0590                                 0.1 * (histos[i].overallhisto->GetMaximum() +
0591                                        histos[i].overallhisto->GetBinError(histos[i].overallhisto->GetMaximumBin()) -
0592                                        histos[i].overallhisto->GetMinimum() +
0593                                        histos[i].overallhisto->GetBinError(histos[i].overallhisto->GetMinimumBin())));
0594   }
0595 
0596   // Applying common y-axis range to each final histo
0597   for (UInt_t i = 0; i < histos.size(); i++) {
0598     histos[i].overallhisto->SetMaximum(overallmax);
0599     histos[i].overallhisto->SetMinimum(overallmin);
0600     if (givenMax != 0)
0601       histos[i].overallhisto->SetMaximum(givenMax);
0602     if (givenMin != 0)
0603       histos[i].overallhisto->SetMinimum(givenMin);
0604     histos[i].overallhisto->DrawClone("axis");
0605 
0606     // set histogram errors to a small value as only errors of the graph should be shown
0607     for (Int_t j = 0; j < histos[i].overallhisto->GetNbinsX(); j++)
0608       histos[i].overallhisto->SetBinError(j + 1, 0.00001);
0609 
0610     // draw final histogram
0611     histos[i].overallhisto->DrawClone("pe1 same");
0612     histos[i].f2->DrawClone("same");
0613     histos[i].overallGraph->DrawClone("|| same");
0614     histos[i].overallGraph->DrawClone("pz same");
0615     histos[i].overallGraph->SetLineStyle(i + 1);
0616     leg->AddEntry(histos[i].overallGraph, histos[i].label + " (" + histos[i].misalignmentfit + ")", "Lp");
0617   }
0618 
0619   leg->Draw();
0620 
0621   // Saving plots
0622   if (variable == "eta")
0623     c->Print("twist_validation." + outputType);
0624   else if (variable == "phi")
0625     c->Print("sagitta_validation_all." + outputType);
0626   else if (variable == "phi1")
0627     c->Print("sagitta_validation_lowEta." + outputType);
0628   else if (variable == "phi2")
0629     c->Print("sagitta_validation_centralEta." + outputType);
0630   else if (variable == "phi3")
0631     c->Print("sagitta_validation_highEta." + outputType);
0632 
0633   delete leg;
0634 }
0635 
0636 // -----------------------------------------------------------------------------
0637 //
0638 //    Auxiliary function : initialize Histograms
0639 //
0640 // -----------------------------------------------------------------------------
0641 void initializeHistograms(ModeType mode, HistoType& histo) {
0642   // -----------------------------------------
0643   //      Initializing Basic Histograms
0644   // -----------------------------------------
0645 
0646   // Allocated memory
0647   histo.selectedTracks = new UInt_t*[histo.eRange.size() - 1];
0648   histo.xAxisBin = new TH1F**[histo.eRange.size() - 1];
0649   histo.negative = new TH1F**[histo.eRange.size() - 1];
0650   histo.positive = new TH1F**[histo.eRange.size() - 1];
0651   histo.Enegative = new TH2F**[histo.eRange.size() - 1];
0652   histo.Epositive = new TH2F**[histo.eRange.size() - 1];
0653   for (unsigned int i = 0; i < (histo.eRange.size() - 1); i++) {
0654     histo.selectedTracks[i] = new UInt_t[histo.etaRange.size() - 1];
0655     histo.xAxisBin[i] = new TH1F*[histo.etaRange.size() - 1];
0656     histo.negative[i] = new TH1F*[histo.etaRange.size() - 1];
0657     histo.positive[i] = new TH1F*[histo.etaRange.size() - 1];
0658     histo.Enegative[i] = new TH2F*[histo.etaRange.size() - 1];
0659     histo.Epositive[i] = new TH2F*[histo.etaRange.size() - 1];
0660   }
0661 
0662   // Loop over energy range
0663   unsigned int index = 0;
0664   for (unsigned int i = 0; i < (histo.eRange.size() - 1); i++) {
0665     // Labels for histo title
0666     Char_t tmpstep[5];
0667     sprintf(tmpstep, "%1.1fGeV", histo.eRange[i]);
0668     TString lowEnergyBorder = tmpstep;
0669     sprintf(tmpstep, "%1.1fGeV", histo.eRange[i + 1]);
0670     TString highEnergyBorder = tmpstep;
0671     TString eTitle = lowEnergyBorder + " #leq ";
0672     if (mode == TRACK_ETA)
0673       eTitle += "E";
0674     else
0675       eTitle += "E_{T}";
0676     eTitle += " < " + highEnergyBorder;
0677 
0678     for (unsigned int j = 0; j < (histo.etaRange.size() - 1); j++) {
0679       // Labels for histo name
0680       Char_t tmpstep[5];
0681       sprintf(tmpstep, "%1.1fGeV", histo.eRange[i]);
0682       TString lowEnergyBorder = tmpstep;
0683       sprintf(tmpstep, "%1.1fGeV", histo.eRange[i + 1]);
0684       TString highEnergyBorder = tmpstep;
0685       TString eTitle = lowEnergyBorder + " #leq ";
0686       if (mode == TRACK_ETA)
0687         eTitle += "E";
0688       else
0689         eTitle += "E_{T}";
0690       eTitle += " < " + highEnergyBorder;
0691 
0692       index++;
0693       Char_t histochar[20];
0694       sprintf(histochar, "%i", index);
0695       TString histotrg = histochar;
0696       histotrg += " (" + histo.label + ")";
0697 
0698       // Creating histograms
0699       histo.negative[i][j] = new TH1F("negative" + histotrg, "negative (" + eTitle + ")", 50, 0, 2);
0700       histo.positive[i][j] = new TH1F("positive" + histotrg, "positive (" + eTitle + ")", 50, 0, 2);
0701       histo.Enegative[i][j] = new TH2F("Enegative" + histotrg, "Enegative (" + eTitle + ")", 5000, 0, 500, 50, 0, 2);
0702       histo.Epositive[i][j] = new TH2F("Epositive" + histotrg, "Epositive (" + eTitle + ")", 5000, 0, 500, 50, 0, 2);
0703       histo.xAxisBin[i][j] = new TH1F("xAxisBin" + histotrg, "xAxisBin (" + eTitle + ")", 1000, -500, 500);
0704 
0705       // Sum of squares of weight for each histogram
0706       histo.Enegative[i][j]->Sumw2();
0707       histo.Epositive[i][j]->Sumw2();
0708       histo.xAxisBin[i][j]->Sumw2();
0709       histo.negative[i][j]->Sumw2();
0710       histo.positive[i][j]->Sumw2();
0711 
0712       //SetDirectory(0)
0713       histo.Enegative[i][j]->SetDirectory(0);
0714       histo.Epositive[i][j]->SetDirectory(0);
0715       histo.xAxisBin[i][j]->SetDirectory(0);
0716       histo.negative[i][j]->SetDirectory(0);
0717       histo.positive[i][j]->SetDirectory(0);
0718 
0719       // Set color
0720       histo.negative[i][j]->SetLineColor(kGreen);
0721       histo.positive[i][j]->SetLineColor(kRed);
0722     }
0723   }
0724 
0725   // -----------------------------------------
0726   //   Initializing Intermediate Histograms
0727   // -----------------------------------------
0728   histo.fit = new TH1F*[histo.etaRange.size() - 1];
0729   histo.combinedxAxisBin = new TH1F*[histo.etaRange.size() - 1];
0730 
0731   for (UInt_t j = 0; j < (histo.etaRange.size() - 1); j++) {
0732     Char_t histochar[20];
0733     sprintf(histochar, "%i", j + 1);
0734     TString histotrg = histochar;
0735     histotrg += "(" + histo.label + ")";
0736     histo.fit[j] = new TH1F("fithisto" + histotrg,       //name
0737                             "fithisto" + histotrg,       //title
0738                             (histo.eRange.size() - 1),   //nbins
0739                             0.,                          //min
0740                             (histo.eRange.size() - 1));  //max
0741     histo.fit[j]->SetDirectory(0);
0742     histo.combinedxAxisBin[j] = new TH1F("combinedxAxisBin" + histotrg,  //name
0743                                          "combinedxAxisBin" + histotrg,  //title
0744                                          (histo.eRange.size() - 1),      //nbins
0745                                          0.,                             //min
0746                                          (histo.eRange.size() - 1));     //max
0747     histo.combinedxAxisBin[j]->SetDirectory(0);
0748   }
0749 
0750   // -----------------------------------------
0751   //       Initializing Final Histograms
0752   // -----------------------------------------
0753   if (mode == TRACK_ETA) {
0754     histo.overallhisto = new TH1F(
0755         "overallhisto (" + histo.label + ")", "overallhisto", (histo.zRange.size() - 1), &histo.zRange.front());
0756   } else {
0757     histo.overallhisto = new TH1F("overallhisto (" + histo.label + ")",
0758                                   "overallhisto",
0759                                   (histo.etaRange.size() - 1),
0760                                   histo.etaRange[0],
0761                                   histo.etaRange[histo.etaRange.size() - 1]);
0762   }
0763   histo.overallhisto->SetDirectory(0);
0764 
0765   histo.overallGraph = new TGraphErrors(histo.overallhisto);
0766 }
0767 
0768 // -----------------------------------------------------------------------------
0769 //
0770 //    Auxiliary function : readTree
0771 //
0772 // -----------------------------------------------------------------------------
0773 void readTree(ModeType mode, TString variable, HistoType& histo, EopElecVariables* track, Double_t radius) {
0774   // Displaying sample name
0775   Long64_t nevent = (Long64_t)histo.tree->GetEntries();
0776   std::cout << "Reading sample labeled by '" << histo.label << "' containing " << nevent << " events ..." << std::endl;
0777 
0778   // Reseting counters
0779   TH1D* NtracksMacro = new TH1D("NtracksMacro", "NtracksMacro", 1, 0, 1);
0780   TH1D* usedtracks = new TH1D("usedtracks", "usedtracks", 1, 0, 1);
0781   TH1D* cut_Mz = new TH1D("cut_Mz", "cut_Mz", 1, 0, 1);
0782   TH1D* cut_eta = new TH1D("cut_eta", "cut_eta", 1, 0, 1);
0783   TH1D* cut_ChargedIso = new TH1D("cut_ChargedIso", "cut_ChargedIso", 1, 0, 1);
0784   TH1D* cut_Ecal = new TH1D("cut_Ecal", "cut_Ecal", 1, 0, 1);
0785   TH1D* cut_HoverE = new TH1D("cut_HoverE", "cut_HoverE", 1, 0, 1);
0786   TH1D* cut_NeutralIso = new TH1D("cut_NeutralIso", "cut_NeutralIso", 1, 0, 1);
0787   TH1D* cut_nHits = new TH1D("cut_nHits", "cut_nHits", 1, 0, 1);
0788   TH1D* cut_nLostHits = new TH1D("cut_nLostHits", "cut_nLostHits", 1, 0, 1);
0789   TH1D* cut_outerRadius = new TH1D("cut_outerRadius", "cut_outerRadius", 1, 0, 1);
0790   TH1D* cut_normalizedChi2 = new TH1D("cut_normalizedChi2", "cut_normalizedChi2", 1, 0, 1);
0791   TH1D* cut_fbrem = new TH1D("cut_fbrem", "cut_fbrem", 1, 0, 1);
0792   TH1D* cut_nCluster = new TH1D("cut_nCluster", "cut_nCluster", 1, 0, 1);
0793   TH1D* cut_EcalRange = new TH1D("cut_EcalRange", "cut_EcalRange", 1, 0, 1);
0794   TH1D* cut_trigger = new TH1D("cut_trigger", "cut_trigger", 1, 0, 1);
0795 
0796   TH1D* P = new TH1D("P", "P", 180, 0, 180);
0797   TH1D* eta = new TH1D("eta", "eta", 100, -5., 5.);
0798   TH1D* Pt = new TH1D("Pt", "Pt", 1000, 0, 1000);
0799   TH1D* Ecal = new TH1D("Ecal", "Ecal", 100, 0, 180);
0800   TH1D* HcalEnergyIn01 = new TH1D("HcalEnergyIn01", "HcalEnergyIn01", 100, 0, 40);
0801   TH1D* HcalEnergyIn02 = new TH1D("HcalEnergyIn02", "HcalEnergyIn02", 100, 0, 40);
0802   TH1D* HcalEnergyIn03 = new TH1D("HcalEnergyIn03", "HcalEnergyIn03", 100, 0, 40);
0803   TH1D* HcalEnergyIn04 = new TH1D("HcalEnergyIn04", "HcalEnergyIn04", 100, 0, 40);
0804   TH1D* HcalEnergyIn05 = new TH1D("HcalEnergyIn05", "HcalEnergyIn05", 100, 0, 40);
0805   TH1D* SumPt = new TH1D("SumPt", "SumPt", 100, 0, 2);
0806   TH1D* HoverE = new TH1D("HoverE", "HoverE", 50, 0, 0.6);
0807   TH1D* distTo1stSC = new TH1D("distTo1stSC", "distTo1stSC", 50, 0, 0.1);
0808   TH1D* distTo2ndSC = new TH1D("distTo2ndSC", "distTo2ndSC", 50, 0, 1.5);
0809   TH1D* fbrem = new TH1D("fbrem", "fbrem", 50, -0.2, 1.);
0810   TH1D* nBasicClus = new TH1D("nBasicClus", "nBasicClus", 12, 0, 12);
0811   TH1D* ScPhiWidth = new TH1D("ScPhiWidth", "ScPhiWidth", 50, 0.005, 0.1);
0812 
0813   TH2D* PhiWidthVSnBC = new TH2D("PhiWidthVSnBC", "PhiWidthVSnBC", 12, 0, 12, 50, 0.005, 0.1);
0814   TH2D* PhiWidthVSfbrem = new TH2D("PhiWidthVSfbrem", "PhiWidthVSfbrem", 100, -0.2, 1., 50, 0.005, 0.1);
0815   TH2D* nBasicClusVSfbrem = new TH2D("nBasicClusVSfbrem", "nBasicClusVSfbrem", 100, -0.2, 1., 12, 0, 12);
0816   TH2D* EopVSfbrem = new TH2D("EopVSfbrem", "EopVSfbrem", 100, -0.2, 1., 100, 0, 3);
0817 
0818   TH1D* EopNegFwd = new TH1D("EopNegFwd", "EopNegFwd", 180, 0, 3);
0819   TH1D* EopNegBwd = new TH1D("EopNegBwd", "EopNegBwd", 180, 0, 3);
0820 
0821   double Mass;
0822   Bool_t MzTag = false;
0823   Bool_t BARREL = true;
0824 
0825   // Loop over tracks
0826   for (Long64_t ientry = 0; ientry < nevent; ientry++) {
0827     // Load the event information
0828     histo.tree->GetEntry(ientry);
0829 
0830     //Control plots
0831     if (track->charge < 0) {
0832       if (track->eta > 0.5)
0833         EopNegFwd->Fill(track->SC_energy / track->p);
0834       if (track->eta < 0.5)
0835         EopNegBwd->Fill(track->SC_energy / track->p);
0836     }
0837 
0838     PhiWidthVSfbrem->Fill(track->fbrem, track->SC_phiWidth);
0839 
0840     // ---- Track Selection ----
0841     NtracksMacro->Fill(0.5);
0842 
0843     //Mz calculation
0844     MzTag = false;
0845     Mass = 0.;
0846     TLorentzVector a(0., 0., 0., 0.);
0847     a.SetXYZM(track->px, track->py, track->pz, 0.511);
0848     TLorentzVector b(0., 0., 0., 0.);
0849     b.SetXYZM(track->px_rejected_track, track->py_rejected_track, track->pz_rejected_track, 0.511);
0850 
0851     Mass = pow((a.E() + b.E()), 2) - pow((a.Px() + b.Px()), 2) - pow((a.Py() + b.Py()), 2) - pow((a.Pz() + b.Pz()), 2);
0852 
0853     Mass = sqrt(Mass);
0854 
0855     if (Mass < 110. && Mass > 70.)
0856       MzTag = true;
0857 
0858     // Z mass window
0859     if (!MzTag)
0860       continue;
0861     cut_Mz->Fill(0.5);
0862 
0863     // -----------        ENDCAP SELECTION    ----------------
0864 
0865     if (!BARREL) {
0866       eta->Fill(track->eta);
0867       if (!track->SC_isEndcap)
0868         continue;
0869       cut_eta->Fill(0.5);
0870 
0871       // Isolation against charged particles
0872       SumPt->Fill(track->SumPtIn05 / track->pt);
0873       if (!track->NoTrackIn0015 || (track->SumPtIn05 / track->pt) > 0.05)
0874         continue;
0875       cut_ChargedIso->Fill(0.5);
0876 
0877       // Ecal energy deposit
0878       Ecal->Fill(track->SC_energy);
0879       if (track->SC_energy < 30.)
0880         continue;
0881       cut_Ecal->Fill(0.5);
0882 
0883       // Hcal over Ecal energy deposit
0884       HoverE->Fill(track->HcalEnergyIn03 / track->SC_energy);
0885       if (track->HcalEnergyIn03 / track->SC_energy > 0.06)
0886         continue;  // before: > 0.08
0887       cut_HoverE->Fill(0.5);
0888 
0889       // Track-SuperCluster matching radius
0890       distTo1stSC->Fill(track->dRto1stSC);
0891       if (track->dRto1stSC > 0.05)
0892         continue;
0893       distTo2ndSC->Fill(track->dRto2ndSC);
0894       if (track->dRto2ndSC < 0.25)
0895         continue;  //min=0.09(1st SC)
0896       cut_NeutralIso->Fill(0.5);
0897 
0898       // a high number of valid hits
0899       if (track->nHits < 13)
0900         continue;
0901       cut_nHits->Fill(0.5);
0902 
0903       // no lost hits
0904       if (track->nLostHits != 0)
0905         continue;
0906       cut_nLostHits->Fill(0.5);
0907 
0908       // outerRadius
0909       //if (track->outerRadius<=99.) continue;
0910       cut_outerRadius->Fill(0.5);
0911 
0912       // good chi2 value
0913       if (track->normalizedChi2 >= 5.)
0914         continue;
0915       cut_normalizedChi2->Fill(0.5);
0916 
0917       // less than 10% bremmstrahlung radiated energy
0918       fbrem->Fill(track->fbrem);
0919       if (track->fbrem > 0.10 || track->fbrem < -0.10)
0920         continue;
0921       cut_fbrem->Fill(0.5);
0922 
0923       // only tracks with associated SuperCluster composed of ONE BasicCluster
0924       nBasicClus->Fill(track->SC_nBasicClus);
0925       if (track->SC_nBasicClus != 1)
0926         continue;
0927       cut_nCluster->Fill(0.5);
0928     }
0929 
0930     // -----------        BARREL SELECTION    ----------------
0931 
0932     if (BARREL) {
0933       eta->Fill(track->eta);
0934       if (!track->SC_isBarrel)
0935         continue;
0936       cut_eta->Fill(0.5);
0937 
0938       // Isolation against charged particles
0939       SumPt->Fill(track->SumPtIn05 / track->pt);
0940       if (!track->NoTrackIn0015 || (track->SumPtIn05 / track->pt) > 0.05)
0941         continue;
0942       cut_ChargedIso->Fill(0.5);
0943 
0944       // Ecal energy deposit
0945       Ecal->Fill(track->SC_energy);
0946       if (track->SC_energy < 25.)
0947         continue;
0948       cut_Ecal->Fill(0.5);
0949 
0950       // Hcal over Ecal energy deposit
0951       HoverE->Fill(track->HcalEnergyIn03 / track->SC_energy);
0952       if (track->HcalEnergyIn03 / track->SC_energy > 0.06)
0953         continue;  // before: > 0.08
0954       cut_HoverE->Fill(0.5);
0955 
0956       // Track-SuperCluster matching radius
0957       distTo1stSC->Fill(track->dRto1stSC);
0958       if (track->dRto1stSC > 0.04)
0959         continue;
0960       distTo2ndSC->Fill(track->dRto2ndSC);
0961       if (track->dRto2ndSC < 0.35)
0962         continue;  //min=0.09(1st SC)
0963       cut_NeutralIso->Fill(0.5);
0964 
0965       // a high number of valid hits
0966       if (track->nHits < 13)
0967         continue;
0968       cut_nHits->Fill(0.5);
0969 
0970       // no lost hits
0971       if (track->nLostHits != 0)
0972         continue;
0973       cut_nLostHits->Fill(0.5);
0974 
0975       // outerRadius
0976       if (track->outerRadius <= 99.)
0977         continue;
0978       cut_outerRadius->Fill(0.5);
0979 
0980       // good chi2 value
0981       if (track->normalizedChi2 >= 5.)
0982         continue;
0983       cut_normalizedChi2->Fill(0.5);
0984 
0985       // less than 10% bremmstrahlung radiated energy
0986       fbrem->Fill(track->fbrem);
0987       if (track->fbrem > 0.1 || track->fbrem < -0.1)
0988         continue;
0989       cut_fbrem->Fill(0.5);
0990 
0991       // only tracks with associated SuperCluster composed of ONE BasicCluster
0992       nBasicClus->Fill(track->SC_nBasicClus);
0993       if (track->SC_nBasicClus != 1)
0994         continue;
0995       cut_nCluster->Fill(0.5);
0996     }
0997 
0998     PhiWidthVSnBC->Fill(track->SC_nBasicClus, track->SC_phiWidth);
0999     nBasicClusVSfbrem->Fill(track->fbrem, track->SC_nBasicClus);
1000     EopVSfbrem->Fill(track->fbrem, track->SC_energy / track->p);
1001 
1002     // only track in studied ECAL range
1003     if (track->SC_energy <= histo.eRange[0] || track->SC_energy >= histo.eRange[histo.eRange.size() - 1])
1004       continue;
1005     cut_EcalRange->Fill(0.5);
1006 
1007     // Loop over eta or phi bins
1008     for (UInt_t j = 0; j < (histo.etaRange.size() - 1); j++) {
1009       // Loop over energy steps
1010       for (UInt_t i = 0; i < (histo.eRange.size() - 1); i++) {
1011         Double_t lowerE = histo.eRange[i];
1012         Double_t upperE = histo.eRange[i + 1];
1013         Double_t lowerEta = histo.etaRange[j];
1014         Double_t upperEta = histo.etaRange[j + 1];
1015 
1016         // Select only tracks in E range
1017         if (!(track->SC_energy > lowerE && track->SC_energy < upperE))
1018           continue;
1019 
1020         // Select only tracks in eta range
1021         if (mode == TRACK_ETA) {
1022           if (!(track->eta >= lowerEta && track->eta < upperEta))
1023             continue;
1024 
1025           usedtracks->Fill(0.5);
1026           histo.selectedTracks[i][j]++;
1027           histo.xAxisBin[i][j]->Fill(radius * 100. / TMath::Tan(2 * TMath::ATan(TMath::Exp(-(track->eta)))));
1028         }
1029 
1030         // Select only tracks in phi range
1031         else if (mode == TRACK_PHI) {
1032           if (!((variable == "phi1" && track->eta < -0.9) || (variable == "phi2" && TMath::Abs(track->eta) < 0.9) ||
1033                 (variable == "phi3" && track->eta > 0.9) || variable == "phi"))
1034             continue;
1035 
1036           if (!(track->phi >= lowerEta && track->phi < upperEta))
1037             continue;
1038 
1039           usedtracks->Fill(0.5);
1040           histo.selectedTracks[i][j]++;
1041           histo.xAxisBin[i][j]->Fill(track->phi);
1042         }
1043 
1044         // e over p
1045         if (track->charge < 0) {
1046           histo.negative[i][j]->Fill((track->SC_energy) / track->p);
1047           histo.Enegative[i][j]->Fill(track->SC_energy * TMath::Sin(track->theta), (track->SC_energy) / track->p);
1048         } else {
1049           histo.positive[i][j]->Fill((track->SC_energy) / track->p);
1050           histo.Epositive[i][j]->Fill(track->SC_energy * TMath::Sin(track->theta), (track->SC_energy) / track->p);
1051         }
1052       }
1053     }
1054   }
1055 
1056   //############ Saving cut table and control histogramms #############
1057 
1058   NtracksMacro->Write();
1059   cut_Ecal->Write();
1060   cut_ChargedIso->Write();
1061   cut_NeutralIso->Write();
1062   cut_HoverE->Write();
1063   cut_nHits->Write();
1064   cut_nLostHits->Write();
1065   cut_outerRadius->Write();
1066   cut_normalizedChi2->Write();
1067   cut_fbrem->Write();
1068   cut_nCluster->Write();
1069   cut_Mz->Write();
1070   cut_EcalRange->Write();
1071   cut_eta->Write();
1072   cut_trigger->Write();
1073   usedtracks->Write();
1074 
1075   Pt->Write();
1076   P->Write();
1077   eta->Write();
1078   Ecal->Write();
1079   HcalEnergyIn01->Write();
1080   HcalEnergyIn02->Write();
1081   HcalEnergyIn03->Write();
1082   HcalEnergyIn04->Write();
1083   HcalEnergyIn05->Write();
1084   SumPt->Write();
1085   HoverE->Write();
1086   distTo1stSC->Write();
1087   distTo2ndSC->Write();
1088   fbrem->Write();
1089   nBasicClus->Write();
1090   ScPhiWidth->Write();
1091   PhiWidthVSnBC->Write();
1092   PhiWidthVSfbrem->Write();
1093   nBasicClusVSfbrem->Write();
1094   EopVSfbrem->Write();
1095 
1096   EopNegFwd->Write();
1097   EopNegBwd->Write();
1098 
1099   //############ CUTS EFFICIENCY ##############
1100 
1101   double NTracksMacro = NtracksMacro->GetBinContent(1);
1102   double usedTracks = usedtracks->GetBinContent(1);
1103   double Cut_Ecal = cut_Ecal->GetBinContent(1);
1104   double Cut_ChargedIso = cut_ChargedIso->GetBinContent(1);
1105   double Cut_NeutralIso = cut_NeutralIso->GetBinContent(1);
1106   double Cut_HoverE = cut_HoverE->GetBinContent(1);
1107   double Cut_nHits = cut_nHits->GetBinContent(1);
1108   double Cut_nLostHits = cut_nLostHits->GetBinContent(1);
1109   double Cut_outerRadius = cut_outerRadius->GetBinContent(1);
1110   double Cut_normalizedChi2 = cut_normalizedChi2->GetBinContent(1);
1111   double Cut_fbrem = cut_fbrem->GetBinContent(1);
1112   double Cut_nCluster = cut_nCluster->GetBinContent(1);
1113   double Cut_Mz = cut_Mz->GetBinContent(1);
1114   double Cut_EcalRange = cut_EcalRange->GetBinContent(1);
1115   double Cut_eta = cut_eta->GetBinContent(1);
1116 
1117   std::cout.setf(std::ios::fixed);
1118   UInt_t precision = std::cout.precision();
1119   std::cout.precision(2);
1120   std::cout << "##################   CUTS EFFICIENCY  ########################" << std::endl;
1121   std::cout << std::endl;
1122   std::cout << "   Cut          | Number of tracks | Subsisting percentage % " << std::endl;
1123   std::cout << "--------------------------------------------------------------" << std::endl;
1124   std::cout << " Trigger selection  | ";
1125   std::cout.width(16);
1126   std::cout << std::right << NTracksMacro << " | ";
1127   std::cout.width(16);
1128   std::cout << NTracksMacro / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1129   std::cout << "--------------------------------------------------------------" << std::endl;
1130   std::cout << " Mz                 | ";
1131   std::cout.width(16);
1132   std::cout << std::right << Cut_Mz << " | ";
1133   std::cout.width(16);
1134   std::cout << Cut_Mz / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1135   std::cout << "--------------------------------------------------------------" << std::endl;
1136   std::cout << " Eta range          | ";
1137   std::cout.width(16);
1138   std::cout << std::right << Cut_eta << " | ";
1139   std::cout.width(16);
1140   std::cout << Cut_eta / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1141   std::cout << "--------------------------------------------------------------" << std::endl;
1142   std::cout << " Iso charged        | ";
1143   std::cout.width(16);
1144   std::cout << std::right << Cut_ChargedIso << " | ";
1145   std::cout.width(16);
1146   std::cout << Cut_ChargedIso / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1147   std::cout << "--------------------------------------------------------------" << std::endl;
1148   std::cout << " Ecal               | ";
1149   std::cout.width(16);
1150   std::cout << std::right << Cut_Ecal << " | ";
1151   std::cout.width(16);
1152   std::cout << Cut_Ecal / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1153   std::cout << "--------------------------------------------------------------" << std::endl;
1154   std::cout << " H over E           | ";
1155   std::cout.width(16);
1156   std::cout << std::right << Cut_HoverE << " | ";
1157   std::cout.width(16);
1158   std::cout << Cut_HoverE / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1159   std::cout << "--------------------------------------------------------------" << std::endl;
1160   std::cout << " Iso Neutral        | ";
1161   std::cout.width(16);
1162   std::cout << std::right << Cut_NeutralIso << " | ";
1163   std::cout.width(16);
1164   std::cout << Cut_NeutralIso / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1165   std::cout << "--------------------------------------------------------------" << std::endl;
1166   std::cout << " nHits              | ";
1167   std::cout.width(16);
1168   std::cout << std::right << Cut_nHits << " | ";
1169   std::cout.width(16);
1170   std::cout << Cut_nHits / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1171   std::cout << "--------------------------------------------------------------" << std::endl;
1172   std::cout << " nLostHits          | ";
1173   std::cout.width(16);
1174   std::cout << std::right << Cut_nLostHits << " | ";
1175   std::cout.width(16);
1176   std::cout << Cut_nLostHits / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1177   std::cout << "--------------------------------------------------------------" << std::endl;
1178   std::cout << " outerRadius        | ";
1179   std::cout.width(16);
1180   std::cout << std::right << Cut_outerRadius << " | ";
1181   std::cout.width(16);
1182   std::cout << Cut_outerRadius / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1183   std::cout << "--------------------------------------------------------------" << std::endl;
1184   std::cout << " normalizedChi2     | ";
1185   std::cout.width(16);
1186   std::cout << std::right << Cut_normalizedChi2 << " | ";
1187   std::cout.width(16);
1188   std::cout << Cut_normalizedChi2 / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1189   std::cout << "--------------------------------------------------------------" << std::endl;
1190   std::cout << " fbrem              | ";
1191   std::cout.width(16);
1192   std::cout << std::right << Cut_fbrem << " | ";
1193   std::cout.width(16);
1194   std::cout << Cut_fbrem / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1195   std::cout << "--------------------------------------------------------------" << std::endl;
1196   std::cout << " nCluster           | ";
1197   std::cout.width(16);
1198   std::cout << std::right << Cut_nCluster << " | ";
1199   std::cout.width(16);
1200   std::cout << Cut_nCluster / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1201   std::cout << "--------------------------------------------------------------" << std::endl;
1202   std::cout << " EcalEnergy Range   | ";
1203   std::cout.width(16);
1204   std::cout << std::right << Cut_EcalRange << " | ";
1205   std::cout.width(16);
1206   std::cout << Cut_EcalRange / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1207   std::cout << "--------------------------------------------------------------" << std::endl;
1208   std::cout << " Used Tracks        | ";
1209   std::cout.width(16);
1210   std::cout << std::right << usedTracks << " | ";
1211   std::cout.width(16);
1212   std::cout << usedTracks / static_cast<Double_t>(NTracksMacro) * 100. << std::endl;
1213   std::cout << std::endl;
1214   std::cout << "##############################################################" << std::endl;
1215   std::cout.unsetf(std::ios::fixed);
1216   std::cout.precision(precision);
1217 }
1218 
1219 // -----------------------------------------------------------------------------
1220 //
1221 //    Auxiliary function : extractgausparams
1222 //
1223 // -----------------------------------------------------------------------------
1224 std::vector<Double_t> extractgausparams(TH1F* histo, TString option1, TString option2) {
1225   // Fitting the histogram with a gaussian function
1226   TF1* f1 = new TF1("f1", "gaus");
1227   f1->SetRange(0., 2.);
1228   histo->Fit("f1", "QR0L");
1229 
1230   // Getting parameters estimated by the fit
1231   double mean = f1->GetParameter(1);
1232   double deviation = f1->GetParameter(2);
1233 
1234   // Iinitializing internal parameters of the iteration algorithms
1235   double lowLim = 0;
1236   double upLim = 0;
1237   double degrade = 0.05;
1238   unsigned int iteration = 0;
1239 
1240   // Iteration procedure
1241   // Iteration is stopped when :
1242   //      more than 10 iterations
1243   //   or fit probability > 0.001
1244   while (iteration == 0 || (f1->GetProb() < 0.001 && iteration < 10)) {
1245     // Computing new bounds for the fit range (2.0 sigma)
1246     lowLim = mean - (2.0 * deviation * (1 - degrade * iteration));
1247     upLim = mean + (2.0 * deviation * (1 - degrade * iteration));
1248     f1->SetRange(lowLim, upLim);
1249 
1250     // Fitting the histo with new bounds
1251     histo->Fit("f1", "QRL0");
1252 
1253     // If the fit succeeds -> extract the new estimated mean value
1254     double newmean = 0;
1255     if (f1->GetParameter(1) < 2 && f1->GetParameter(1) > 0)
1256       newmean = f1->GetParameter(1);
1257     // Else -> keep the previous mean value
1258     else
1259       newmean = mean;
1260 
1261     // Computing new bounds for the fit with the new mean value (1.5 sigma)
1262     lowLim = newmean - (1.5 * deviation);  //*(1-degrade*iteration));
1263     upLim = newmean + (1.5 * deviation);   //*(1-degrade*iteration));
1264     f1->SetRange(lowLim, upLim);
1265     f1->SetLineWidth(1);
1266 
1267     // Fitting the histo with new bounds
1268     histo->Fit("f1", "QRL+i" + option1, option2);
1269 
1270     // if the fit succeeds -> extract the new estimated mean value
1271     if (f1->GetParameter(1) < 2 && f1->GetParameter(1) > 0)
1272       mean = f1->GetParameter(1);
1273 
1274     // Computing new deviation
1275     deviation = f1->GetParameter(2);
1276 
1277     // Next iteration
1278     iteration++;
1279   }
1280 
1281   // return the mean value + its error
1282   std::vector<Double_t> params;
1283   params.push_back(f1->GetParameter(1));
1284   params.push_back(f1->GetParError(1));
1285   params.push_back(lowLim);
1286   params.push_back(upLim);
1287   delete f1;
1288   return params;
1289 }
1290 
1291 // -----------------------------------------------------------------------------
1292 //
1293 //    Auxiliary function : CheckArguments
1294 //
1295 // -----------------------------------------------------------------------------
1296 Bool_t checkArguments(TString variable,  //input
1297                       TString path,
1298                       TString alignmentWithLabel,
1299                       TString outputType,
1300                       Double_t radius,
1301                       Bool_t verbose,
1302                       Double_t givenMin,
1303                       Double_t givenMax,
1304                       ModeType& mode,  // ouput
1305                       std::vector<TFile*>& files,
1306                       std::vector<TString>& labels) {
1307   // Determining mode
1308   if (variable == "eta")
1309     mode = TRACK_ETA;
1310   else if (variable == "phi" || variable == "phi1" || variable == "phi2" || variable == "phi3")
1311     mode = TRACK_PHI;
1312   else {
1313     std::cout << "variable can be eta or phi" << std::endl;
1314     std::cout << "phi1, phi2 and phi3 are for phi in different eta bins" << std::endl;
1315     return false;
1316   }
1317 
1318   // Checking outputType
1319   TString allowed[] = {"eps", "pdf", "svg", "gif", "xpm", "png", "jpg", "tiff", "C", "cxx", "root", "xml"};
1320   bool find = false;
1321   for (unsigned int i = 0; i < 12; i++) {
1322     if (outputType == allowed[i]) {
1323       find = true;
1324       break;
1325     }
1326   }
1327   if (!find) {
1328     std::cout << "ERROR: output type called '" + outputType + "' is unknown !" << std::endl
1329               << "The available output types are : " << std::endl;
1330     for (unsigned int i = 0; i < 12; i++) {
1331       std::cout << " " << allowed[i];
1332     }
1333     std::cout << std::endl;
1334     return false;
1335   }
1336 
1337   // Checking radius
1338   if (radius <= 0) {
1339     std::cout << "ERROR: The radius cannot be null or negative" << std::endl;
1340     return false;
1341   }
1342 
1343   // Checking bounds
1344   if (givenMin > givenMax) {
1345     std::cout << "ERROR: Min value is greater than Max value" << std::endl;
1346     return false;
1347   }
1348 
1349   // Reading the list of files
1350   TObjArray* fileLabelPairs = alignmentWithLabel.Tokenize("\\");
1351   std::cout << "Reading input ROOT files ..." << std::endl;
1352   for (Int_t i = 0; i < fileLabelPairs->GetEntries(); i++) {
1353     TObjArray* singleFileLabelPair = TString(fileLabelPairs->At(i)->GetName()).Tokenize("=");
1354 
1355     if (singleFileLabelPair->GetEntries() == 2) {
1356       TFile* theFile = new TFile(path + (TString(singleFileLabelPair->At(0)->GetName())));
1357       if (!theFile->IsOpen()) {
1358         std::cout << "ERROR : the file '" << theFile->GetName() << "' is not found" << std::endl;
1359         return false;
1360       }
1361       files.push_back(theFile);
1362       labels.push_back(singleFileLabelPair->At(1)->GetName());
1363     } else {
1364       std::cout << "ERROR : Please give file name and legend entry in the following form:\n"
1365                 << " filename1=legendentry1\\filename2=legendentry2\\..." << std::endl;
1366       return false;
1367     }
1368   }
1369 
1370   // Check there is at least of file to analyze
1371   if (files.size() == 0) {
1372     std::cout << "-> No file to analyze" << std::endl;
1373     return false;
1374   }
1375 
1376   // Check labels are unique
1377   for (unsigned int i = 0; i < labels.size(); i++)
1378     for (unsigned int j = i + 1; j < labels.size(); j++) {
1379       if (labels[i] == labels[j]) {
1380         std::cout << "the label '" << labels[i] << "' is twice" << std::endl;
1381         return false;
1382       }
1383     }
1384 
1385   return true;
1386 }
1387 
1388 // -----------------------------------------------------------------------------
1389 //
1390 //    Main (only for stand alone compiled program)
1391 //
1392 // -----------------------------------------------------------------------------
1393 void configureROOTstyle(Bool_t verbose) {
1394   // Configuring style
1395   gROOT->cd();
1396   gROOT->SetStyle("Plain");
1397   if (verbose) {
1398     std::cout << "\n all formerly produced control plots in ./controlEOP/ deleted.\n" << std::endl;
1399     gROOT->ProcessLine(".!if [ -d controlEOP ]; then rm controlEOP/control_eop*; else mkdir controlEOP; fi");
1400   }
1401   gStyle->SetPadBorderMode(0);
1402   gStyle->SetOptStat(0);
1403   gStyle->SetTitleFont(42, "XYZ");
1404   gStyle->SetLabelFont(42, "XYZ");
1405   gStyle->SetEndErrorSize(5);
1406   gStyle->SetLineStyleString(2, "80 20");
1407   gStyle->SetLineStyleString(3, "40 18");
1408   gStyle->SetLineStyleString(4, "20 16");
1409 }
1410 
1411 // -----------------------------------------------------------------------------
1412 //
1413 //    Initialize Tree
1414 //
1415 // -----------------------------------------------------------------------------
1416 Bool_t initializeTree(std::vector<TFile*>& files, std::vector<TTree*>& trees, EopElecVariables* track) {
1417   // resize the tree size
1418   trees.resize(files.size());
1419 
1420   // Loop over the different ROOT files
1421   for (unsigned int i = 0; i < files.size(); i++) {
1422     // Displaying the file name
1423     std::cout << "Opening the file : " << files[i]->GetName() << std::endl;
1424 
1425     // Extracting the TTree from the ROOT file
1426     TTree* theTree = dynamic_cast<TTree*>(files[i]->Get("energyOverMomentumTree/EopTree"));
1427 
1428     // Checking if the TTree is found
1429     if (theTree == 0) {
1430       std::cout << "ERROR : the tree 'EopTree' is not found ! This file is skipped !" << std::endl;
1431       return false;
1432     }
1433 
1434     // Storing the TTree in the container
1435     trees[i] = theTree;
1436 
1437     // Connecting the branches of the TTree to the event object
1438     trees[i]->SetMakeClass(1);
1439     trees[i]->SetBranchAddress("EopElecVariables", &track);
1440     trees[i]->SetBranchAddress("charge", &track->charge);
1441     trees[i]->SetBranchAddress("nHits", &track->nHits);
1442     trees[i]->SetBranchAddress("nLostHits", &track->nLostHits);
1443     trees[i]->SetBranchAddress("innerOk", &track->innerOk);
1444     trees[i]->SetBranchAddress("outerRadius", &track->outerRadius);
1445     trees[i]->SetBranchAddress("chi2", &track->chi2);
1446     trees[i]->SetBranchAddress("normalizedChi2", &track->normalizedChi2);
1447     trees[i]->SetBranchAddress("px_rejected_track", &track->px_rejected_track);
1448     trees[i]->SetBranchAddress("py_rejected_track", &track->py_rejected_track);
1449     trees[i]->SetBranchAddress("pz_rejected_track", &track->pz_rejected_track);
1450     trees[i]->SetBranchAddress("px", &track->px);
1451     trees[i]->SetBranchAddress("py", &track->py);
1452     trees[i]->SetBranchAddress("pz", &track->pz);
1453     trees[i]->SetBranchAddress("p", &track->p);
1454     trees[i]->SetBranchAddress("pIn", &track->pIn);
1455     trees[i]->SetBranchAddress("etaIn", &track->etaIn);
1456     trees[i]->SetBranchAddress("phiIn", &track->phiIn);
1457     trees[i]->SetBranchAddress("pOut", &track->pOut);
1458     trees[i]->SetBranchAddress("etaOut", &track->etaOut);
1459     trees[i]->SetBranchAddress("phiOut", &track->phiOut);
1460     trees[i]->SetBranchAddress("pt", &track->pt);
1461     trees[i]->SetBranchAddress("ptError", &track->ptError);
1462     trees[i]->SetBranchAddress("theta", &track->theta);
1463     trees[i]->SetBranchAddress("eta", &track->eta);
1464     trees[i]->SetBranchAddress("phi", &track->phi);
1465     trees[i]->SetBranchAddress("fbrem", &track->fbrem);
1466     trees[i]->SetBranchAddress("MaxPtIn01", &track->MaxPtIn01);
1467     trees[i]->SetBranchAddress("SumPtIn01", &track->SumPtIn01);
1468     trees[i]->SetBranchAddress("NoTrackIn0015", &track->NoTrackIn0015);
1469     trees[i]->SetBranchAddress("MaxPtIn02", &track->MaxPtIn02);
1470     trees[i]->SetBranchAddress("SumPtIn02", &track->SumPtIn02);
1471     trees[i]->SetBranchAddress("NoTrackIn0020", &track->NoTrackIn0020);
1472     trees[i]->SetBranchAddress("MaxPtIn03", &track->MaxPtIn03);
1473     trees[i]->SetBranchAddress("SumPtIn03", &track->SumPtIn03);
1474     trees[i]->SetBranchAddress("NoTrackIn0025", &track->NoTrackIn0025);
1475     trees[i]->SetBranchAddress("MaxPtIn04", &track->MaxPtIn04);
1476     trees[i]->SetBranchAddress("SumPtIn04", &track->SumPtIn04);
1477     trees[i]->SetBranchAddress("NoTrackIn0030", &track->NoTrackIn0030);
1478     trees[i]->SetBranchAddress("MaxPtIn05", &track->MaxPtIn05);
1479     trees[i]->SetBranchAddress("SumPtIn05", &track->SumPtIn05);
1480     trees[i]->SetBranchAddress("NoTrackIn0035", &track->NoTrackIn0035);
1481     trees[i]->SetBranchAddress("NoTrackIn0040", &track->NoTrackIn0040);
1482     trees[i]->SetBranchAddress("SC_algoID", &track->SC_algoID);
1483     trees[i]->SetBranchAddress("SC_energy", &track->SC_energy);
1484     trees[i]->SetBranchAddress("SC_nBasicClus", &track->SC_nBasicClus);
1485     trees[i]->SetBranchAddress("SC_etaWidth", &track->SC_etaWidth);
1486     trees[i]->SetBranchAddress("SC_phiWidth", &track->SC_phiWidth);
1487     trees[i]->SetBranchAddress("SC_eta", &track->SC_eta);
1488     trees[i]->SetBranchAddress("SC_phi", &track->SC_phi);
1489     trees[i]->SetBranchAddress("SC_isBarrel", &track->SC_isBarrel);
1490     trees[i]->SetBranchAddress("SC_isEndcap", &track->SC_isEndcap);
1491     trees[i]->SetBranchAddress("dRto1stSC", &track->dRto1stSC);
1492     trees[i]->SetBranchAddress("dRto2ndSC", &track->dRto2ndSC);
1493     trees[i]->SetBranchAddress("HcalEnergyIn01", &track->HcalEnergyIn01);
1494     trees[i]->SetBranchAddress("HcalEnergyIn02", &track->HcalEnergyIn02);
1495     trees[i]->SetBranchAddress("HcalEnergyIn03", &track->HcalEnergyIn03);
1496     trees[i]->SetBranchAddress("HcalEnergyIn04", &track->HcalEnergyIn04);
1497     trees[i]->SetBranchAddress("HcalEnergyIn05", &track->HcalEnergyIn05);
1498     trees[i]->SetBranchAddress("isEcalDriven", &track->isEcalDriven);
1499     trees[i]->SetBranchAddress("isTrackerDriven", &track->isTrackerDriven);
1500     trees[i]->SetBranchAddress("RunNumber", &track->RunNumber);
1501     trees[i]->SetBranchAddress("EvtNumber", &track->EvtNumber);
1502   }
1503   return true;
1504 }
1505 
1506 // -----------------------------------------------------------------------------
1507 //
1508 //    Main (only for stand alone compiled program)
1509 //
1510 // -----------------------------------------------------------------------------
1511 int main() {
1512   //momentumBiasValidation("eta","./","EopTree.root=Fit : ","root",1.,true);
1513   momentumBiasValidation("eta",
1514                          "/opt/sbg/data/data1/cms/cgoetzma/Thesis/Eop_withElectrons/Outputs/GsfTracks/",
1515                          "QCD120to170.root=Chris1\\QCD15to30.root=Chris2\\QCD170to300.root=Chris3\\QCD300to470.root="
1516                          "Chris4\\QCD30to50.root=Chris5\\QCD470to600.root=Chris6\\QCD50to80.root=Chris7\\QCD600to800."
1517                          "root=Chris8\\QCD800to1000.root=Chris9\\QCD80to120.root=Chris10\\DyToEE.root=Chris11",
1518                          "root",
1519                          1.,
1520                          true);
1521   return 0;
1522 }