Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-12-16 08:21:18 

0001 #include <TStyle.h>
0002 #include <TCanvas.h>
0003 #include <TTree.h>
0004 
0005 #include <iostream>
0006 #include "TString.h"
0007 #include "TAxis.h"
0008 #include "TProfile.h"
0009 #include "TF1.h"
0010 #include "TH1.h"
0011 #include "TH2F.h"
0012 #include "TGraphErrors.h"
0013 #include "TROOT.h"
0014 #include "TDirectory.h"
0015 #include "TFile.h"
0016 #include "TDirectoryFile.h"
0017 #include "TLegend.h"
0018 #include "TPaveLabel.h"
0019 #include "TChain.h"
0020 #include "TMath.h"
0021 #include "TLatex.h"
0022 #include "TVirtualFitter.h"
0023 #include "TMatrixD.h"
0024 #include "../interface/EopVariables.h"
0025 
0026 using namespace std;
0027 
0028 namespace eop {
0029 
0030   vector<Double_t> extractgausparams(TH1F *histo, TString option1 = "0", TString option2 = "");
0031   vector<Double_t> extractgausparams(TH1F *histo, TString option1, TString option2) {
0032     if (histo->GetEntries() < 10) {
0033       return {0., 0., 0., 0.};
0034     }
0035 
0036     TF1 *f1 = new TF1("f1", "gaus");
0037     f1->SetRange(0., 2.);
0038     histo->Fit("f1", "QR0L");
0039 
0040     double mean = f1->GetParameter(1);
0041     double deviation = f1->GetParameter(2);
0042 
0043     double lowLim;
0044     double upLim;
0045     double newmean;
0046 
0047     double degrade = 0.05;
0048     unsigned int iteration = 0;
0049 
0050     while (iteration == 0 || (f1->GetProb() < 0.001 && iteration < 10)) {
0051       lowLim = mean - (2.0 * deviation * (1 - degrade * iteration));
0052       upLim = mean + (2.0 * deviation * (1 - degrade * iteration));
0053 
0054       f1->SetRange(lowLim, upLim);
0055       histo->Fit("f1", "QRL0");
0056 
0057       newmean = mean;
0058       if (f1->GetParameter(1) < 2 && f1->GetParameter(1) > 0)
0059         newmean = f1->GetParameter(1);
0060       lowLim = newmean - (1.5 * deviation);  //*(1-degrade*iteration));
0061       upLim = newmean + (1.5 * deviation);   //*(1-degrade*iteration));
0062       f1->SetRange(lowLim, upLim);
0063       f1->SetLineWidth(1);
0064       histo->Fit("f1", "QRL+i" + option1, option2);
0065       if (f1->GetParameter(1) < 2 && f1->GetParameter(1) > 0)
0066         mean = f1->GetParameter(1);
0067       deviation = f1->GetParameter(2);
0068       iteration++;
0069     }
0070     vector<Double_t> params;
0071     params.push_back(f1->GetParameter(1));
0072     params.push_back(f1->GetParError(1));
0073     params.push_back(lowLim);
0074     params.push_back(upLim);
0075     return params;
0076   }
0077 
0078   void momentumBiasValidation(
0079       TString variable = "eta",
0080       TString path = "/scratch/hh/current/cms/user/henderle/",
0081       TString alignmentWithLabel =
0082           "EOP_TBDkbMinBias_2011.root=Summer 2011\\EOP_TBD_2011.root=no bows\\EOP_TBDfrom0T_2011.root=from 0T, no "
0083           "bows\\EOP_plain_2011.root=no mass constraint",
0084       bool verbose = false,
0085       Double_t givenMin = 0.,
0086       Double_t givenMax = 0.) {
0087     time_t start = time(0);
0088 
0089     // give radius at which the misalignment is calculated (arbitrary)
0090     Double_t radius = 1.;
0091 
0092     // configure style
0093     gROOT->cd();
0094     gROOT->SetStyle("Plain");
0095     if (verbose) {
0096       std::cout << "\n all formerly produced control plots in ./controlEOP/ deleted.\n" << std::endl;
0097       gROOT->ProcessLine(".!if [ -d controlEOP ]; then rm controlEOP/control_eop*.eps; else mkdir controlEOP; fi");
0098     }
0099     gStyle->SetPadBorderMode(0);
0100     gStyle->SetOptStat(0);
0101     gStyle->SetTitleFont(42, "XYZ");
0102     gStyle->SetLabelFont(42, "XYZ");
0103     gStyle->SetEndErrorSize(5);
0104     gStyle->SetLineStyleString(2, "80 20");
0105     gStyle->SetLineStyleString(3, "40 18");
0106     gStyle->SetLineStyleString(4, "20 16");
0107 
0108     // create canvas
0109     TCanvas *c = new TCanvas("c", "Canvas", 0, 0, 1150, 880);
0110     TCanvas *ccontrol = new TCanvas("ccontrol", "controlCanvas", 0, 0, 600, 300);
0111 
0112     // create angular binning
0113     TString binVar;
0114     if (variable == "eta")
0115       binVar = "track_eta";
0116     else if (variable == "phi" || variable == "phi1" || variable == "phi2" || variable == "phi3")
0117       binVar = "track_phi";
0118     else {
0119       std::cout << "variable can be eta or phi" << std::endl;
0120       std::cout << "phi1, phi2 and phi3 are for phi in different eta bins" << std::endl;
0121     }
0122 
0123     Int_t numberOfBins = 0;
0124     Double_t startbin = 0.;
0125     Double_t lastbin = 0.;
0126     if (binVar == "track_eta") {
0127       // tracks beyond abs(eta)=1.65 do not have radii > 99cm (cut value).
0128       startbin = -1.65;
0129       lastbin = 1.65;
0130       numberOfBins = 18;  // ~ outermost TEC
0131       //startbin = -.917; lastbin = .917; numberOfBins = 10;// ~ outermost TOB
0132       //startbin = -1.28; lastbin = 1.28; numberOfBins = 14;// ~ innermost TOB
0133     }
0134     if (binVar == "track_phi") {
0135       startbin = -TMath::Pi();
0136       lastbin = TMath::Pi();
0137       numberOfBins = 18;
0138     }
0139     Double_t binsize = (lastbin - startbin) / numberOfBins;
0140     vector<Double_t> binningstrvec;
0141     vector<Double_t> zBinning_;
0142     for (Int_t i = 0; i <= numberOfBins; i++) {
0143       binningstrvec.push_back(startbin + i * binsize);
0144       zBinning_.push_back(radius * 100. / TMath::Tan(2 * TMath::ATan(TMath::Exp(-(startbin + i * binsize)))));
0145     }
0146 
0147     // create energy binning
0148     Int_t startStep = 0;
0149     const Int_t NumberOFSteps = 2;
0150     Double_t steparray[NumberOFSteps + 1] = {50, 58, 80};  //data 2 steps
0151     //const Int_t NumberOFSteps = 6;Double_t steparray[NumberOFSteps+1] = {50,55,60,65,70,75,80};//mc 6 steps
0152 
0153     vector<Double_t> steps;
0154     vector<TString> stepstrg;
0155     Char_t tmpstep[10];
0156     for (Int_t ii = 0; ii <= NumberOFSteps; ii++) {
0157       steps.push_back(steparray[ii]);
0158       sprintf(tmpstep, "%1.1fGeV", steparray[ii]);
0159       stepstrg.push_back(tmpstep);
0160     }
0161 
0162     // read files
0163     vector<TFile *> file;
0164     vector<TString> label;
0165 
0166     TObjArray *fileLabelPairs = alignmentWithLabel.Tokenize("\\");
0167     for (Int_t i = 0; i < fileLabelPairs->GetEntries(); ++i) {
0168       TObjArray *singleFileLabelPair = TString(fileLabelPairs->At(i)->GetName()).Tokenize("=");
0169 
0170       if (singleFileLabelPair->GetEntries() == 2) {
0171         file.push_back(new TFile(path + (TString(singleFileLabelPair->At(0)->GetName()))));
0172         label.push_back(singleFileLabelPair->At(1)->GetName());
0173       } else {
0174         std::cout << "Please give file name and legend entry in the following form:\n"
0175                   << " filename1=legendentry1\\filename2=legendentry2\\..." << std::endl;
0176       }
0177     }
0178     cout << "number of files: " << file.size() << endl;
0179 
0180     // create trees
0181     vector<TTree *> tree;
0182     EopVariables *track = new EopVariables();
0183 
0184     for (unsigned int i = 0; i < file.size(); i++) {
0185       tree.push_back((TTree *)file[i]->Get("energyOverMomentumTree/EopTree"));
0186       tree[i]->SetMakeClass(1);
0187       tree[i]->SetBranchAddress("EopVariables", &track);
0188       tree[i]->SetBranchAddress("track_outerRadius", &(track->track_outerRadius));
0189       tree[i]->SetBranchAddress("track_chi2", &track->track_chi2);
0190       tree[i]->SetBranchAddress("track_normalizedChi2", &track->track_normalizedChi2);
0191       tree[i]->SetBranchAddress("track_p", &track->track_p);
0192       tree[i]->SetBranchAddress("track_pt", &track->track_pt);
0193       tree[i]->SetBranchAddress("track_ptError", &track->track_ptError);
0194       tree[i]->SetBranchAddress("track_theta", &track->track_theta);
0195       tree[i]->SetBranchAddress("track_eta", &track->track_eta);
0196       tree[i]->SetBranchAddress("track_phi", &track->track_phi);
0197       tree[i]->SetBranchAddress("track_emc1", &track->track_emc1);
0198       tree[i]->SetBranchAddress("track_emc3", &track->track_emc3);
0199       tree[i]->SetBranchAddress("track_emc5", &track->track_emc5);
0200       tree[i]->SetBranchAddress("track_hac1", &track->track_hac1);
0201       tree[i]->SetBranchAddress("track_hac3", &track->track_hac3);
0202       tree[i]->SetBranchAddress("track_hac5", &track->track_hac5);
0203       tree[i]->SetBranchAddress("track_maxPNearby", &track->track_maxPNearby);
0204       tree[i]->SetBranchAddress("track_EnergyIn", &track->track_EnergyIn);
0205       tree[i]->SetBranchAddress("track_EnergyOut", &track->track_EnergyOut);
0206       tree[i]->SetBranchAddress("distofmax", &track->distofmax);
0207       tree[i]->SetBranchAddress("track_charge", &track->track_charge);
0208       tree[i]->SetBranchAddress("track_nHits", &track->track_nHits);
0209       tree[i]->SetBranchAddress("track_nLostHits", &track->track_nLostHits);
0210       tree[i]->SetBranchAddress("track_innerOk", &track->track_innerOk);
0211     }
0212 
0213     // create histogram vectors
0214 
0215     // basic histograms
0216     vector<TH1F *> histonegative;
0217     vector<TH1F *> histopositive;
0218     vector<TH2F *> Energynegative;
0219     vector<TH2F *> Energypositive;
0220     vector<TH1F *> xAxisBin;
0221     vector<Int_t> selectedTracks;
0222     // intermediate histograms
0223     vector<TH1F *> fithisto;
0224     vector<TH1F *> combinedxAxisBin;
0225     // final histograms
0226     vector<TGraphErrors *> overallGraph;
0227     vector<TH1F *> overallhisto;
0228 
0229     // book histograms
0230     // basic histograms
0231     for (UInt_t i = 0; i < NumberOFSteps * numberOfBins * file.size(); i++) {
0232       Char_t histochar[20];
0233       sprintf(histochar, "%i", i + 1);
0234       TString histostrg = histochar;
0235       TString lowEnergyBorder = stepstrg[i % NumberOFSteps];
0236       TString highEnergyBorder = stepstrg[i % NumberOFSteps + 1];
0237       if (binVar == "track_eta")
0238         histonegative.push_back(
0239             new TH1F("histonegative" + histostrg, lowEnergyBorder + " #leq E < " + highEnergyBorder, 50, 0, 2));
0240       else
0241         histonegative.push_back(
0242             new TH1F("histonegative" + histostrg, lowEnergyBorder + " #leq E_{T} < " + highEnergyBorder, 50, 0, 2));
0243       histopositive.push_back(new TH1F("histopositive" + histostrg, "histopositive" + histostrg, 50, 0, 2));
0244       Energynegative.push_back(
0245           new TH2F("Energynegative" + histostrg, "Energynegative" + histostrg, 5000, 0, 500, 50, 0, 2));
0246       Energypositive.push_back(
0247           new TH2F("Energypositive" + histostrg, "Energypositive" + histostrg, 5000, 0, 500, 50, 0, 2));
0248       xAxisBin.push_back(new TH1F("xAxisBin" + histostrg, "xAxisBin" + histostrg, 1000, -500, 500));
0249       selectedTracks.push_back(0);
0250       Energynegative[i]->Sumw2();
0251       Energypositive[i]->Sumw2();
0252       xAxisBin[i]->Sumw2();
0253       histonegative[i]->Sumw2();
0254       histopositive[i]->Sumw2();
0255       histonegative[i]->SetLineColor(kGreen);
0256       histopositive[i]->SetLineColor(kRed);
0257     }
0258     // intermediate histograms
0259     for (UInt_t i = 0; i < numberOfBins * file.size(); i++) {
0260       Char_t histochar[20];
0261       sprintf(histochar, "%i", i + 1);
0262       TString histostrg = histochar;
0263       fithisto.push_back(new TH1F("fithisto" + histostrg, "fithisto" + histostrg, NumberOFSteps, 0, NumberOFSteps));
0264       combinedxAxisBin.push_back(
0265           new TH1F("combinedxAxisBin" + histostrg, "combinedxAxisBin" + histostrg, NumberOFSteps, 0, NumberOFSteps));
0266     }
0267     // final histograms
0268     for (UInt_t i = 0; i < file.size(); i++) {
0269       TString cnt;
0270       cnt += i;
0271       if (binVar == "track_eta") {
0272         overallhisto.push_back(new TH1F("overallhisto" + cnt, "overallhisto" + cnt, numberOfBins, &zBinning_.front()));
0273         overallGraph.push_back(new TGraphErrors(overallhisto.back()));
0274       } else {
0275         overallhisto.push_back(new TH1F(
0276             "overallhisto" + cnt, "overallhisto" + cnt, numberOfBins, startbin, startbin + numberOfBins * binsize));
0277         overallGraph.push_back(new TGraphErrors(overallhisto.back()));
0278       }
0279     }
0280 
0281     // fill basic histograms
0282     Long64_t nevent;
0283     Int_t usedtracks;
0284     for (UInt_t isample = 0; isample < file.size(); isample++) {
0285       usedtracks = 0;
0286       nevent = (Long64_t)tree[isample]->GetEntries();
0287       cout << nevent << " tracks in sample " << isample + 1 << endl;
0288       for (Long64_t ientry = 0; ientry < nevent; ++ientry) {
0289         tree[isample]->GetEntry(ientry);
0290         for (Int_t i = 0; i < numberOfBins; i++) {
0291           for (Int_t iStep = startStep; iStep < NumberOFSteps; iStep++) {
0292             Double_t lowerEcut = steps[iStep];
0293             Double_t upperEcut = steps[iStep + 1];
0294             Double_t lowcut = binningstrvec[i];
0295             Double_t highcut = binningstrvec[i + 1];
0296             // track selection
0297             if (track->track_nHits >= 13 && track->track_nLostHits == 0 && track->track_outerRadius > 99. &&
0298                 track->track_normalizedChi2 < 5. && track->track_EnergyIn < 1. && track->track_EnergyOut < 8. &&
0299                 track->track_hac3 > lowerEcut && track->track_hac3 < upperEcut) {
0300               if (binVar == "track_eta") {
0301                 if (track->track_eta >= lowcut && track->track_eta < highcut) {
0302                   usedtracks++;
0303                   selectedTracks[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]++;
0304                   xAxisBin[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->Fill(
0305                       radius * 100. / TMath::Tan(2 * TMath::ATan(TMath::Exp(-(track->track_eta)))));
0306                   if (track->track_charge < 0) {
0307                     histonegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->Fill(
0308                         (track->track_hac3) / track->track_p);
0309                     Energynegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->Fill(
0310                         track->track_hac3 * TMath::Sin(track->track_theta), (track->track_hac3) / track->track_p);
0311                   }
0312                   if (track->track_charge > 0) {
0313                     histopositive[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->Fill(
0314                         (track->track_hac3) / track->track_p);
0315                     Energypositive[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->Fill(
0316                         track->track_hac3 * TMath::Sin(track->track_theta), (track->track_hac3) / track->track_p);
0317                   }
0318                 }
0319               } else if (binVar == "track_phi") {
0320                 if ((variable == "phi1" && track->track_eta < -0.9) ||
0321                     (variable == "phi2" && TMath::Abs(track->track_eta) < 0.9) ||
0322                     (variable == "phi3" && track->track_eta > 0.9) || variable == "phi")
0323                   if (track->track_phi >= lowcut && track->track_phi < highcut) {
0324                     usedtracks++;
0325                     selectedTracks[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]++;
0326                     xAxisBin[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->Fill(
0327                         track->track_phi);
0328                     if (track->track_charge < 0) {
0329                       histonegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->Fill(
0330                           (track->track_hac3) / track->track_p);
0331                       Energynegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->Fill(
0332                           track->track_hac3 * TMath::Sin(track->track_theta), (track->track_hac3) / track->track_p);
0333                     }
0334                     if (track->track_charge > 0) {
0335                       histopositive[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->Fill(
0336                           (track->track_hac3) / track->track_p);
0337                       Energypositive[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->Fill(
0338                           track->track_hac3 * TMath::Sin(track->track_theta), (track->track_hac3) / track->track_p);
0339                     }
0340                   }
0341               }
0342             }
0343           }
0344         }
0345       }
0346       cout << "number of used tracks in this sample: " << usedtracks << endl;
0347     }
0348     // calculate misalignment per bin
0349     Double_t misalignment;
0350     Double_t misaliUncert;
0351     vector<TString> misalignmentfit;
0352     vector<TF1 *> f2;
0353 
0354     for (UInt_t isample = 0; isample < file.size(); isample++) {
0355       for (Int_t i = 0; i < numberOfBins; i++) {
0356         if (verbose)
0357           cout << binningstrvec[i] << " < " + binVar + " < " << binningstrvec[i + 1] << endl;
0358         for (Int_t iStep = startStep; iStep < NumberOFSteps; iStep++) {
0359           vector<Double_t> curvNeg;
0360           vector<Double_t> curvPos;
0361           if (verbose) {
0362             TString controlName = "controlEOP/control_eop_sample";
0363             controlName += isample;
0364             controlName += "_bin";
0365             controlName += i;
0366             controlName += "_energy";
0367             controlName += iStep;
0368             controlName += ".eps";
0369             ccontrol->cd();
0370             Double_t posMax =
0371                 histopositive[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->GetMaximum();
0372             Double_t negMax =
0373                 histonegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->GetMaximum();
0374             if (posMax > negMax)
0375               histonegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->SetMaximum(1.1 *
0376                                                                                                               posMax);
0377             histonegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->DrawClone();
0378             histopositive[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->DrawClone("same");
0379             curvNeg = eop::extractgausparams(
0380                 histonegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)], "", "");
0381             curvPos = eop::extractgausparams(
0382                 histopositive[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)], "", "same");
0383             ccontrol->Print(controlName);
0384           } else {
0385             curvNeg = eop::extractgausparams(
0386                 histonegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]);
0387             curvPos = eop::extractgausparams(
0388                 histopositive[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]);
0389           }
0390 
0391           misalignment = 0.;
0392           misaliUncert = 1000.;
0393 
0394           if (selectedTracks[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)] != 0) {
0395             Energynegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]
0396                 ->GetYaxis()
0397                 ->SetRangeUser(curvNeg[2], curvNeg[3]);
0398             Energypositive[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]
0399                 ->GetYaxis()
0400                 ->SetRangeUser(curvPos[2], curvPos[3]);
0401 
0402             Double_t meanEnergyNeg =
0403                 Energynegative[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->GetMean();
0404             Double_t meanEnergyPos =
0405                 Energypositive[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->GetMean();
0406             Double_t meanEnergy = (meanEnergyNeg + meanEnergyPos) /
0407                                   2.;  // use mean of positive and negative tracks to reduce energy dependence
0408             if (verbose)
0409               std::cout << "difference in energy between positive and negative tracks: "
0410                         << meanEnergyNeg - meanEnergyPos << std::endl;
0411 
0412             if (binVar == "track_eta") {
0413               misalignment = 1000000. * 0.5 *
0414                              (-TMath::ASin((0.57 * radius / meanEnergy) * curvNeg[0]) +
0415                               TMath::ASin((0.57 * radius / meanEnergy) * curvPos[0]));
0416               misaliUncert =
0417                   1000000. * 0.5 *
0418                   (TMath::Sqrt(((0.57 * 0.57 * radius * radius * curvNeg[1] * curvNeg[1]) /
0419                                 (meanEnergy * meanEnergy - 0.57 * 0.57 * radius * radius * curvNeg[0] * curvNeg[0])) +
0420                                ((0.57 * 0.57 * radius * radius * curvPos[1] * curvPos[1]) /
0421                                 (meanEnergy * meanEnergy - 0.57 * 0.57 * radius * radius * curvPos[0] * curvPos[0]))));
0422             }
0423             if (binVar == "track_phi") {
0424               misalignment = 1000. * (curvPos[0] - curvNeg[0]) / (curvPos[0] + curvNeg[0]);
0425               misaliUncert = 1000. * 2 / ((curvPos[0] + curvNeg[0]) * (curvPos[0] + curvNeg[0])) *
0426                              TMath::Sqrt((curvPos[0] * curvPos[0] * curvPos[1] * curvPos[1]) +
0427                                          (curvNeg[0] * curvNeg[0] * curvNeg[1] * curvNeg[1]));
0428             }
0429           }
0430           if (verbose)
0431             cout << "misalignment: " << misalignment << "+-" << misaliUncert << endl << endl;
0432           // fill intermediate histograms
0433           fithisto[i + isample * numberOfBins]->SetBinContent(iStep + 1, misalignment);
0434           fithisto[i + isample * numberOfBins]->SetBinError(iStep + 1, misaliUncert);
0435           Double_t xBinCentre =
0436               xAxisBin[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->GetMean();
0437           Double_t xBinCenUnc =
0438               xAxisBin[i * NumberOFSteps + iStep + isample * (NumberOFSteps * numberOfBins)]->GetMeanError();
0439           combinedxAxisBin[i + isample * numberOfBins]->SetBinContent(iStep + 1, xBinCentre);
0440           combinedxAxisBin[i + isample * numberOfBins]->SetBinError(iStep + 1, xBinCenUnc);
0441         }
0442 
0443         Double_t overallmisalignment;
0444         Double_t overallmisaliUncert;
0445         Double_t overallxBin{0.f};
0446         Double_t overallxBinUncert{0.f};
0447         // calculate mean of different energy bins
0448         if (NumberOFSteps > 1 && (fithisto[i + isample * numberOfBins]->Integral() > 0)) {
0449           TF1 *fit = new TF1("fit", "pol0", startStep, NumberOFSteps);
0450           TF1 *fit2 = new TF1("fit2", "pol0", startStep, NumberOFSteps);
0451           if (verbose)
0452             fithisto[i + isample * numberOfBins]->Fit("fit", "0");
0453           else
0454             fithisto[i + isample * numberOfBins]->Fit("fit", "Q0");
0455           combinedxAxisBin[i + isample * numberOfBins]->Fit("fit2", "Q0");
0456           overallmisalignment = fit->GetParameter(0);
0457           overallmisaliUncert = fit->GetParError(0);
0458           overallxBin = fit2->GetParameter(0);
0459           overallxBinUncert = fit2->GetParError(0);
0460           fit->Delete();
0461         } else {
0462           overallmisalignment = misalignment;
0463           overallmisaliUncert = misaliUncert;
0464         }
0465         // fill final histograms
0466         overallhisto[isample]->SetBinContent(i + 1, overallmisalignment);
0467         overallhisto[isample]->SetBinError(i + 1, overallmisaliUncert);
0468         overallGraph[isample]->SetPoint(i, overallxBin, overallmisalignment);
0469         overallGraph[isample]->SetPointError(i, overallxBinUncert, overallmisaliUncert);
0470       }
0471     }
0472 
0473     // create legend
0474     TLegend *leg = new TLegend(0.13, 0.74, 0.98, 0.94);
0475     leg->SetFillColor(10);
0476     leg->SetTextFont(42);
0477     leg->SetTextSize(0.038);
0478     if (variable == "phi1")
0479       leg->SetHeader("low #eta tracks (#eta < -0.9)");
0480     if (variable == "phi2")
0481       leg->SetHeader("central tracks (|#eta| < 0.9)");
0482     if (variable == "phi3")
0483       leg->SetHeader("high #eta tracks (#eta > 0.9)");
0484 
0485     for (UInt_t isample = 0; isample < file.size(); isample++) {
0486       // configure final histogram
0487       if (binVar == "track_eta") {
0488         overallhisto[isample]->GetXaxis()->SetTitle("z [cm]");
0489         overallhisto[isample]->GetYaxis()->SetTitle("#Delta#phi [#murad]");
0490       }
0491       if (binVar == "track_phi") {
0492         overallhisto[isample]->GetXaxis()->SetTitle("#phi");
0493         overallhisto[isample]->GetYaxis()->SetTitle("#Delta#phi [a.u.]");
0494       }
0495       overallhisto[isample]->GetYaxis()->SetTitleOffset(1.05);
0496       overallhisto[isample]->GetYaxis()->SetTitleSize(0.065);
0497       overallhisto[isample]->GetYaxis()->SetLabelSize(0.065);
0498       overallhisto[isample]->GetXaxis()->SetTitleOffset(0.8);
0499       overallhisto[isample]->GetXaxis()->SetTitleSize(0.065);
0500       overallhisto[isample]->GetXaxis()->SetLabelSize(0.065);
0501       overallhisto[isample]->SetLineWidth(2);
0502       overallhisto[isample]->SetLineColor(isample + 1);
0503       overallhisto[isample]->SetMarkerColor(isample + 1);
0504       overallGraph[isample]->SetLineWidth(2);
0505       overallGraph[isample]->SetLineColor(isample + 1);
0506       overallGraph[isample]->SetMarkerColor(isample + 1);
0507       if (isample == 2) {
0508         overallhisto[isample]->SetLineColor(kGreen + 3);
0509         overallhisto[isample]->SetMarkerColor(kGreen + 3);
0510         overallGraph[isample]->SetLineColor(kGreen + 3);
0511         overallGraph[isample]->SetMarkerColor(kGreen + 3);
0512       }
0513       overallGraph[isample]->SetMarkerStyle(isample + 20);
0514       overallGraph[isample]->SetMarkerSize(2);
0515 
0516       // fit to final histogram
0517       Char_t funchar[10];
0518       sprintf(funchar, "func%i", isample + 1);
0519       TString func = funchar;
0520       if (binVar == "track_eta")
0521         f2.push_back(new TF1(func, "[0]+[1]*x/100.", -500, 500));  //Divide by 100. cm->m
0522       if (binVar == "track_phi")
0523         f2.push_back(new TF1(func, "[0]+[1]*TMath::Cos(x+[2])", -500, 500));
0524 
0525       f2[isample]->SetLineColor(isample + 1);
0526       if (isample == 2)
0527         f2[isample]->SetLineColor(kGreen + 3);
0528       f2[isample]->SetLineStyle(isample + 1);
0529       f2[isample]->SetLineWidth(2);
0530 
0531       if (verbose) {
0532         if (overallGraph[isample]->Integral() != 0.f) {
0533           overallGraph[isample]->Fit(func, "mR0+");
0534 
0535           cout << "Covariance Matrix:" << endl;
0536           TVirtualFitter *fitter = TVirtualFitter::GetFitter();
0537           TMatrixD matrix(2, 2, fitter->GetCovarianceMatrix());
0538           Double_t oneOne = fitter->GetCovarianceMatrixElement(0, 0);
0539           Double_t oneTwo = fitter->GetCovarianceMatrixElement(0, 1);
0540           Double_t twoOne = fitter->GetCovarianceMatrixElement(1, 0);
0541           Double_t twoTwo = fitter->GetCovarianceMatrixElement(1, 1);
0542 
0543           cout << "( " << oneOne << ", " << twoOne << ")" << endl;
0544           cout << "( " << oneTwo << ", " << twoTwo << ")" << endl;
0545         }
0546       } else {
0547         if (overallGraph[isample]->Integral() != 0) {
0548           overallGraph[isample]->Fit(func, "QmR0+");
0549         }
0550       }
0551 
0552       // print fit parameter
0553       if (binVar == "track_eta")
0554         cout << "const: " << f2[isample]->GetParameter(0) << "+-" << f2[isample]->GetParError(0)
0555              << ", slope: " << f2[isample]->GetParameter(1) << "+-" << f2[isample]->GetParError(1) << endl;
0556       if (binVar == "track_phi")
0557         cout << "const: " << f2[isample]->GetParameter(0) << "+-" << f2[isample]->GetParError(0)
0558              << ", amplitude: " << f2[isample]->GetParameter(1) << "+-" << f2[isample]->GetParError(1)
0559              << ", shift: " << f2[isample]->GetParameter(2) << "+-" << f2[isample]->GetParError(2) << endl;
0560       cout << "fit probability: " << f2[isample]->GetProb() << endl;
0561       if (verbose)
0562         cout << "chi^2/Ndof: " << f2[isample]->GetChisquare() / f2[isample]->GetNDF() << endl;
0563       // write fit function to legend
0564       Char_t misalignmentfitchar[20];
0565       sprintf(misalignmentfitchar, "%1.f", f2[isample]->GetParameter(0));
0566       misalignmentfit.push_back("(");
0567       misalignmentfit[isample] += misalignmentfitchar;
0568       misalignmentfit[isample] += "#pm";
0569       sprintf(misalignmentfitchar, "%1.f", f2[isample]->GetParError(0));
0570       misalignmentfit[isample] += misalignmentfitchar;
0571       if (variable == "eta") {
0572         misalignmentfit[isample] += ")#murad #upoint r[m] + (";
0573         sprintf(misalignmentfitchar, "%1.f", f2[isample]->GetParameter(1));
0574         misalignmentfit[isample] += misalignmentfitchar;
0575         misalignmentfit[isample] += "#pm";
0576         sprintf(misalignmentfitchar, "%1.f", f2[isample]->GetParError(1));
0577         misalignmentfit[isample] += misalignmentfitchar;
0578         misalignmentfit[isample] += ")#murad #upoint z[m]";
0579       } else if (variable.Contains("phi")) {
0580         misalignmentfit[isample] += ") + (";
0581         sprintf(misalignmentfitchar, "%1.f", f2[isample]->GetParameter(1));
0582         misalignmentfit[isample] += misalignmentfitchar;
0583         misalignmentfit[isample] += "#pm";
0584         sprintf(misalignmentfitchar, "%1.f", f2[isample]->GetParError(1));
0585         misalignmentfit[isample] += misalignmentfitchar;
0586         misalignmentfit[isample] += ") #upoint cos(#phi";
0587         if (f2[isample]->GetParameter(2) > 0.)
0588           misalignmentfit[isample] += "+";
0589         sprintf(misalignmentfitchar, "%1.1f", f2[isample]->GetParameter(2));
0590         misalignmentfit[isample] += misalignmentfitchar;
0591         misalignmentfit[isample] += "#pm";
0592         sprintf(misalignmentfitchar, "%1.1f", f2[isample]->GetParError(2));
0593         misalignmentfit[isample] += misalignmentfitchar;
0594         misalignmentfit[isample] += ")";
0595       }
0596 
0597       // set pad margins
0598       c->cd();
0599       gPad->SetTopMargin(0.06);
0600       gPad->SetBottomMargin(0.12);
0601       gPad->SetLeftMargin(0.13);
0602       gPad->SetRightMargin(0.02);
0603 
0604       // determine resonable y-axis range
0605       Double_t overallmax = 0.;
0606       Double_t overallmin = 0.;
0607       if (isample == 0) {
0608         for (UInt_t i = 0; i < file.size(); i++) {
0609           overallmax = max(
0610               overallhisto[i]->GetMaximum() + overallhisto[i]->GetBinError(overallhisto[i]->GetMaximumBin()) +
0611                   0.55 *
0612                       (overallhisto[i]->GetMaximum() + overallhisto[i]->GetBinError(overallhisto[i]->GetMaximumBin()) -
0613                        overallhisto[i]->GetMinimum() + overallhisto[i]->GetBinError(overallhisto[i]->GetMinimumBin())),
0614               overallmax);
0615           overallmin = min(
0616               overallhisto[i]->GetMinimum() - fabs(overallhisto[i]->GetBinError(overallhisto[i]->GetMinimumBin())) -
0617                   0.1 *
0618                       (overallhisto[i]->GetMaximum() + overallhisto[i]->GetBinError(overallhisto[i]->GetMaximumBin()) -
0619                        overallhisto[i]->GetMinimum() + overallhisto[i]->GetBinError(overallhisto[i]->GetMinimumBin())),
0620               overallmin);
0621         }
0622         overallhisto[isample]->SetMaximum(overallmax);
0623         overallhisto[isample]->SetMinimum(overallmin);
0624         if (givenMax)
0625           overallhisto[isample]->SetMaximum(givenMax);
0626         if (givenMin)
0627           overallhisto[isample]->SetMinimum(givenMin);
0628         overallhisto[isample]->DrawClone("axis");
0629       }
0630       // set histogram errors to a small value as only errors of the graph should be shown
0631       for (int i = 0; i < overallhisto[isample]->GetNbinsX(); i++) {
0632         overallhisto[isample]->SetBinError(i + 1, 0.00001);
0633       }
0634       // draw final histogram
0635       overallhisto[isample]->DrawClone("pe1 same");
0636       f2[isample]->DrawClone("same");
0637       overallGraph[isample]->DrawClone("|| same");
0638       overallGraph[isample]->DrawClone("pz same");
0639       overallGraph[isample]->SetLineStyle(isample + 1);
0640       leg->AddEntry(overallGraph[isample], label[isample] + " (" + misalignmentfit[isample] + ")", "Lp");
0641     }
0642     leg->Draw();
0643 
0644     TPaveLabel *CMSlabel = new TPaveLabel(0.13, 0.94, 0.5, 1., "CMS preliminary 2012", "br NDC");
0645     CMSlabel->SetFillStyle(0);
0646     CMSlabel->SetBorderSize(0);
0647     CMSlabel->SetTextSize(0.95);
0648     CMSlabel->SetTextAlign(12);
0649     CMSlabel->SetTextFont(42);
0650     CMSlabel->Draw("same");
0651 
0652     // print plots
0653     if (variable == "eta") {
0654       c->Print("twist_validation.eps");
0655       if (verbose)
0656         c->Print("twist_validation.root");
0657     } else if (variable == "phi") {
0658       c->Print("sagitta_validation_all.eps");
0659       if (verbose)
0660         c->Print("sagitta_validation_all.root");
0661     } else if (variable == "phi1") {
0662       c->Print("sagitta_validation_lowEta.eps");
0663       if (verbose)
0664         c->Print("sagitta_validation_lowEta.root");
0665     } else if (variable == "phi2") {
0666       c->Print("sagitta_validation_centralEta.eps");
0667       if (verbose)
0668         c->Print("sagitta_validation_centralEta.root");
0669     } else if (variable == "phi3") {
0670       c->Print("sagitta_validation_highEta.eps");
0671       if (verbose)
0672         c->Print("sagitta_validation_highEta.root");
0673     }
0674 
0675     time_t end = time(0);
0676     cout << "Done in " << int(difftime(end, start)) / 60 << " min and " << int(difftime(end, start)) % 60 << " sec."
0677          << endl;
0678   }
0679 
0680 }  // namespace eop
[ 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!