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