Back to home page

Project CMSSW displayed by LXR
 
 

    

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