Back to home page

Project CMSSW displayed by LXR
 
 

    

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   // give radius at which the misalignment is calculated (arbitrary)
0035   Double_t radius = 1.;
0036 
0037   // configure style
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   // create canvas
0054   TCanvas *c = new TCanvas("c","Canvas",0,0,1150,880);
0055   TCanvas* ccontrol = new TCanvas("ccontrol","controlCanvas",0,0,600,300);
0056 
0057   // create angular binning
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     // tracks beyond abs(eta)=1.65 do not have radii > 99cm (cut value).
0070     startbin = -1.65; lastbin = 1.65; numberOfBins = 18;// ~ outermost TEC
0071     //startbin = -.917; lastbin = .917; numberOfBins = 10;// ~ outermost TOB
0072     //startbin = -1.28; lastbin = 1.28; numberOfBins = 14;// ~ innermost TOB
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   // create energy binning
0087   Int_t startStep = 0;
0088   const Int_t NumberOFSteps = 2;Double_t steparray[NumberOFSteps+1] = {50,58,80};//data 2 steps
0089   //const Int_t NumberOFSteps = 6;Double_t steparray[NumberOFSteps+1] = {50,55,60,65,70,75,80};//mc 6 steps
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   // read files
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   // create trees
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   // create histogram vectors
0153 
0154   // basic histograms
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   // intermediate histograms
0162   vector<TH1F*> fithisto;
0163   vector<TH1F*> combinedxAxisBin;
0164   // final histograms
0165   vector<TGraphErrors*> overallGraph;
0166   vector<TH1F*> overallhisto;
0167 
0168   // book histograms
0169   // basic histograms
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   // intermediate histograms
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   // final histograms
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   // fill basic histograms
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       // track selection
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   // calculate misalignment per bin
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.;// use mean of positive and negative tracks to reduce energy dependence
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       // fill intermediate histograms
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     // calculate mean of different energy bins
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     // fill final histograms
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   // create legend
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     // configure final histogram
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     // fit to final histogram
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));//Divide by 100. cm->m
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     // print fit parameter
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     // write fit function to legend
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     // set pad margins
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     // determine resonable y-axis range
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     // set histogram errors to a small value as only errors of the graph should be shown
0495     for(int i=0; i<overallhisto[isample]->GetNbinsX(); i++){
0496       overallhisto[isample]->SetBinError(i+1,0.00001);
0497     }
0498     // draw final histogram
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   // print plots
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);//*(1-degrade*iteration));
0569       upLim  = newmean + (1.5*deviation);//*(1-degrade*iteration));
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 }
[ 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!