Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:57:00 

0001 /* 
0002  Macro to make PV Validation plots
0003 
0004  usage:
0005  $ root -l
0006  $ root [0] .L PlotPVValidation.C++
0007  $ root [1] PlotPVValidation("filename1.root=legendentry1,filename2.root=legendentry2",2,"<your TCut>",Bool_t useBS=false,Bool_t isNormal_=false,Bool_t verbose=false);
0008 
0009  - YourCut is a TCut (default ="")
0010  
0011  // M. Musich - INFN Torino
0012 
0013 */
0014 
0015 #include <string>
0016 #include <sstream>
0017 #include <vector>
0018 #include <Riostream.h>
0019 #include <iomanip>
0020 #include "TFile.h"
0021 #include "TPaveStats.h"
0022 #include "TList.h"
0023 #include "TNtuple.h"
0024 #include "TTree.h"
0025 #include "TH1.h"
0026 #include "TH2.h"
0027 #include "THStack.h"
0028 #include "TStyle.h"
0029 #include "TLegendEntry.h"
0030 #include <TPaveText.h>
0031 #include "TCut.h"
0032 #include "TCanvas.h"
0033 #include "TLegend.h"
0034 #include "TGraphErrors.h"
0035 #include "TF1.h"
0036 #include "TMath.h"
0037 #include "TVectorD.h"
0038 #include <string>
0039 #include <sstream>
0040 #include <vector>
0041 #include <Riostream.h>
0042 #include <iomanip>
0043 #include "TFile.h"
0044 #include "TPaveStats.h"
0045 #include "THStack.h"
0046 #include "TCut.h"
0047 #include "TCanvas.h"
0048 #include "TColor.h"
0049 #include "TLegend.h"
0050 #include "TGraphErrors.h"
0051 #include "TF1.h"
0052 #include "TString.h"
0053 #include "TMath.h"
0054 #include <stdio.h>
0055 #include <iostream>
0056 #include <fstream>
0057 #include <string>
0058 #include <stdio.h>
0059 #include <iostream>
0060 #include <fstream>
0061 #include <string>
0062 #include <TDatime.h>
0063 #include <TSpectrum.h>
0064 
0065 // *------- constants ------* //
0066 const Int_t nplots = 24; 
0067 const Int_t nbins_ = 200;
0068 Int_t _divs = nplots/2.;
0069 Float_t _boundMin   = -0.5;
0070 Float_t _boundSx    = (nplots/4.)-0.5;
0071 Float_t _boundDx    = 3*(nplots/4.)-0.5;
0072 Float_t _boundMax   = nplots-0.5;
0073 
0074 ofstream outfile("FittedDeltaZ.txt");
0075 
0076 TList *FileList;
0077 TList *LabelList;
0078 TFile *Target;
0079 
0080 void MainLoop(TString namesandlabels,Int_t nOfFiles);
0081 
0082 // *------- main function ------* //
0083 void PlotPVValidation(TString namesandlabels,Int_t nOfFiles,TCut theCut="",Bool_t useBS=false,Bool_t isNormal_=false,Bool_t verbose=false);
0084 void fitResiduals(TH1 *hist);
0085 void doubleGausFitResiduals(TH1F *hist);
0086 
0087 Double_t fULine(Double_t *x, Double_t *par);
0088 Double_t fDLine(Double_t *x, Double_t *par);
0089 void FitULine(TH1 *hist);
0090 void FitDLine(TH1 *hist);
0091 
0092 Double_t myCosine(Double_t *x, Double_t *par);
0093 Double_t simplesine(Double_t *x, Double_t *par);
0094 Double_t splitfunc(Double_t *x, Double_t *par);
0095 Double_t splitfunconly(Double_t *x, Double_t *par);
0096 
0097 void MakeNiceTF1Style(TF1 *f1,int color);
0098 void MakeNiceHistoStyle(TH1 *hist);
0099 void MakeNiceTrendPlotStyle(TH1 *hist);
0100 
0101 Float_t calcFWHM(TH1F *hist);
0102 
0103 /*--------------------------------------------------------------------*/
0104 void MainLoop(TString namesandlabels,Int_t nOfFiles)
0105 /*--------------------------------------------------------------------*/
0106 {
0107 
0108   TObjArray *nameandlabelpairs = namesandlabels.Tokenize(",");
0109   TObjArray *aFileLegPair = TString(nameandlabelpairs->At(0)->GetName()).Tokenize("=");  
0110   
0111   TFile *f = TFile::Open(aFileLegPair->At(0)->GetName());
0112   TTree * evTree = (TTree *) f->Get("tree");
0113   unsigned int minRunNumber_=evTree->GetMinimum("RunNumber");
0114   unsigned int maxRunNumber_=evTree->GetMaximum("RunNumber");
0115 
0116   for(UInt_t i=minRunNumber_;i<maxRunNumber_;i++){
0117     char theRunCutstring[128];
0118     sprintf(theRunCutstring,"(RunNumber>=%.1i) && (RunNumber<%.1i)",i,i+1);
0119     std::cout<<theRunCutstring<<std::endl;
0120     TCut theRunCut=theRunCutstring;
0121     PlotPVValidation(namesandlabels,nOfFiles,theRunCut,false,false,false);
0122   }
0123 }
0124 
0125 /*--------------------------------------------------------------------*/
0126 void PlotPVValidation(TString namesandlabels,Int_t nOfFiles, TCut theCut,Bool_t useBS,Bool_t isNormal_,Bool_t verbose)
0127 /*--------------------------------------------------------------------*/
0128 {
0129 
0130   TCut theDefaultCut= "hasRecVertex==1&&isGoodTrack==1";
0131 
0132   TH1::StatOverflows(kTRUE);
0133   gStyle->SetOptTitle(1);
0134   gStyle->SetOptStat("e");
0135   gStyle->SetPadTopMargin(0.12);
0136   gStyle->SetPadBottomMargin(0.15);
0137   gStyle->SetPadLeftMargin(0.17);
0138   gStyle->SetPadRightMargin(0.02);
0139   gStyle->SetPadBorderMode(0);
0140   gStyle->SetTitleFillColor(10);
0141   gStyle->SetTitleFont(42);
0142   gStyle->SetTitleColor(1);
0143   gStyle->SetTitleTextColor(1);
0144   gStyle->SetTitleFontSize(0.06);
0145   gStyle->SetTitleBorderSize(0);
0146   gStyle->SetStatColor(kWhite);
0147   gStyle->SetStatFont(42);
0148   gStyle->SetStatFontSize(0.05);///---> gStyle->SetStatFontSize(0.025);
0149   gStyle->SetStatTextColor(1);
0150   gStyle->SetStatFormat("6.4g");
0151   gStyle->SetStatBorderSize(1);
0152   gStyle->SetPadTickX(1);  // To get tick marks on the opposite side of the frame
0153   gStyle->SetPadTickY(1);
0154   gStyle->SetPadBorderMode(0);
0155   gStyle->SetOptFit(1);
0156   gStyle->SetNdivisions(510);
0157 
0158   using namespace std;
0159   using namespace TMath;
0160 
0161   FileList = new TList();
0162   LabelList = new TList();
0163   
0164   TObjArray *nameandlabelpairs = namesandlabels.Tokenize(",");
0165   for (Int_t i = 0; i < nameandlabelpairs->GetEntries(); ++i) {
0166     TObjArray *aFileLegPair = TString(nameandlabelpairs->At(i)->GetName()).Tokenize("=");
0167     
0168     if(aFileLegPair->GetEntries() == 2) {
0169       FileList->Add( TFile::Open(aFileLegPair->At(0)->GetName())  );  // 2
0170       LabelList->Add( aFileLegPair->At(1) );
0171     }
0172     else {
0173       std::cout << "Please give file name and legend entry in the following form:\n" 
0174         << " filename1=legendentry1,filename2=legendentry2\n";
0175       
0176     }    
0177   }
0178   
0179   Int_t NOfFiles =  FileList->GetSize();
0180   if(NOfFiles!=nOfFiles){
0181     cout<<"&MSG-e: NOfFiles = "<<nOfFiles<<" != "<<NOfFiles<<" given as input"<<endl;
0182     return;
0183   }  
0184   
0185   TTree *trees[nOfFiles]; 
0186   TString LegLabels[nOfFiles];  
0187 
0188   for(Int_t j=0; j < nOfFiles; j++) {
0189     
0190     // Retrieve files
0191 
0192     TFile *fin = (TFile*)FileList->At(j);    
0193     trees[j] = (TTree*)fin->Get("tree");
0194 
0195     // Retrieve labels
0196 
0197     TObjString* legend = (TObjString*)LabelList->At(j);
0198     LegLabels[j] = legend->String();
0199     cout<<"PlotPVValidation() label["<<j<<"]"<<LegLabels[j]<<endl;
0200     
0201   }
0202 
0203   TString label = LegLabels[nOfFiles-1];
0204   label.ReplaceAll(" ","_");
0205   
0206   
0207   /*----------------------------------*/
0208   TString _app;
0209   if(isNormal_) _app = "normal_";
0210   else _app="";
0211   /*----------------------------------*/
0212 
0213   TFile fileout("histos"+_app+label+".root", "new");
0214   
0215   //#################### Canvases #####################
0216 
0217   TCanvas *c1 = new TCanvas("c1","c1",800,550);
0218   c1->SetFillColor(10);  
0219   
0220   TCanvas *c2 = new TCanvas("c2","c2",1460,600);
0221   c2->SetFillColor(10);  
0222   c2->Divide(nplots/2,2);
0223   
0224   TCanvas *c3 = new TCanvas("c3","c3",1460,600);
0225   c3->SetFillColor(10);  
0226   c3->Divide(nplots/2,2);
0227 
0228   TCanvas *c2Eta = new TCanvas("c2Eta","c2Eta",1460,600);
0229   c2Eta->SetFillColor(10);  
0230   c2Eta->Divide(nplots/2,2);
0231   
0232   TCanvas *c3Eta = new TCanvas("c3Eta","c3Eta",1460,600);
0233   c3Eta->SetFillColor(10);  
0234   c3Eta->Divide(nplots/2,2);
0235    
0236   TCanvas *c4 = new TCanvas("c4","c4",1460,550);
0237   c4->SetFillColor(10);  
0238   c4->Divide(2,1);
0239   c4->cd(1)->SetBottomMargin(0.12);
0240   c4->cd(1)->SetLeftMargin(0.13);
0241   c4->cd(2)->SetBottomMargin(0.12);
0242   c4->cd(2)->SetLeftMargin(0.13);
0243 
0244   TCanvas *c4monly  = new TCanvas("c4monly","c4monly",730,550);
0245   c4monly->cd()->SetBottomMargin(0.12);
0246   c4monly->cd()->SetLeftMargin(0.13);
0247   c4monly->SetFillColor(10); 
0248   
0249   TCanvas *c5 = new TCanvas("c5","c5",1460,550);
0250   c5->SetFillColor(10);  
0251   c5->Divide(2,1);
0252   c5->cd(1)->SetBottomMargin(0.12);
0253   c5->cd(1)->SetLeftMargin(0.13);
0254   c5->cd(2)->SetBottomMargin(0.12);
0255   c5->cd(2)->SetLeftMargin(0.13);
0256 
0257   TCanvas *c5monly  = new TCanvas("c5monly","c5monly",730,550);
0258   c5monly->cd()->SetBottomMargin(0.12);
0259   c5monly->cd()->SetLeftMargin(0.13);
0260   c5monly->SetFillColor(10); 
0261 
0262   TCanvas *c4Eta = new TCanvas("c4Eta","c4Eta",1460,550);
0263   c4Eta->SetFillColor(10);  
0264   c4Eta->Divide(2,1);
0265   c4Eta->cd(1)->SetBottomMargin(0.12);
0266   c4Eta->cd(1)->SetLeftMargin(0.13);
0267   c4Eta->cd(2)->SetBottomMargin(0.12);
0268   c4Eta->cd(2)->SetLeftMargin(0.13);
0269   
0270   TCanvas *c5Eta = new TCanvas("c5Eta","c5Eta",1460,550);
0271   c5Eta->SetFillColor(10); 
0272   c5Eta->Divide(2,1);
0273   c5Eta->cd(1)->SetBottomMargin(0.12);
0274   c5Eta->cd(1)->SetLeftMargin(0.13);
0275   c5Eta->cd(2)->SetBottomMargin(0.12);
0276   c5Eta->cd(2)->SetLeftMargin(0.13);
0277                
0278   TLegend *lego = new TLegend(0.15,0.75,0.45,0.89);
0279   lego->SetFillColor(10);
0280   
0281   Int_t colors[8]={kBlack,kRed,kBlue,9,kGreen,5,8};
0282   //Int_t colors[8]={1,2,4,3,8,9,5};
0283   //Int_t colors[8]={6,2,4,1,3,8,9,5};
0284   //Int_t styles[8]={20,4,21,22,25,28,4,29}; 
0285   Int_t styles[8]={20,21,22,23,25,28,4,29}; 
0286 
0287   TH1F *histosdxy[nOfFiles][nplots];
0288   TH1F *histosdz[nOfFiles][nplots];
0289 
0290   TH1F *histosEtadxy[nOfFiles][nplots];
0291   TH1F *histosEtadz[nOfFiles][nplots];
0292 
0293   //############################
0294   // Canvas of bin residuals
0295   //############################
0296   
0297   float phipitch = (2*TMath::Pi())/nplots;
0298   float etapitch = 5./nplots;
0299   //cout<<"etapitch= "<<etapitch;
0300   
0301   for(Int_t i=0; i<nplots; i++){
0302     for(Int_t j=0; j < nOfFiles; j++){
0303 
0304       //for 12 bins
0305       //float phiF=(-TMath::Pi()+i*0.5235)*(180/TMath::Pi());
0306       //float phiL=(-TMath::Pi()+(i+1)*0.5235)*(180/TMath::Pi());
0307       //for 6 bins
0308       //float phiF=(-TMath::Pi()+i*1.047)*(180/TMath::Pi());
0309       //float phiL=(-TMath::Pi()+(i+1)*1.047)*(180/TMath::Pi());
0310 
0311       //for general number of bins
0312   
0313       float phiF=(-TMath::Pi()+i*phipitch)*(180/TMath::Pi());
0314       float phiL=(-TMath::Pi()+(i+1)*phipitch)*(180/TMath::Pi());
0315 
0316       float etaF=-2.5+i*etapitch;
0317       float etaL=-2.5+(i+1)*etapitch;
0318       
0319       float dxymax_phi,dzmax_phi, dxymax_eta, dzmax_eta;
0320       
0321       if(!isNormal_){
0322     if(useBS){
0323       dxymax_phi=1500;
0324       dzmax_phi=300000;
0325       dxymax_eta=1500;
0326       dzmax_eta=300000;
0327     } else {
0328       dxymax_phi=3000;
0329       dzmax_phi=3000;
0330       dxymax_eta=3000;
0331       dzmax_eta=3000;
0332     }
0333       } else {
0334     dxymax_phi=30;
0335     dzmax_phi=30;
0336     dxymax_eta=30;
0337     dzmax_eta=30;
0338       }
0339 
0340       //###################### vs phi #####################
0341 
0342       char histotitledxy[129];
0343       if(!isNormal_) sprintf(histotitledxy,"%.2f #circ<#varphi<%.2f #circ;d_{xy} [#mum];tracks",phiF,phiL);
0344       else sprintf(histotitledxy,"%.2f #circ<#varphi<%.2f #circ;d_{xy}/#sigma_{d_{xy}};tracks",phiF,phiL);
0345       char treeleg1[129];
0346       sprintf(treeleg1,"histodxyfile%i_plot%i",j,i);
0347       histosdxy[j][i] = new TH1F(treeleg1,histotitledxy,nbins_,-dxymax_phi,dxymax_phi);
0348       histosdxy[j][i]->SetLineColor(colors[j]);
0349       MakeNiceHistoStyle(histosdxy[j][i]);
0350     
0351       char histotitledz[129];
0352       if(!isNormal_) sprintf(histotitledz,"%.2f #circ<#varphi<%.2f #circ;d_{z} [#mum];tracks",phiF,phiL);
0353       else sprintf(histotitledz,"%.2f #circ<#varphi<%.2f #circ;d_{z}/#sigma_{d_{z}};tracks",phiF,phiL);
0354       char treeleg2[129];
0355       sprintf(treeleg2,"histodzfile%i_plot%i",j,i);
0356       histosdz[j][i] = new TH1F(treeleg2,histotitledz,nbins_,-dzmax_phi,dzmax_phi);
0357       histosdz[j][i]->SetLineColor(colors[j]);
0358       MakeNiceHistoStyle(histosdz[j][i]);
0359       
0360       //###################### vs eta ##################### 
0361 
0362       char histotitledxyeta[129];
0363       if(!isNormal_) sprintf(histotitledxyeta,"%.2f <#eta<%.2f ;d_{xy} [#mum];tracks",etaF,etaL);
0364       else sprintf(histotitledxyeta,"%.2f <#eta<%.2f;d_{xy}/#sigma_{d_{xy}};tracks",etaF,etaL);
0365       char treeleg1eta[129];
0366       sprintf(treeleg1eta,"histodxyEtafile%i_plot%i",j,i);
0367       histosEtadxy[j][i] = new TH1F(treeleg1eta,histotitledxyeta,nbins_,-dxymax_eta,dxymax_eta);
0368       histosEtadxy[j][i]->SetLineColor(colors[j]);
0369       MakeNiceHistoStyle(histosEtadxy[j][i]);
0370       
0371       char histotitledzeta[129];
0372       if(!isNormal_) sprintf(histotitledzeta,"%.2f <#eta<%.2f ;d_{z} [#mum];tracks",etaF,etaL);
0373       else sprintf(histotitledzeta,"%.2f <#eta<%.2f;d_{z}/#sigma_{d_{z}};tracks",etaF,etaL);
0374       char treeleg2eta[129];
0375       sprintf(treeleg2eta,"histodzEtafile%i_plot%i",j,i);
0376       histosEtadz[j][i] = new TH1F(treeleg2eta,histotitledzeta,nbins_,-dzmax_eta,dzmax_eta);
0377       histosEtadz[j][i]->SetLineColor(colors[j]);
0378       MakeNiceHistoStyle(histosEtadz[j][i]);
0379     }
0380   }
0381   
0382   //################ Phi #####################
0383 
0384   TObject *statObjsdxy[nOfFiles][nplots];
0385   TPaveStats *statsdxy[nOfFiles][nplots];
0386   TObject *statObjsdz[nOfFiles][nplots];
0387   TPaveStats *statsdz[nOfFiles][nplots];
0388   
0389   TH1F *histoMeansdxy[nOfFiles];  
0390   TH1F *histoMeansdz[nOfFiles];     
0391   TH1F *histoSigmasdxy[nOfFiles];   
0392   TH1F *histoSigmasdz[nOfFiles];    
0393   
0394   TF1 *tmpsdxy[nOfFiles][nplots];
0395   TF1 *tmpsdz[nOfFiles][nplots];
0396   
0397   TF1 *FitDzUp[nOfFiles];
0398   TF1 *FitDzDown[nOfFiles];
0399   
0400   TF1 *fleft[nOfFiles]; 
0401   TF1 *fright[nOfFiles];
0402   TF1 *fall[nOfFiles];  
0403 
0404   Double_t meansdxy[nOfFiles][nplots];
0405   Double_t meansdxyError[nOfFiles][nplots];
0406   //Double_t mediansdxy[nOfFiles][nplots];
0407   Double_t sigmasdxy[nOfFiles][nplots];
0408   Double_t sigmasdxyError[nOfFiles][nplots];
0409 
0410   Double_t meansdz[nOfFiles][nplots];
0411   Double_t meansdzError[nOfFiles][nplots];
0412   //Double_t mediansdz[nOfFiles][nplots];
0413   Double_t sigmasdz[nOfFiles][nplots];
0414   Double_t sigmasdzError[nOfFiles][nplots];
0415 
0416   //################ Eta #####################
0417 
0418   TObject *statObjsdxyEta[nOfFiles][nplots];
0419   TPaveStats *statsdxyEta[nOfFiles][nplots];
0420   TObject *statObjsdzEta[nOfFiles][nplots];
0421   TPaveStats *statsdzEta[nOfFiles][nplots]; 
0422 
0423   TH1F *histoMeansdxyEta[nOfFiles];  
0424   TH1F *histoMeansdzEta[nOfFiles];     
0425   TH1F *histoSigmasdxyEta[nOfFiles];   
0426   TH1F *histoSigmasdzEta[nOfFiles];    
0427   
0428   TF1 *tmpsdxyEta[nOfFiles][nplots];
0429   TF1 *tmpsdzEta[nOfFiles][nplots];
0430   
0431   Double_t meansdxyEta[nOfFiles][nplots];
0432   Double_t meansdxyEtaError[nOfFiles][nplots];
0433   //Double_t mediansdxy[nOfFiles][nplots];
0434   Double_t sigmasdxyEta[nOfFiles][nplots];
0435   Double_t sigmasdxyEtaError[nOfFiles][nplots];
0436 
0437   Double_t meansdzEta[nOfFiles][nplots];
0438   Double_t meansdzEtaError[nOfFiles][nplots];
0439   //Double_t mediansdz[nOfFiles][nplots];
0440   Double_t sigmasdzEta[nOfFiles][nplots];
0441   Double_t sigmasdzEtaError[nOfFiles][nplots];
0442 
0443   Double_t highedge=nplots-0.5;
0444   Double_t lowedge=-0.5;
0445 
0446   //#################### initialization limits ################################
0447   
0448   Double_t dxymean_phi, dxymean_eta, dxysigma_phi, dxysigma_eta, dzmean_phi, dzmean_eta, dzsigma_phi, dzsigma_eta ;
0449 
0450   if(!isNormal_){
0451     if(useBS){
0452       dxymean_phi=40; 
0453       dxymean_eta=40; 
0454       dxysigma_phi=310; 
0455       dxysigma_eta=510; 
0456       dzmean_phi=2000; 
0457       dzmean_eta=2000; 
0458       dzsigma_phi=100000; 
0459       dzsigma_eta=100000;
0460     } else {
0461       dxymean_phi=40; 
0462       dxymean_eta=40; 
0463       dxysigma_phi=260; 
0464       dxysigma_eta=260; 
0465       dzmean_phi=50; 
0466       dzmean_eta=160; 
0467       dzsigma_phi=410; 
0468       dzsigma_eta=1000;
0469     }
0470   } else {
0471     dxymean_phi=0.5; 
0472     dxymean_eta=0.5; 
0473     dxysigma_phi=2.; 
0474     dxysigma_eta=2.; 
0475     dzmean_phi=0.5; 
0476     dzmean_eta=0.5; 
0477     dzsigma_phi=2.; 
0478     dzsigma_eta=2.;
0479   }
0480 
0481   for(Int_t j=0; j < nOfFiles; j++) { 
0482 
0483     //################ phi #########################
0484 
0485     TString temptitle;
0486     if(!isNormal_) temptitle = "#LT d_{xy} #GT vs #varphi sector";
0487     else temptitle = "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #varphi sector";
0488 
0489     char meanslegdxy[129];
0490     sprintf(meanslegdxy,"histomeansdxy_file%i",j);
0491     histoMeansdxy[j]  = new TH1F(meanslegdxy,temptitle,nplots,lowedge,highedge);
0492     histoMeansdxy[j]->SetLineColor(colors[j]);
0493     histoMeansdxy[j]->GetYaxis()->SetRangeUser(-dxymean_phi,dxymean_phi);
0494     histoMeansdxy[j]->SetMarkerStyle(styles[j]);
0495     histoMeansdxy[j]->SetMarkerColor(colors[j]);
0496     histoMeansdxy[j]->GetXaxis()->SetTitle("#varphi (sector) [degrees]");
0497     if(!isNormal_) histoMeansdxy[j]->GetYaxis()->SetTitle ("#LT d_{xy} #GT [#mum]");
0498     else histoMeansdxy[j]->GetYaxis()->SetTitle ("#LT d_{xy}/#sigma_{d_{xy}} #GT");
0499     MakeNiceTrendPlotStyle(histoMeansdxy[j]);
0500     
0501     temptitle.Clear();
0502     if(!isNormal_) temptitle = "#LT d_{z} #GT vs #varphi sector";
0503     else temptitle = "#LT d_{z}/#sigma_{d_{z}} #GT vs #varphi sector";
0504 
0505     char meanslegdz[129];
0506     sprintf(meanslegdz,"histomeansdz_file%i",j);
0507     histoMeansdz[j]   = new TH1F(meanslegdz,temptitle,nplots,lowedge,highedge); 
0508     histoMeansdz[j]->SetLineColor(colors[j]);
0509     histoMeansdz[j]->GetYaxis()->SetRangeUser(-dzmean_phi,dzmean_phi);
0510     histoMeansdz[j]->SetMarkerStyle(styles[j]);
0511     histoMeansdz[j]->SetMarkerColor(colors[j]);
0512     histoMeansdz[j]->GetXaxis()->SetTitle("#varphi (sector) [degrees]");
0513     if(!isNormal_) histoMeansdz[j]->GetYaxis()->SetTitle ("#LT d_{z} #GT [#mum]");
0514     else histoMeansdz[j]->GetYaxis()->SetTitle ("#LT d_{z}/#sigma_{d_{z}} #GT");
0515     MakeNiceTrendPlotStyle(histoMeansdz[j]);
0516 
0517     temptitle.Clear();
0518     if(!isNormal_) temptitle = "#sigma_{d_{xy}} vs #varphi sector";
0519     else temptitle = "width(d_{xy}/#sigma_{d_{xy}}) vs #varphi sector";
0520 
0521     char sigmaslegdxy[129];
0522     sprintf(sigmaslegdxy,"histosigmasdxy_file%i",j);
0523     histoSigmasdxy[j] = new TH1F(sigmaslegdxy,temptitle,nplots,lowedge,highedge);
0524     histoSigmasdxy[j]->SetLineColor(colors[j]);
0525     histoSigmasdxy[j]->GetYaxis()->SetRangeUser(0,dxysigma_phi);
0526     histoSigmasdxy[j]->SetMarkerStyle(styles[j]);
0527     histoSigmasdxy[j]->SetMarkerColor(colors[j]);
0528     histoSigmasdxy[j]->GetXaxis()->SetTitle("#varphi (sector) [degrees]");
0529     if(!isNormal_) histoSigmasdxy[j]->GetYaxis()->SetTitle ("#sigma_{d_{xy}} [#mum]");
0530     else histoSigmasdxy[j]->GetYaxis()->SetTitle ("width(d_{xy}/#sigma_{d_{xy}})");
0531     MakeNiceTrendPlotStyle(histoSigmasdxy[j]);
0532 
0533     temptitle.Clear();
0534     if(!isNormal_) temptitle = "#sigma_{d_{z}} vs #varphi sector";
0535     else temptitle = "width(d_{z}/#sigma_{d_{z}}) vs #varphi sector";
0536 
0537     char sigmaslegdz[129];
0538     sprintf(sigmaslegdz,"histosigmasdz_file%i",j);
0539     histoSigmasdz[j]  = new TH1F(sigmaslegdz,temptitle,nplots,lowedge,highedge);
0540     histoSigmasdz[j]->SetLineColor(colors[j]);
0541     histoSigmasdz[j]->GetYaxis()->SetRangeUser(0,dzsigma_phi);
0542     histoSigmasdz[j]->SetMarkerStyle(styles[j]);
0543     histoSigmasdz[j]->SetMarkerColor(colors[j]);
0544     histoSigmasdz[j]->GetXaxis()->SetTitle("#varphi (sector) [degrees]");
0545     if(!isNormal_) histoSigmasdz[j]->GetYaxis()->SetTitle ("#sigma_{d_{z}} [#mum]");
0546     else histoSigmasdz[j]->GetYaxis()->SetTitle ("width(d_{z}/#sigma_{d_{z}})");
0547     MakeNiceTrendPlotStyle(histoSigmasdz[j]);
0548 
0549     //################ eta #########################
0550   
0551     temptitle.Clear();
0552     if(!isNormal_) temptitle = "#LT d_{xy} #GT vs #eta sector"; 
0553     else temptitle = "#LT  d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector";
0554 
0555     char meanslegdxyEta[129];
0556     sprintf(meanslegdxyEta,"histomeansEtadxy_file%i",j);
0557     histoMeansdxyEta[j]  = new TH1F(meanslegdxyEta,temptitle,nplots,lowedge,highedge);
0558     histoMeansdxyEta[j]->SetLineColor(colors[j]);
0559     histoMeansdxyEta[j]->GetYaxis()->SetRangeUser(-dxymean_eta,dxymean_eta);
0560     histoMeansdxyEta[j]->SetMarkerStyle(styles[j]);
0561     histoMeansdxyEta[j]->SetMarkerColor(colors[j]);
0562     histoMeansdxyEta[j]->GetXaxis()->SetTitle("#eta (sector)");
0563     if(!isNormal_) histoMeansdxyEta[j]->GetYaxis()->SetTitle ("#LT d_{xy} #GT [#mum]");
0564     else histoMeansdxyEta[j]->GetYaxis()->SetTitle ("#LT d_{xy}/#sigma_{d_{xy}} #GT");
0565     MakeNiceTrendPlotStyle(histoMeansdxyEta[j]);
0566    
0567     temptitle.Clear();
0568     if(!isNormal_) temptitle = "#LT d_{z} #GT vs #eta sector"; 
0569     else temptitle = "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector";
0570 
0571     char meanslegdzEta[129];
0572     sprintf(meanslegdzEta,"histomeansEtadz_file%i",j);
0573     histoMeansdzEta[j]   = new TH1F(meanslegdzEta,temptitle,nplots,lowedge,highedge); 
0574     histoMeansdzEta[j]->SetLineColor(colors[j]);
0575     histoMeansdzEta[j]->GetYaxis()->SetRangeUser(-dzmean_eta,dzmean_eta);
0576     histoMeansdzEta[j]->SetMarkerStyle(styles[j]);
0577     histoMeansdzEta[j]->SetMarkerColor(colors[j]);
0578     histoMeansdzEta[j]->GetXaxis()->SetTitle("#eta (sector)");
0579     if(!isNormal_) histoMeansdzEta[j]->GetYaxis()->SetTitle ("#LT d_{z} #GT [#mum]");
0580     else histoMeansdzEta[j]->GetYaxis()->SetTitle ("#LT d_{z}/#sigma_{d_{z}} #GT");
0581     MakeNiceTrendPlotStyle(histoMeansdzEta[j]);
0582 
0583     temptitle.Clear();
0584     if(!isNormal_) temptitle = "#sigma_{d_{xy}} vs #eta sector"; 
0585     else temptitle = "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector";
0586 
0587     char sigmaslegdxyEta[129];
0588     sprintf(sigmaslegdxyEta,"histosigmasEtadxy_file%i",j);
0589     histoSigmasdxyEta[j] = new TH1F(sigmaslegdxyEta,temptitle,nplots,lowedge,highedge);
0590     histoSigmasdxyEta[j]->SetLineColor(colors[j]);
0591     histoSigmasdxyEta[j]->GetYaxis()->SetRangeUser(0,dxysigma_eta);
0592     histoSigmasdxyEta[j]->SetMarkerStyle(styles[j]);
0593     histoSigmasdxyEta[j]->SetMarkerColor(colors[j]);
0594     histoSigmasdxyEta[j]->GetXaxis()->SetTitle("#eta (sector)");
0595     if(!isNormal_) histoSigmasdxyEta[j]->GetYaxis()->SetTitle ("#sigma_{d_{xy}} [#mum]");
0596     else histoSigmasdxyEta[j]->GetYaxis()->SetTitle ("width(d_{xy}/#sigma_{d_{xy}})");
0597     MakeNiceTrendPlotStyle(histoSigmasdxyEta[j]);
0598 
0599     temptitle.Clear();
0600     if(!isNormal_) temptitle = "#sigma_{d_{z}} vs #eta sector"; 
0601     else temptitle = "width(d_{z}/#sigma_{d_{z}}) vs #eta sector";
0602 
0603     char sigmaslegdzEta[129];
0604     sprintf(sigmaslegdzEta,"histosigmasEtadz_file%i",j);
0605     histoSigmasdzEta[j]  = new TH1F(sigmaslegdzEta,temptitle,nplots,lowedge,highedge);
0606     histoSigmasdzEta[j]->SetLineColor(colors[j]);
0607     histoSigmasdzEta[j]->GetYaxis()->SetRangeUser(0,dzsigma_eta);
0608     histoSigmasdzEta[j]->SetMarkerStyle(styles[j]);
0609     histoSigmasdzEta[j]->SetMarkerColor(colors[j]);
0610     histoSigmasdzEta[j]->GetXaxis()->SetTitle("#eta (sector)");
0611     if(!isNormal_) histoSigmasdzEta[j]->GetYaxis()->SetTitle ("#sigma_{d_{z}} [#mum]");
0612     else histoSigmasdzEta[j]->GetYaxis()->SetTitle ("width(d_{z}/#sigma_{d_{z}})");
0613     MakeNiceTrendPlotStyle(histoSigmasdzEta[j]);
0614   }
0615 
0616   
0617    //Histos maxima
0618    
0619    double hmaxDxyphi=0.;
0620    double hmaxDxyeta=0.;
0621    double hmaxDzphi=0.;
0622    double hmaxDzeta=0.;
0623 
0624    //######################## dxy Residuals #################################
0625    
0626    for(Int_t j=0; j < nOfFiles; j++) { 
0627      std::cout<<"PlotPVValidation() Running on file: "<<LegLabels[j]<<std::endl;
0628      for(Int_t i=0; i<nplots; i++){
0629        if(verbose){ 
0630      std::cout<<"PlotPVValidation() Fitting the : "<<i<<"/"<<nplots<<" bin ";
0631      TDatime *time = new TDatime();
0632      time->Print();
0633      //std::cout<<std::endl;
0634        }
0635 
0636        //################### phi ########################
0637        
0638        char thePhiCutstring[128];
0639        sprintf(thePhiCutstring,"(phi>-TMath::Pi()+%.1i*(%.3f))&&(phi<-TMath::Pi()+%.1i*(%.3f))",i,phipitch,i+1,phipitch);
0640        TCut thePhiCut = thePhiCutstring;
0641        //cout<<thePhiCutstring<<endl;
0642        
0643        char treeleg3[129];
0644 
0645        if(!isNormal_){
0646      if(useBS){
0647        sprintf(treeleg3,"dxyBs*10000>>histodxyfile%i_plot%i",j,i);
0648      } else {
0649        sprintf(treeleg3,"dxyFromMyVertex*10000>>histodxyfile%i_plot%i",j,i);
0650      } 
0651        } else {
0652      sprintf(treeleg3,"IPTsigFromMyVertex>>histodxyfile%i_plot%i",j,i);
0653        }
0654 
0655        char phipositionString[129];
0656        float_t phiInterval = (360.)/nplots;
0657        float phiposition = (-180+i*phiInterval)+(phiInterval/2);
0658        sprintf(phipositionString,"%.f",phiposition);
0659        //cout<<"phiposition: "<<phiposition<<" phipositionString: "<<phipositionString<<endl;
0660 
0661        c2->cd(i+1);
0662        if(j==0){
0663      trees[j]->Draw(treeleg3,theCut&&theDefaultCut&&thePhiCut);  
0664      hmaxDxyphi=histosdxy[0][i]->GetMaximum();
0665        }
0666        else{ 
0667      trees[j]->Draw(treeleg3,theCut&&theDefaultCut&&thePhiCut,"sames"); 
0668      
0669      if(histosdxy[j][i]->GetMaximum() >  hmaxDxyphi){
0670        hmaxDxyphi = histosdxy[j][i]->GetMaximum();
0671      }
0672        }
0673        
0674        histosdxy[0][i]->GetYaxis()->SetRangeUser(0.01,hmaxDxyphi*1.10);
0675        histosdxy[0][i]->Draw("sames");
0676        
0677        if(j!=0){
0678      histosdxy[j][i]->Draw("sames");
0679        }
0680 
0681        histosdxy[j][i]->Write();
0682 
0683        fitResiduals(histosdxy[j][i]);   
0684        tmpsdxy[j][i] = (TF1*)histosdxy[j][i]->GetListOfFunctions()->FindObject("tmp");   
0685        if(tmpsdxy[j][i])  {
0686      tmpsdxy[j][i]->SetLineColor(colors[j]);
0687      tmpsdxy[j][i]->SetLineWidth(2); 
0688      tmpsdxy[j][i]->Draw("same"); 
0689        
0690      meansdxy[j][i]=(tmpsdxy[j][i])->GetParameter(1);
0691      sigmasdxy[j][i]=(tmpsdxy[j][i])->GetParameter(2);
0692      
0693      meansdxyError[j][i]=(tmpsdxy[j][i])->GetParError(1);
0694      sigmasdxyError[j][i]=(tmpsdxy[j][i])->GetParError(2); 
0695      
0696      histoMeansdxy[j]->SetBinContent(i+1,meansdxy[j][i]); 
0697      histoMeansdxy[j]->SetBinError(i+1,meansdxyError[j][i]); 
0698      histoMeansdxy[j]->GetXaxis()->SetBinLabel(i+1,phipositionString);
0699      histoSigmasdxy[j]->SetBinContent(i+1,sigmasdxy[j][i]);
0700      histoSigmasdxy[j]->SetBinError(i+1,sigmasdxyError[j][i]); 
0701      histoSigmasdxy[j]->GetXaxis()->SetBinLabel(i+1,phipositionString);
0702        }
0703 
0704        else{
0705      histoMeansdxy[j]->SetBinContent(i+1,0.); 
0706      histoMeansdxy[j]->SetBinError(i+1,0.);
0707      histoMeansdxy[j]->GetXaxis()->SetBinLabel(i+1,phipositionString);
0708      histoSigmasdxy[j]->SetBinContent(i+1,0.);
0709      histoSigmasdxy[j]->SetBinError(i+1,0.); 
0710      histoSigmasdxy[j]->GetXaxis()->SetBinLabel(i+1,phipositionString); 
0711        }
0712        
0713        c2->Draw("sames");   
0714        c2->cd(i+1);
0715        
0716        statObjsdxy[j][i] = histosdxy[j][i]->GetListOfFunctions()->FindObject("stats");
0717        statsdxy[j][i] = static_cast<TPaveStats*>(statObjsdxy[j][i]);
0718        statsdxy[j][i]->SetLineColor(colors[j]);
0719        statsdxy[j][i]->SetTextColor(colors[j]);
0720        statsdxy[j][i]->SetShadowColor(10);
0721        statsdxy[j][i]->SetFillColor(10);
0722        statsdxy[j][i]->SetX1NDC(0.66);
0723        statsdxy[j][i]->SetY1NDC(0.80-j*0.19);
0724        statsdxy[j][i]->SetX2NDC(0.99);
0725        statsdxy[j][i]->SetY2NDC(0.99-j*0.19);
0726        statsdxy[j][i]->Draw("same");   
0727   
0728        //################### eta ########################
0729        
0730        char theEtaCutstring[128];
0731        sprintf(theEtaCutstring,"(eta>-2.5+%.1i*(%.3f))&&(eta<-2.5+%.1i*(%.3f))",i,etapitch,i+1,etapitch);
0732        TCut theEtaCut = theEtaCutstring;
0733        //cout<<theEtaCutstring<<endl;
0734        
0735        char treeleg3Eta[129];
0736 
0737        if(!isNormal_){
0738      if(useBS){
0739        sprintf(treeleg3Eta,"dxyBs*10000>>histodxyEtafile%i_plot%i",j,i);
0740      } else {
0741        sprintf(treeleg3Eta,"dxyFromMyVertex*10000>>histodxyEtafile%i_plot%i",j,i);
0742      }
0743        } else {
0744      sprintf(treeleg3Eta,"IPTsigFromMyVertex>>histodxyEtafile%i_plot%i",j,i);
0745        }      
0746 
0747        char etapositionString[129];
0748        float etaposition = (-2.5+i*(5./nplots))+((5./nplots)/2);
0749        sprintf(etapositionString,"%.1f",etaposition);
0750        //cout<<"etaposition: "<<etaposition<<" etapositionString: "<<etapositionString<<endl;
0751        
0752        c2Eta->cd(i+1);
0753        if(j==0){
0754      trees[j]->Draw(treeleg3Eta,theCut&&theDefaultCut&&theEtaCut);
0755      hmaxDxyeta=histosEtadxy[0][i]->GetMaximum();  
0756        }
0757        else{ 
0758      trees[j]->Draw(treeleg3Eta,theCut&&theDefaultCut&&theEtaCut,"sames"); 
0759      if(histosEtadxy[j][i]->GetMaximum() >  hmaxDxyeta){
0760        hmaxDxyeta = histosEtadxy[j][i]->GetMaximum();
0761      }
0762        }
0763 
0764        histosEtadxy[0][i]->GetYaxis()->SetRangeUser(0.01,hmaxDxyeta*1.10);
0765        histosEtadxy[0][i]->Draw("sames");
0766        
0767        if(j!=0){
0768      histosEtadxy[j][i]->Draw("sames");
0769        }
0770 
0771        histosEtadxy[j][i]->Write();
0772 
0773        fitResiduals(histosEtadxy[j][i]);   
0774        tmpsdxyEta[j][i] = (TF1*)histosEtadxy[j][i]->GetListOfFunctions()->FindObject("tmp");   
0775        if(tmpsdxyEta[j][i])  {
0776      tmpsdxyEta[j][i]->SetLineColor(colors[j]);
0777      tmpsdxyEta[j][i]->SetLineWidth(2); 
0778      
0779      meansdxyEta[j][i]=(tmpsdxyEta[j][i])->GetParameter(1);
0780      sigmasdxyEta[j][i]=(tmpsdxyEta[j][i])->GetParameter(2);
0781      
0782      meansdxyEtaError[j][i]=(tmpsdxyEta[j][i])->GetParError(1);
0783      sigmasdxyEtaError[j][i]=(tmpsdxyEta[j][i])->GetParError(2);
0784      
0785      histoMeansdxyEta[j]->SetBinContent(i+1,meansdxyEta[j][i]); 
0786      histoMeansdxyEta[j]->SetBinError(i+1,meansdxyEtaError[j][i]); 
0787      histoMeansdxyEta[j]->GetXaxis()->SetBinLabel(i+1,etapositionString);
0788      histoSigmasdxyEta[j]->SetBinContent(i+1,sigmasdxyEta[j][i]);
0789      histoSigmasdxyEta[j]->SetBinError(i+1,sigmasdxyEtaError[j][i]); 
0790      histoSigmasdxyEta[j]->GetXaxis()->SetBinLabel(i+1,etapositionString);
0791        } else{
0792      histoMeansdxyEta[j]->SetBinContent(i+1,0.); 
0793      histoMeansdxyEta[j]->SetBinError(i+1,0.);   
0794      histoMeansdxyEta[j]->GetXaxis()->SetBinLabel(i+1,etapositionString);
0795      histoSigmasdxyEta[j]->SetBinContent(i+1,0.);
0796      histoSigmasdxyEta[j]->SetBinError(i+1,0.);
0797      histoSigmasdxyEta[j]->GetXaxis()->SetBinLabel(i+1,etapositionString);
0798        }
0799 
0800        c2Eta->Draw("sames");   
0801        c2Eta->cd(i+1);
0802        
0803        statObjsdxyEta[j][i] = histosEtadxy[j][i]->GetListOfFunctions()->FindObject("stats");
0804        statsdxyEta[j][i] = static_cast<TPaveStats*>(statObjsdxyEta[j][i]);
0805        statsdxyEta[j][i]->SetLineColor(colors[j]);
0806        statsdxyEta[j][i]->SetTextColor(colors[j]);
0807        statsdxyEta[j][i]->SetShadowColor(10);
0808        statsdxyEta[j][i]->SetFillColor(10);
0809        statsdxyEta[j][i]->SetX1NDC(0.66);
0810        statsdxyEta[j][i]->SetY1NDC(0.80-j*0.18);
0811        statsdxyEta[j][i]->SetX2NDC(0.99);
0812        statsdxyEta[j][i]->SetY2NDC(0.99-j*0.19);
0813        statsdxyEta[j][i]->Draw("same"); 
0814        
0815        //######################## dz Residuals #################################
0816        
0817        char treeleg4[129];
0818        
0819        if(!isNormal_){
0820      if(useBS){
0821        sprintf(treeleg4,"dzBs*10000>>histodzfile%i_plot%i",j,i);
0822      } else {
0823        sprintf(treeleg4,"dzFromMyVertex*10000>>histodzfile%i_plot%i",j,i);
0824      }
0825        } else {
0826      sprintf(treeleg4,"IPLsigFromMyVertex>>histodzfile%i_plot%i",j,i);
0827        }        
0828 
0829        c3->cd(i+1);
0830        if(j==0){
0831      trees[j]->Draw(treeleg4,theCut&&theDefaultCut&&thePhiCut);
0832      hmaxDzphi=histosdz[0][i]->GetMaximum();  
0833        }
0834        
0835        else{ 
0836      trees[j]->Draw(treeleg4,theCut&&theDefaultCut&&thePhiCut,"sames"); 
0837      if(histosdz[j][i]->GetMaximum() >  hmaxDzphi){
0838        hmaxDzphi = histosdz[j][i]->GetMaximum();
0839      }
0840        }
0841        
0842        histosdz[0][i]->GetYaxis()->SetRangeUser(0.01,hmaxDzphi*1.10);
0843        histosdz[0][i]->Draw("sames");
0844        
0845        if(j!=0){
0846      histosdz[j][i]->Draw("sames");
0847        }
0848        
0849        histosdz[j][i]->Write();
0850 
0851        fitResiduals(histosdz[j][i]);   
0852        tmpsdz[j][i] = (TF1*)histosdz[j][i]->GetListOfFunctions()->FindObject("tmp");   
0853        if(tmpsdz[j][i])  {
0854      tmpsdz[j][i]->SetLineColor(colors[j]);
0855      tmpsdz[j][i]->SetLineWidth(2);
0856      tmpsdz[j][i]->Draw("same");  
0857          
0858      meansdz[j][i]=(tmpsdz[j][i])->GetParameter(1);
0859      sigmasdz[j][i]=(tmpsdz[j][i])->GetParameter(2);
0860      
0861      meansdzError[j][i]=(tmpsdz[j][i])->GetParError(1);
0862      sigmasdzError[j][i]=(tmpsdz[j][i])->GetParError(2);
0863      
0864      histoMeansdz[j]->SetBinContent(i+1,meansdz[j][i]); 
0865      histoMeansdz[j]->SetBinError(i+1,meansdzError[j][i]); 
0866      histoMeansdz[j]->GetXaxis()->SetBinLabel(i+1,phipositionString);
0867      histoSigmasdz[j]->SetBinContent(i+1,sigmasdz[j][i]);
0868      histoSigmasdz[j]->SetBinError(i+1,sigmasdzError[j][i]);
0869      histoSigmasdz[j]->GetXaxis()->SetBinLabel(i+1,phipositionString);
0870        }
0871        
0872        else{
0873      histoMeansdz[j]->SetBinContent(i+1,0.); 
0874      histoMeansdz[j]->SetBinError(i+1,0.); 
0875      histoMeansdz[j]->GetXaxis()->SetBinLabel(i+1,phipositionString);
0876      histoSigmasdz[j]->SetBinContent(i+1,0.);
0877      histoSigmasdz[j]->SetBinError(i+1,0.);
0878      histoSigmasdz[j]->GetXaxis()->SetBinLabel(i+1,phipositionString);
0879        }
0880        
0881        c3->Draw("sames");
0882        c3->cd(i+1);
0883        
0884        statObjsdz[j][i] = histosdz[j][i]->GetListOfFunctions()->FindObject("stats");
0885        statsdz[j][i] = static_cast<TPaveStats*>(statObjsdz[j][i]);     
0886        statsdz[j][i]->SetLineColor(colors[j]);
0887        statsdz[j][i]->SetTextColor(colors[j]);
0888        statsdz[j][i]->SetFillColor(10);
0889        statsdz[j][i]->SetShadowColor(10);
0890        statsdz[j][i]->SetX1NDC(0.66);
0891        statsdz[j][i]->SetY1NDC(0.80-j*0.18);
0892        statsdz[j][i]->SetX2NDC(0.99);
0893        statsdz[j][i]->SetY2NDC(0.99-j*0.19);
0894        statsdz[j][i]->Draw("same");
0895        
0896        //################### eta ########################
0897               
0898        char treeleg4Eta[129];
0899 
0900        if(!isNormal_){
0901      if(useBS){
0902        sprintf(treeleg4Eta,"dzBs*10000>>histodzEtafile%i_plot%i",j,i);  
0903      } else {
0904        sprintf(treeleg4Eta,"dzFromMyVertex*10000>>histodzEtafile%i_plot%i",j,i);
0905      } 
0906        } else {
0907      sprintf(treeleg4Eta,"IPLsigFromMyVertex>>histodzEtafile%i_plot%i",j,i);
0908        }      
0909        
0910        c3Eta->cd(i+1);
0911        if(j==0){
0912      trees[j]->Draw(treeleg4Eta,theCut&&theDefaultCut&&theEtaCut);
0913      hmaxDzeta=histosEtadz[0][i]->GetMaximum();    
0914        }
0915        else{ trees[j]->Draw(treeleg4Eta,theCut&&theDefaultCut&&theEtaCut,"sames");
0916      if(histosEtadz[j][i]->GetMaximum() >  hmaxDzeta){
0917        hmaxDzeta = histosdz[j][i]->GetMaximum();
0918      }
0919        }
0920        
0921        histosEtadz[0][i]->GetYaxis()->SetRangeUser(0.01,hmaxDzeta*1.10);
0922        histosEtadz[0][i]->Draw("sames");
0923        
0924        if(j!=0){
0925      histosEtadz[j][i]->Draw("sames");
0926        }
0927        
0928        histosEtadz[j][i]->Write();
0929 
0930        fitResiduals(histosEtadz[j][i]);   
0931        tmpsdzEta[j][i] = (TF1*)histosEtadz[j][i]->GetListOfFunctions()->FindObject("tmp");   
0932        if(tmpsdzEta[j][i])  {
0933      tmpsdzEta[j][i]->SetLineColor(colors[j]);
0934      tmpsdzEta[j][i]->SetLineWidth(2); 
0935      
0936      meansdzEta[j][i]=(tmpsdzEta[j][i])->GetParameter(1);
0937      sigmasdzEta[j][i]=(tmpsdzEta[j][i])->GetParameter(2);
0938      
0939      meansdzEtaError[j][i]=(tmpsdzEta[j][i])->GetParError(1);
0940      sigmasdzEtaError[j][i]=(tmpsdzEta[j][i])->GetParError(2);
0941      
0942      histoMeansdzEta[j]->SetBinContent(i+1,meansdzEta[j][i]); 
0943      histoMeansdzEta[j]->SetBinError(i+1,meansdzEtaError[j][i]); 
0944      histoMeansdzEta[j]->GetXaxis()->SetBinLabel(i+1,etapositionString);
0945      histoSigmasdzEta[j]->SetBinContent(i+1,sigmasdzEta[j][i]);
0946      histoSigmasdzEta[j]->SetBinError(i+1,sigmasdzEtaError[j][i]);
0947      histoSigmasdzEta[j]->GetXaxis()->SetBinLabel(i+1,etapositionString);
0948      
0949        }  else {
0950      histoMeansdzEta[j]->SetBinContent(i+1,0.); 
0951      histoMeansdzEta[j]->SetBinError(i+1,0.); 
0952      histoMeansdzEta[j]->GetXaxis()->SetBinLabel(i+1,etapositionString);
0953      histoSigmasdzEta[j]->SetBinContent(i+1,0.);
0954      histoSigmasdzEta[j]->SetBinError(i+1,0.);
0955      histoSigmasdzEta[j]->GetXaxis()->SetBinLabel(i+1,etapositionString);
0956        }
0957        
0958        c3Eta->Draw("sames");   
0959        c3Eta->cd(i+1);
0960        
0961        statObjsdzEta[j][i] = histosEtadz[j][i]->GetListOfFunctions()->FindObject("stats");
0962        statsdzEta[j][i] = static_cast<TPaveStats*>(statObjsdzEta[j][i]);
0963        statsdzEta[j][i]->SetLineColor(colors[j]);
0964        statsdzEta[j][i]->SetTextColor(colors[j]);
0965        statsdzEta[j][i]->SetShadowColor(10);
0966        statsdzEta[j][i]->SetFillColor(10);
0967        statsdzEta[j][i]->SetX1NDC(0.66);
0968        statsdzEta[j][i]->SetY1NDC(0.80-j*0.18);
0969        statsdzEta[j][i]->SetX2NDC(0.99);
0970        statsdzEta[j][i]->SetY2NDC(0.99-j*0.19);
0971        statsdzEta[j][i]->Draw("same"); 
0972    
0973      }
0974    }
0975    
0976    TString Canvas2Title =_app+"VertexDxyResidualsPhiBin_"+label;
0977    TString Canvas2formatpng=Canvas2Title+".png"; 
0978    TString Canvas2formateps=Canvas2Title+".pdf"; 
0979    c2->SaveAs(Canvas2formatpng);
0980    c2->SaveAs(Canvas2formateps);
0981    c2->Write();
0982   
0983    TString Canvas2TitleEta =_app+"VertexDxyResidualsEtaBin_"+label;
0984    TString Canvas2formatEtapng=Canvas2TitleEta+".png"; 
0985    TString Canvas2formatEtaeps=Canvas2TitleEta+".pdf"; 
0986    c2Eta->SaveAs(Canvas2formatEtapng);
0987    c2Eta->SaveAs(Canvas2formatEtaeps);
0988    c2Eta->Write();
0989   
0990    TString Canvas3Title =_app+"VertexDzResidualsPhiBin_"+label;
0991    TString Canvas3formatpng=Canvas3Title+".png";
0992    TString Canvas3formateps=Canvas3Title+".pdf";  
0993    c3->SaveAs(Canvas3formatpng);
0994    c3->SaveAs(Canvas3formateps);
0995    c3->Write();
0996  
0997    TString Canvas3TitleEta =_app+"VertexDzResidualsEtaBin_"+label;
0998    TString Canvas3formatEtapng=Canvas3TitleEta+".png";
0999    TString Canvas3formatEtaeps=Canvas3TitleEta+".pdf";
1000    c3Eta->SaveAs(Canvas3formatEtapng);
1001    c3Eta->SaveAs(Canvas3formatEtaeps);
1002    c3Eta->Write();
1003    
1004    TPaveText *pt = new TPaveText(0.63,0.73,0.8,0.83,"NDC");
1005    pt->SetFillColor(10);
1006    pt->SetTextColor(1);
1007    //pt->SetTextSize(0.05);
1008    pt->SetTextFont(42);
1009    pt->SetTextAlign(11);
1010    TText *text1 = pt->AddText("Tk Alignment 2012");
1011    text1->SetTextSize(0.05);
1012    TText *text2 = pt->AddText("#sqrt{s}=8 TeV");
1013    text2->SetTextSize(0.05); 
1014    
1015    //######################################################
1016    //  Histograms Means and Widths
1017    //######################################################
1018    
1019    TLegend *legoHisto = new TLegend(0.20,0.65,0.40,0.85); 
1020    legoHisto->SetLineColor(10);
1021    legoHisto->SetTextSize(0.035);
1022    legoHisto->SetTextFont(42);
1023    
1024    TLegend *legoHistoSeparation = new TLegend(0.14,0.67,0.58,0.87); 
1025    legoHistoSeparation-> SetNColumns(2);
1026    legoHistoSeparation->SetLineColor(10);
1027    legoHistoSeparation->SetTextFont(42);
1028    legoHistoSeparation->SetFillColor(10);
1029    legoHistoSeparation->SetTextSize(0.04);
1030    legoHistoSeparation->SetShadowColor(10);
1031    
1032    TLegend *legoHistoSine = new TLegend(0.14,0.16,0.58,0.34); 
1033    legoHistoSine-> SetNColumns(2);
1034    legoHistoSine->SetLineColor(10);
1035    legoHistoSine->SetTextFont(42);
1036    legoHistoSine->SetFillColor(10);
1037    legoHistoSine->SetTextSize(0.03);
1038    legoHistoSine->SetShadowColor(10);
1039 
1040    TF1 *func[nOfFiles];
1041    TF1 *func2[nOfFiles];
1042 
1043    for(Int_t j=0; j < nOfFiles; j++) { 
1044      
1045      //###########################
1046      // Means
1047      //###########################
1048      
1049     //##################### phi ################## 
1050      
1051      c4->cd(1); 
1052      if(j==0){
1053        histoMeansdxy[j]->Draw("e1");
1054      }
1055      else  histoMeansdxy[j]->Draw("e1sames");
1056      
1057      legoHisto->AddEntry(histoMeansdxy[j],LegLabels[j]);
1058      legoHisto->SetFillColor(10);
1059      legoHisto->SetShadowColor(10);
1060      legoHisto->Draw();
1061      pt->Draw("sames");
1062 
1063     c5->cd(1);
1064     if(j==0){
1065       histoMeansdz[j]->Draw("e1");
1066     }
1067     else  histoMeansdz[j]->Draw("e1sames");
1068     legoHisto->Draw();
1069     pt->Draw("sames");
1070 
1071     //#########################################
1072     // 
1073     // Plots with means only
1074     //
1075     //#########################################
1076 
1077     c4monly->cd();
1078 
1079     if(j==0){
1080       histoMeansdxy[j]->SetMarkerSize(1.1);
1081       histoMeansdxy[j]->Draw("e1");
1082     }
1083     else  histoMeansdxy[j]->Draw("e1sames");
1084 
1085     TCanvas *theNewCanvas3 = new TCanvas("NewCanvas3","Fitting Canvas 3",800,600);
1086     TH1F *hnewDxy = (TH1F*)histoMeansdxy[j]->Clone("hnewDxy");
1087     
1088     func[j] = new TF1(Form("myfunc%i",j),splitfunc,_boundMin,_boundMax,4);
1089     func[j]->SetParameters(5,5,5,5);
1090     func[j]->SetParLimits(0,-80.,80.);
1091     func[j]->SetParLimits(1,-80.,80.);
1092     func[j]->SetParLimits(2,-80.,80.);
1093     func[j]->SetParLimits(3,-80.,80.);
1094     func[j]->SetParNames("A","B","k_{+}","k_{-}");
1095     func[j]->SetLineColor(colors[j]);
1096     func[j]->SetNpx(1000);
1097     
1098     func2[j] = new TF1(Form("myfunc2%i",j),simplesine,_boundMin,_boundMax,2);
1099     //func2[j] = new TF1(Form("myfunc2%i",j),myCosine,_boundMin,_boundMax,2);
1100     func2[j]->SetParameters(5,5);
1101     func2[j]->SetParLimits(0,-80.,80.);
1102     func2[j]->SetParLimits(1,-80.,80.);
1103     func2[j]->SetParNames("A","B");
1104     func2[j]->SetLineColor(colors[j]);
1105     func2[j]->SetNpx(1000);
1106      
1107     theNewCanvas3->cd();
1108     hnewDxy->Draw();
1109     hnewDxy->Fit(func2[j]);
1110     
1111     c4monly->cd();
1112     func2[j]->Draw("same");
1113     histoMeansdxy[j]->Draw("e1sames");
1114     
1115     TString COUT1 = Form("#splitline{"+LegLabels[j]+"}{#splitline{d_{x}=%.1f #pm %.1f #mum}{d_{y}=%.1f #pm %.1f #mum}}",func2[j]->GetParameter(0),func2[j]->GetParError(0),func2[j]->GetParameter(1),func2[j]->GetParError(0));
1116     
1117     legoHistoSine->AddEntry(histoMeansdxy[j],COUT1);
1118     legoHistoSine->Draw();
1119 
1120     TPaveText *pi = new TPaveText(0.14,0.71,0.54,0.84,"NDC");
1121     pi->SetFillColor(10);
1122     pi->SetTextColor(1);
1123     pi->SetTextSize(0.04);
1124     pi->SetTextFont(42);
1125     pi->SetTextAlign(11);
1126     TText *textPi1 = pi->AddText("Fitting function:");
1127     textPi1->SetTextSize(0.05);
1128     // TText *textPi2 = pi->AddText("d_{xy} (#varphi) = #ltbar #splitline{|#varphi|<#pi/2     A cos #varphi + B sin #varphi + k_{+} (1- sin #varphi) }{|#varphi|>#pi/2       - A cos #varphi  - B sin #varphi - k_{-} (1- sin#varphi) } ");
1129     //   TText *textPi2 = pi->AddText(Form(" #splitline{d_{xy} (#varphi) = A sin #varphi + B cos #varphi }{A=%.1f  B=%.1f}",func2[j]->GetParameter(0),func2[j]->GetParameter(1)));
1130     TText *textPi2 = pi->AddText("d_{xy} (#varphi) = d_{x} sin #varphi + d_{y} cos #varphi");
1131     textPi2->SetTextSize(0.04); 
1132     textPi2->SetTextSize(0.04); 
1133 
1134     pi->Draw("sames");
1135 
1136     //legoHisto->Draw();
1137     pt->Draw("sames");
1138 
1139     c5monly->cd();
1140     
1141     Double_t deltaZ(0);
1142     Double_t sigmadeltaZ(-1);
1143     
1144     if(j==0){
1145       histoMeansdz[j]->SetMarkerSize(1.1);
1146       histoMeansdz[j]->Draw("e1");
1147     }
1148     else  histoMeansdz[j]->Draw("e1sames");
1149    
1150     TCanvas *theNewCanvas2 = new TCanvas("NewCanvas2","Fitting Canvas 2",800,600);
1151     theNewCanvas2->Divide(2,1);
1152 
1153     TH1F *hnewUp = (TH1F*)histoMeansdz[j]->Clone("hnewUp");
1154     TH1F *hnewDown = (TH1F*)histoMeansdz[j]->Clone("hnewDown");
1155     
1156     fleft[j] = new TF1(Form("fleft_%i",j),fULine,_boundMin,_boundSx,1);
1157     fright[j] = new TF1(Form("fright_%i",j),fULine,_boundDx,_boundMax,1);
1158     fall[j] = new TF1(Form("fall_%i",j),fDLine,_boundSx,_boundDx,1);
1159     
1160     FitULine(hnewUp);  
1161     FitDzUp[j]   = (TF1*)hnewUp->GetListOfFunctions()->FindObject("lineUp"); 
1162     if(FitDzUp[j]){
1163       fleft[j]->SetParameters(FitDzUp[j]->GetParameters());
1164       fleft[j]->SetParErrors(FitDzUp[j]->GetParErrors());
1165       hnewUp->GetListOfFunctions()->Add(fleft[j]);
1166       fright[j]->SetParameters(FitDzUp[j]->GetParameters());
1167       fright[j]->SetParErrors(FitDzUp[j]->GetParErrors());
1168       hnewUp->GetListOfFunctions()->Add(fright[j]);
1169       FitDzUp[j]->Delete();
1170 
1171       theNewCanvas2->cd(1);
1172       MakeNiceTF1Style(fright[j],colors[j]);
1173       MakeNiceTF1Style(fleft[j],colors[j]);
1174       fright[j]->Draw("same");
1175       fleft[j]->Draw("same");
1176     }
1177     
1178     FitDLine(hnewDown);  
1179     FitDzDown[j] = (TF1*)hnewDown->GetListOfFunctions()->FindObject("lineDown");    
1180     
1181     if(FitDzDown[j]){
1182       fall[j]->SetParameters(FitDzDown[j]->GetParameters());
1183       fall[j]->SetParErrors(FitDzDown[j]->GetParErrors());
1184       hnewDown->GetListOfFunctions()->Add(fall[j]);
1185       FitDzDown[j]->Delete();
1186       
1187       theNewCanvas2->cd(2);
1188       MakeNiceTF1Style(fall[j],colors[j]);
1189       fall[j]->Draw("same");
1190       c5monly->cd();
1191       fright[j]->Draw("sames");
1192       fleft[j]->Draw("same");
1193       fall[j]->Draw("same");
1194     }
1195     
1196     if(j==nOfFiles-1){
1197       theNewCanvas2->Close();
1198     }
1199 
1200     deltaZ=(fright[j]->GetParameter(0) - fall[j]->GetParameter(0));
1201     sigmadeltaZ=TMath::Sqrt(fright[j]->GetParError(0)*fright[j]->GetParError(0) + fall[j]->GetParError(0)*fall[j]->GetParError(0));
1202     TString COUT2 = Form("#splitline{"+LegLabels[j]+"}{#Delta z = %.f #pm %.f #mum}",deltaZ,sigmadeltaZ);
1203     
1204     if(j==nOfFiles-1){ 
1205       outfile <<deltaZ<<"|"<<sigmadeltaZ<<endl;
1206     }
1207 
1208     legoHistoSeparation->AddEntry(histoMeansdxy[j],COUT2);
1209     legoHistoSeparation->Draw();
1210 
1211     pt->Draw("sames");
1212 
1213     //##################### Eta ################## 
1214 
1215     c4Eta->cd(1); 
1216     if(j==0){
1217       histoMeansdxyEta[j]->Draw("e1");
1218     }
1219     else  histoMeansdxyEta[j]->Draw("e1sames");
1220     
1221     legoHisto->Draw();
1222     pt->Draw("sames");
1223 
1224     c5Eta->cd(1);
1225     if(j==0){
1226       histoMeansdzEta[j]->Draw("e1");
1227     }
1228     else  histoMeansdzEta[j]->Draw("e1sames");
1229     legoHisto->Draw();
1230     pt->Draw("sames");
1231 
1232     //###########################
1233     // Sigmas
1234     //###########################
1235     
1236     //##################### Phi ################## 
1237 
1238     c4->cd(2);
1239     if(j==0){
1240       histoSigmasdxy[j]->Draw("e1");
1241     }
1242     else histoSigmasdxy[j]->Draw("e1sames");
1243     legoHisto->Draw();
1244     pt->Draw("sames");
1245 
1246     c5->cd(2);
1247     if(j==0){
1248       histoSigmasdz[j]->Draw("e1");
1249     }
1250     else histoSigmasdz[j]->Draw("e1sames");
1251     legoHisto->Draw();
1252     pt->Draw("sames");
1253 
1254     //##################### Eta ##################  
1255 
1256     c4Eta->cd(2);
1257     if(j==0){
1258       histoSigmasdxyEta[j]->Draw("e1");
1259     }
1260     else histoSigmasdxyEta[j]->Draw("e1sames");
1261     legoHisto->Draw();
1262     pt->Draw("sames");
1263 
1264     c5Eta->cd(2);
1265     if(j==0){
1266       histoSigmasdzEta[j]->Draw("e1");
1267     }
1268     else histoSigmasdzEta[j]->Draw("e1sames");
1269     legoHisto->Draw();
1270     pt->Draw("sames");
1271       
1272   }
1273   
1274   TString Canvas4Title =_app+"HistoMeansSigmasDxy_"+label;
1275   TString Canvas4formatpng=Canvas4Title+".png"; 
1276   TString Canvas4formatpdf=Canvas4Title+".pdf"; 
1277   c4->SaveAs(Canvas4formatpng);
1278   c4->SaveAs(Canvas4formatpdf);
1279   c4->Write();
1280   
1281   TString Canvas4Titlemonly =_app+"HistoMeansDxy_"+label;
1282   TString Canvas4formatpngmonly=Canvas4Titlemonly+".png"; 
1283   TString Canvas4formatpdfmonly=Canvas4Titlemonly+".pdf"; 
1284   TString Canvas4formatrootmonly=Canvas4Titlemonly+".root"; 
1285   c4monly->SaveAs(Canvas4formatpngmonly);
1286   c4monly->SaveAs(Canvas4formatpdfmonly);
1287   c4monly->Write();
1288  
1289   TString Canvas5Title =_app+"HistoMeansSigmasDz_"+label;
1290   TString Canvas5formatpng=Canvas5Title+".png"; 
1291   TString Canvas5formatpdf=Canvas5Title+".pdf";
1292   c5->SaveAs(Canvas5formatpng);
1293   c5->Write();
1294   c5->SaveAs(Canvas5formatpdf);
1295 
1296   TString Canvas5Titlemonly =_app+"HistoMeansDzFit_"+label;
1297   TString Canvas5formatpngmonly=Canvas5Titlemonly+".png"; 
1298   TString Canvas5formatpdfmonly=Canvas5Titlemonly+".pdf"; 
1299   TString Canvas5formatrootmonly=Canvas5Titlemonly+".root"; 
1300   c5monly->SaveAs(Canvas5formatpngmonly);
1301   c5monly->SaveAs(Canvas5formatpdfmonly);
1302   c5monly->Write();
1303  
1304   TString Canvas4TitleEta =_app+"HistoMeansSigmasDxyEta_"+label;
1305   TString Canvas4formatEtapng=Canvas4TitleEta+".png"; 
1306   TString Canvas4formatEtapdf=Canvas4TitleEta+".pdf"; 
1307   c4Eta->SaveAs(Canvas4formatEtapng);
1308   c4Eta->SaveAs(Canvas4formatEtapdf);
1309   c4Eta->Write();
1310   
1311   TString Canvas5TitleEta =_app+"HistoMeansSigmasDzEta_"+label;
1312   TString Canvas5formatEtapng=Canvas5TitleEta+".png";
1313   TString Canvas5formatEtapdf=Canvas5TitleEta+".pdf";
1314   c5Eta->SaveAs(Canvas5formatEtapng);
1315   c5Eta->SaveAs(Canvas5formatEtapdf);
1316   c5Eta->Write();
1317   
1318   outfile.close();
1319 
1320 }
1321 
1322 //##########################################
1323 // Fitting Function
1324 //##########################################
1325 
1326 /*--------------------------------------------------------------------*/
1327 void fitResiduals(TH1 *hist)
1328 /*--------------------------------------------------------------------*/
1329 {
1330   //float fitResult(9999);
1331   if (!hist || hist->GetEntries() < 20) return;
1332   
1333   float mean  = hist->GetMean();
1334   float sigma = hist->GetRMS();
1335   
1336   TF1 func("tmp", "gaus", mean - 1.5*sigma, mean + 1.5*sigma); 
1337   if (0 == hist->Fit(&func,"QNR")) { // N: do not blow up file by storing fit!
1338     mean  = func.GetParameter(1);
1339     sigma = func.GetParameter(2);
1340     // second fit: three sigma of first fit around mean of first fit
1341     func.SetRange(mean - 2*sigma, mean + 2*sigma);
1342       // I: integral gives more correct results if binning is too wide
1343       // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1344     if (0 == hist->Fit(&func, "Q0LR")) {
1345       if (hist->GetFunction(func.GetName())) { // Take care that it is later on drawn:
1346     hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
1347       }
1348     }
1349   }
1350   return;
1351 }
1352 
1353 /*--------------------------------------------------------------------*/
1354 void olddoubleGausFitResiduals(TH1 *hist)
1355 /*--------------------------------------------------------------------*/
1356 {
1357   //float fitResult(9999);
1358   // if (!hist || hist->GetEntries() < 20) return;
1359   
1360   float mean  = hist->GetMean();
1361   float sigma = hist->GetRMS();
1362   float interval = (hist->GetXaxis()->GetBinLowEdge(hist->GetNbinsX()+1) - hist->GetXaxis()->GetBinLowEdge(0))/2.;
1363 
1364   Double_t par[6];
1365   TF1 func1("gaus1", "gaus", mean - 1.5*sigma, mean + 1.5*sigma); 
1366   // TF1 func2("gaus2", "gaus", mean - 2.5*sigma, mean + 2.5*sigma); 
1367   TF1 func2("gaus2", "gaus", mean - interval, mean + interval); 
1368   TF1 func("tmp", "gaus(0)+gaus(3)", mean - 3*sigma, mean + 3*sigma); 
1369   
1370   hist->Fit(&func1,"QNR");
1371   hist->Fit(&func2,"QNR+");
1372   func1.GetParameters(&par[0]);
1373   func2.GetParameters(&par[3]);
1374   // cout<<"partials fit done!"<<endl;
1375   
1376   if(hist->GetEntries()>20) {
1377     
1378     // cout<<"histo entries: "<<hist->GetEntries()<<endl; 
1379 
1380     func.SetParameters(par);
1381     func.SetParLimits(1,par[1] + 1*par[2],par[1] + 1*par[2]);
1382     func.SetParLimits(3,par[1] + 1*par[2],par[1] + 1*par[2]);
1383     
1384     if (0 == hist->Fit(&func,"QNR")) { // N: do not blow up file by storing fit!
1385       
1386       // cout<<"before extracting parameters"<<endl;
1387       
1388       mean  = func.GetParameter(1);
1389       sigma = func.GetParameter(5);
1390       
1391       // cout<<"first total fit done!"<<endl;
1392       
1393       // second fit: three sigma of first fit around mean of first fit
1394       
1395       func.SetRange(mean - 3*sigma, mean + 3*sigma);
1396       // I: integral gives more correct results if binning is too wide
1397       // L: Likelihood can treat empty bins correctly (if hist not weighted...)
1398       if (0 == hist->Fit(&func, "Q0LR")) {
1399     if (hist->GetFunction(func.GetName())) { // Take care that it is later on drawn:
1400       hist->GetFunction(func.GetName())->ResetBit(TF1::kNotDraw);
1401     }
1402  
1403     //cout<<"fit done!"<<endl;
1404     
1405       }
1406     }
1407     
1408   } else {
1409     cout<<"Unable to perform double gaussian fit "<<endl;
1410     // func.SetParameters(0);
1411     hist->Fit(&func, "Q0LR");
1412   }
1413   
1414   return;
1415 }
1416 
1417 /*--------------------------------------------------------------------*/
1418 Double_t fDLine(Double_t *x, Double_t *par)
1419 /*--------------------------------------------------------------------*/
1420 {
1421   if (x[0] < _boundSx && x[0] > _boundDx) {
1422     TF1::RejectPoint();
1423     return 0;
1424   }
1425   return par[0];
1426 }
1427 
1428 /*--------------------------------------------------------------------*/
1429 Double_t fULine(Double_t *x, Double_t *par)
1430 /*--------------------------------------------------------------------*/
1431 {
1432   if (x[0] >= _boundSx && x[0] <= _boundDx) {
1433     TF1::RejectPoint();
1434     return 0;
1435   }
1436   return par[0];
1437 }
1438 
1439 /*--------------------------------------------------------------------*/
1440 void FitULine(TH1 *hist)
1441 /*--------------------------------------------------------------------*/
1442 { 
1443   // define fitting function
1444   TF1 func1("lineUp",fULine,_boundMin,_boundMax,1);
1445   //TF1 func1("lineUp","pol0",-0.5,11.5);
1446   
1447   if (0 == hist->Fit(&func1,"QR")) {
1448     if (hist->GetFunction(func1.GetName())) { // Take care that it is later on drawn:
1449       hist->GetFunction(func1.GetName())->ResetBit(TF1::kNotDraw);
1450     }
1451     cout<<"fit Up done!"<<endl;
1452   }
1453   
1454 }
1455 
1456 /*--------------------------------------------------------------------*/
1457 void FitDLine(TH1 *hist)
1458 /*--------------------------------------------------------------------*/
1459 {
1460   // define fitting function
1461   // TF1 func1("lineDown",fDLine,-0.5,11.5,1);
1462   
1463   TF1 func2("lineDown","pol0",_boundSx,_boundDx);
1464   func2.SetRange(_boundSx,_boundDx);
1465   
1466   if (0 == hist->Fit(&func2,"QR")) {
1467     if (hist->GetFunction(func2.GetName())) { // Take care that it is later on drawn:
1468       hist->GetFunction(func2.GetName())->ResetBit(TF1::kNotDraw);
1469     }
1470     cout<<"fit Down done!"<<endl;
1471   } 
1472 }
1473 
1474 /*--------------------------------------------------------------------*/
1475 void MakeNiceTF1Style(TF1 *f1,int color)
1476 /*--------------------------------------------------------------------*/
1477 {
1478   f1->SetLineColor(color);
1479   f1->SetLineWidth(3);
1480   f1->SetLineStyle(2);
1481 }
1482 
1483 /*--------------------------------------------------------------------*/
1484 void MakeNiceHistoStyle(TH1 *hist)
1485 /*--------------------------------------------------------------------*/
1486 {
1487   
1488   //hist->SetTitleSize(0.09); 
1489   hist->GetYaxis()->SetTitleOffset(1.2);
1490   hist->GetXaxis()->SetTitleOffset(0.9);
1491   hist->GetYaxis()->SetTitleSize(0.08);
1492   hist->GetXaxis()->SetTitleSize(0.07);
1493   hist->GetXaxis()->SetTitleFont(42); 
1494   hist->GetYaxis()->SetTitleFont(42);  
1495   hist->GetXaxis()->SetLabelSize(0.05); 
1496   hist->GetYaxis()->SetLabelSize(0.05);
1497   hist->GetXaxis()->SetLabelFont(42);
1498   hist->GetYaxis()->SetLabelFont(42);
1499   hist->SetMarkerSize(1.2);
1500   
1501 }
1502 /*--------------------------------------------------------------------*/
1503 void MakeNiceTrendPlotStyle(TH1 *hist)
1504 /*--------------------------------------------------------------------*/
1505 { 
1506   hist->SetStats(kFALSE);  
1507   hist->SetLineWidth(2);
1508   hist->GetXaxis()->CenterTitle(true);
1509   hist->GetYaxis()->CenterTitle(true);
1510   hist->GetXaxis()->SetTitleFont(42); 
1511   hist->GetYaxis()->SetTitleFont(42);  
1512   hist->GetXaxis()->SetTitleSize(0.06);
1513   hist->GetYaxis()->SetTitleSize(0.06);
1514   hist->GetXaxis()->SetTitleOffset(1.0);
1515   hist->GetYaxis()->SetTitleOffset(1.0);
1516   hist->GetXaxis()->SetLabelFont(42);
1517   hist->GetYaxis()->SetLabelFont(42);
1518   hist->GetYaxis()->SetLabelSize(.06);
1519   hist->GetXaxis()->SetLabelSize(.05);
1520   hist->SetMarkerSize(1.3);
1521 }
1522 
1523 //
1524 // Sinusoidal functions to plot
1525 //
1526 /*--------------------------------------------------------------------*/
1527 Double_t simplesine(Double_t *x, Double_t *par)
1528 /*--------------------------------------------------------------------*/
1529 {
1530   Float_t xx = x[0];
1531   Double_t fitval =-par[1]*cos((TMath::Pi()/_divs)*(xx-_divs))+par[0]*sin((TMath::Pi()/_divs)*(xx-_divs));
1532   return fitval;  
1533 }
1534 
1535 /*--------------------------------------------------------------------*/
1536 Double_t splitfunc(Double_t *x, Double_t *par)
1537 /*--------------------------------------------------------------------*/
1538 {
1539   Double_t xx = x[0];
1540   Double_t dxhs = -par[1]*cos((TMath::Pi()/_divs)*(xx-_divs))+par[0]*sin((TMath::Pi()/_divs)*(xx-_divs)) + par[2]*(1 - sin((TMath::Pi()/_divs)*(xx-_divs)));
1541   Double_t sxhs = -par[1]*cos((TMath::Pi()/_divs)*(xx-_divs))+par[0]*sin((TMath::Pi()/_divs)*(xx-_divs)) - par[3]*(1 - sin((TMath::Pi()/_divs)*(xx-_divs)));
1542                            
1543   if (xx > _boundSx && xx < _boundDx) {
1544     return dxhs;
1545   } else if (xx> _boundMin && xx < _boundSx) {
1546     return sxhs;
1547   } else if (xx == _boundSx && xx == _boundDx) {
1548     return par[0];
1549   } else if (xx > _boundDx && xx < _boundMax)
1550     return sxhs;
1551   else {
1552     std::cout<<"shouldn't ever be the case"<<std::endl;
1553     return 0;
1554   }
1555 }
1556 
1557 /*--------------------------------------------------------------------*/
1558 Double_t splitfunconly(Double_t *x, Double_t *par)
1559 /*--------------------------------------------------------------------*/
1560 {
1561   Double_t xx = x[0];
1562   Double_t dxhs =  par[0]*(1 - sin((TMath::Pi()/_divs)*(xx-_divs)));
1563   Double_t sxhs = -par[1]*(1 - sin((TMath::Pi()/_divs)*(xx-_divs)));
1564                            
1565   if (xx >= 2.50 && xx <= 8.50) {
1566     return dxhs;
1567   } else {
1568     return sxhs;
1569   }
1570 }
1571 
1572 /*--------------------------------------------------------------------*/
1573 Double_t myCosine(Double_t *x, Double_t *par)
1574 /*--------------------------------------------------------------------*/
1575 {
1576   Float_t xx = x[0];
1577   Double_t f = par[0]*cos((TMath::Pi()/_divs)*(xx-_divs) + par[2]);
1578   return f;
1579 }
1580 
1581 /*--------------------------------------------------------------------*/
1582 void doubleGausFitResiduals(TH1F *hist)
1583 /*--------------------------------------------------------------------*/
1584 {
1585 
1586   //TH1F *hist2 = (TH1F*)hist->Clone("h2");
1587   //Use TSpectrum to find the peak candidates
1588   TSpectrum *s = new TSpectrum(2);
1589   Int_t nfound = s->Search(hist,1,"new");
1590   printf("Found %d candidate peaks to fit \n",nfound);
1591   Float_t *xpeaks = s->GetPositionX();
1592   
1593   float mean   = hist->GetMean();
1594   float sigma  = hist->GetRMS();
1595   float x_peak = xpeaks[0];
1596   float fwhm   = calcFWHM(hist);
1597 
1598   std::cout<<" mean: "<<mean 
1599        <<" x_peak: "<<x_peak 
1600        <<" sigma: "<<sigma
1601        <<" fwhm: "<<fwhm<<std::endl; 
1602   
1603   float min1_ = x_peak-fwhm; 
1604   float max1_ = x_peak+fwhm;
1605               
1606   float min2_ = x_peak-2*fwhm;
1607   float max2_ = x_peak+2*fwhm;
1608    
1609   Double_t par[6];
1610   TF1 *g1    = new TF1("g1","gaus",min1_,max1_);
1611   TF1 *g2    = new TF1("g2","gaus",min2_,max2_);
1612   TF1 *func = new TF1("tmp","gaus(0)+gaus(3)",min2_,max2_);
1613   hist->Fit(g1,"QNR");
1614   hist->Fit(g2,"QNR+");
1615 
1616   g1->GetParameters(&par[0]);
1617   g2->GetParameters(&par[3]);
1618   func->SetParameters(par);
1619   // func->FixParameter(1,x_peak);
1620   // func->FixParameter(3,x_peak); 
1621   
1622   hist->Fit(func,"Q0LR+");
1623   
1624   hist->Draw();
1625   func->Draw("same");
1626 }
1627 
1628 
1629 /*--------------------------------------------------------------------*/
1630 Float_t calcFWHM(TH1F *hist)
1631 /*--------------------------------------------------------------------*/
1632 {
1633 
1634   Float_t FWHM(999.);
1635 
1636   Int_t bin_max = hist->GetMaximumBin();
1637   Int_t the_max = hist->GetMaximum();
1638   Int_t halfmax_bin_sx(-999999);
1639   Int_t halfmax_bin_dx(999999);  
1640   
1641   for(Int_t i=1; i<bin_max; i++){
1642     if(hist->GetBinContent(i) > the_max/2. ){
1643       halfmax_bin_sx = i-1;
1644       break;
1645     }  
1646   }
1647   
1648   for(Int_t j=bin_max; j<hist->GetNbinsX(); j++){
1649     //std::cout<<"bin: "<<j<<" bincontent:"<<hist->GetBinContent(j)<<std::endl;
1650     if(hist->GetBinContent(j) < the_max/2. ){
1651       halfmax_bin_dx = j+1;
1652       break;
1653     }  
1654   } 
1655   
1656   // std::cout<<"themax: "<<the_max<<" binmax: "<<bin_max<<" halfmax_bin_sx: "<< halfmax_bin_sx << "  halfmax_bin_dx: "<< halfmax_bin_dx<<std::endl;
1657 
1658   FWHM = hist->GetBinCenter(halfmax_bin_dx) - hist->GetBinCenter(halfmax_bin_sx);
1659   if(FWHM>3000.){
1660     std::cout<<"hist name:"<<hist->GetName()<<std::endl;
1661     std::cout<<"themax: "<<the_max<<" binmax: "<<bin_max<<" halfmax_bin_sx: "<< halfmax_bin_sx << "  halfmax_bin_dx: "<< halfmax_bin_dx<<std::endl;
1662   }
1663   // std::cout<<"FWHM="<<FWHM<<std::endl;
1664   return FWHM;
1665   
1666 }
[ 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!