root/BxAnalysis/NtupleChecker.C

0001 //-----------------------------------------------------
0002 //  Authors: Lorenzo Uplegger  : uplegger@cern.ch
0003 //           Sushil s. Chauhan : sushil@fnal.gov
0004 //
0005 //   Last modified: lm 07 July 2011
0006 //------------------------------------------------------
0007 // This script read the pv  ntuples from Beam Spot
0008 // area and make plots for each bx using pv Fitter.
0009 // It also create  a canvas with all the bx plotted for
0010 // each beam spot variable in the output file.
0011 //-----------------------------------------------------
0012
0013 #include "RecoVertex/BeamSpotProducer/interface/BeamSpotTreeData.h"
0014 #include "RecoVertex/BeamSpotProducer/interface/FcnBeamSpotFitPV.h"
0015 #include "BSFitData.h"
0016
0017 #include "TFitterMinuit.h"
0018 #include "Minuit2/FCNBase.h"
0019 #include <TTree.h>
0020 #include <TH1F.h>
0021 #include <TH2F.h>
0022 #include <TFile.h>
0023 #include <TString.h>
0024 #include <TList.h>
0025 #include <TSystemFile.h>
0026 #include <TSystemDirectory.h>
0027 #include <TCanvas.h>
0028 #include <TLegend.h>
0029 #include <iostream>
0030 #include <fstream>
0031
0032 using namespace std;
0033
0034 //--------Global definiton of some funtions
0035 bool runPVFitter(std::map< int, std::vector<BeamSpotFitPVData> > bxMap_, int NFits, int tmpLumiCounter);
0036 void PlotHistoBX(TString OFN, TString ODN);
0037 void DefineHistStyle(TH1F *h1, int bx);
0038 void FillHist(TH1F* h, int fitN, double pos, double posError, TString lab);
0039 void FillnPVInfo(TH1F* h2,int fitN, float npv, TString lab );
0040 void PlotAllBunches(TCanvas* myCanvas,TH1F* hist[],int bunchN, std::map<int,int> tmpMap_);
0041
0042
0043 std::map<int, std::map< int, std::vector<BeamSpotFitPVData> > > StoreRunLSPVdata_;
0044 std::map<int, map<int, std::vector<BSFitData> > > FitDone_Results_;
0045 std::map<int, TString > FitDone_RunLumiRange_;
0046 std::map<int, map<int, float> > bxMapPVSize_;
0047
0048 ofstream outdata;
0049
0050
0051
0052 //-----Here is the main funtion---------------------------
0053 void NtupleChecker(){
0054
0055   //----------------------------------------------//
0056   //                Input parameters              //
0057   //----------------------------------------------//
0058   Int_t beginRunNumber = -1;           //give -1 if do not want to set lower run limit
0059   Int_t endRunNumber   = -1;               //give -1 if do not want to set upper run limit
0060
0061   Int_t beginLSNumber  = -1;
0062   Int_t endLSNumber    = -1;
0063
0064   Int_t FitNLumi       = 100;
0065
0066   // TString OutPutFileName ="BxAnalysis_Fill_1718.root";
0067   //TString OutPutDir      ="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_BSPOT/burkett/BxNtuples/";
0068
0069
0070   // TString OutPutFileName ="BxAnalysis_Fill_1815_test.root";
0071   // TString OutPutDir      ="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_GLOBAL/malgeri/cmssw/CMSSW_4_2_2/src/RecoVertex/BeamSpotProducer/scripts/root/BxAnalysis/mytestntuple";
0072   //TString OutPutFileName ="BxAnalysis_Fill_1815.root";
0073   //TString OutPutDir      ="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_GLOBAL/malgeri/cmssw/CMSSW_4_2_2/src/RecoVertex/BeamSpotProducer/scripts/root/BxAnalysis/";
0074
0075   TString OutPutFileName ="BxAnalysis_Fill_1901_JetPD.root";
0076   TString OutPutDir      ="/afs/cern.ch/cms/CAF/CMSCOMM/COMM_BSPOT/burkett/BxNtuples/Fill_1901_JetPD/";
0077
0078
0079   //---------------------------------------------//
0080
0081
0082
0083   //save few things to ouput log file
0084   outdata.open("LogFile.dat");
0085   cout<<"-----------Step - 1 : Storing Info from Root file to a vector-----------------------"<<endl;
0086   outdata<<" -----------Step - 1 : Storing Info from Root file to a vector-----------------------"<<endl;
0087
0088   //clear the map before storing
0089   StoreRunLSPVdata_.clear();
0090   FitDone_RunLumiRange_.clear();
0091   bxMapPVSize_.clear();
0092   FitDone_Results_.clear();
0093
0094   //set direcotry structure
0095   //  TString path = "/afs/cern.ch/cms/CAF/CMSCOMM/COMM_GLOBAL/malgeri/cmssw/CMSSW_4_2_2/src/RecoVertex/BeamSpotProducer/scripts/root/BxAnalysis/mytestntuple/";
0096
0097   TString path = "/afs/cern.ch/cms/CAF/CMSCOMM/COMM_BSPOT/burkett/BxNtuples/Fill_1901_JetPD/";
0098
0099   TSystemDirectory sourceDir("hi",path);
0100   TList* fileList = sourceDir.GetListOfFiles();
0101   TIter next(fileList);
0102   TSystemFile* fileName;
0103
0104   //define file numbers
0105   int fileNumber = 1;
0106   int maxFiles = -1;
0107   BeamSpotTreeData aData;
0108
0109
0110   //-------------------Store all input in a map-------------------------
0111   //check filename
0112   while ((fileName = (TSystemFile*)next()) && fileNumber >  maxFiles){
0113
0114     if(TString(fileName->GetName()) == "." || TString(fileName->GetName()) == ".."  ){
0115       continue;
0116     }
0117
0118
0119     //set tree
0120     TTree* aTree = 0;
0121     TFile file(path+fileName->GetName(),"READ");//STARTUP
0122     cout << "Opening file: " << path+fileName->GetName() << endl;
0123     file.cd();
0124     aTree = (TTree*)file.Get("PrimaryVertices");
0125
0126     cout << (100*fileNumber)/(fileList->GetSize()-2) << "% of files done." << endl;
0127     ++fileNumber;
0128     if(aTree == 0){
0129       cout << "Can't find the tree" << endl;
0130       continue;
0131     }//check if file is not found
0132
0133
0134     aData.setBranchAddress(aTree);
0135     //ok start reading the tree one b; one
0136     for(unsigned int entry=0; entry<aTree->GetEntries(); entry++){
0137       aTree->GetEntry(entry);
0138
0139       //for all runs and all LS
0140       if( beginRunNumber == -1 && endRunNumber == -1 ){
0141         if(beginLSNumber== -1 && endLSNumber == -1){
0142       BeamSpotFitPVData Input=aData.getPvData();
0143       StoreRunLSPVdata_[aData.getRun()][aData.getLumi()].push_back(Input);
0144         }}
0145
0146
0147       //for Selected Runs and Slected or All LS
0148       if((aData.getRun()>= beginRunNumber && aData.getRun()<= endRunNumber) ||  (aData.getRun()>= beginRunNumber && endRunNumber == -1) || (beginRunNumber==-1 && aData.getRun()<= endRunNumber) ){
0149         if((aData.getLumi()<= beginLSNumber && aData.getLumi()>= endLSNumber) || (beginLSNumber ==-1 && endLSNumber ==-1) ){
0150       BeamSpotFitPVData Input=aData.getPvData();
0151       StoreRunLSPVdata_[aData.getRun()][aData.getLumi()].push_back(Input);
0152         }}
0153
0154     }//loop over entries of tree
0155
0156     file.Close();
0157
0158   }//loop over all the files
0159
0160   //All the runs
0161   for( std::map<int, std::map< int, std::vector<BeamSpotFitPVData> > >::iterator p=StoreRunLSPVdata_.begin(); p!=StoreRunLSPVdata_.end(); ++p)
0162     {  for( std::map< int, std::vector<BeamSpotFitPVData> >::iterator q=p->second.begin(); q!=p->second.end(); ++q)
0163         {  cout<<" Run Number= "<<p->first<<"   LS Number= "<<q->first<<endl;
0164       outdata<<" Run Number= "<<p->first<<"   LS Number= "<<q->first<<endl;
0165         }
0166     }
0167
0168
0169
0170   if(StoreRunLSPVdata_.size()==0)outdata<<" EIther the file is  missing or  it does not contain any data!!!!! Please check  "<<endl;
0171
0172   cout<<"-----------Step - 2 : Now Running the PVFitter for each Bunch Xrossing-----------------------"<<endl;
0173   outdata<<"-----------Step - 2 : Now Running the PVFitter for each Bunch Xrossing-----------------------"<<endl;
0174
0175   //-----------------Run the fitter and store the reulsts for plotting ----------------------------
0176
0177   Int_t Fit_Done = 0;
0178   std::vector<BeamSpotFitPVData> pvStore_;
0179   std::map< int, std::vector<BeamSpotFitPVData> > bxMap_;    //this map for each bx for N lumi
0180   bxMap_.clear();
0181
0182   int Lumi_lo, Lumi_up,RunNumber;
0183   bool pv_fit_Ok;
0184   Int_t LastRun=0;
0185   Int_t LastLumi,LumiCounter,tmpLastLumi;
0186   LumiCounter=0;
0187
0188
0189   for( std::map<int, std::map< int, std::vector<BeamSpotFitPVData> > >::iterator RunIt=StoreRunLSPVdata_.begin(); RunIt!=StoreRunLSPVdata_.end(); ++RunIt)
0190     {
0191       if(StoreRunLSPVdata_.size()==0) continue;
0192
0193       if(LastRun==0)LastRun=RunIt->first;
0194       LumiCounter=0;
0195       RunNumber=RunIt->first;
0196
0197       std::map< int, std::vector<BeamSpotFitPVData> >::iterator LastLumiIt;
0198
0199       for( std::map< int, std::vector<BeamSpotFitPVData> >::iterator LumiIt=RunIt->second.begin(); LumiIt!=RunIt->second.end(); ++LumiIt)
0200     {
0201           LastLumiIt = RunIt->second.end();
0202           LastLumiIt--;
0203
0204           if(LumiCounter==0)Lumi_lo=LumiIt->first;
0205           LumiCounter++;
0206
0207           bool RemainingLSFit=false;
0208
0209           if((LumiCounter == FitNLumi)  || (LumiIt->first == LastLumiIt->first  && LumiCounter != FitNLumi && LumiCounter >0) ){
0210         RemainingLSFit=true;
0211         Lumi_up=LumiIt->first;}
0212
0213
0214       //now loop over pv vectors hence: ->Run->LS->PV/Bx
0215       for(size_t pvIt=0; pvIt < (LumiIt->second).size(); pvIt++)
0216         {
0217           int bx = (LumiIt->second)[pvIt].bunchCrossing;
0218           bxMap_[bx].push_back((LumiIt->second)[pvIt]);
0219           //if(RunNumber ==140123 && LumiIt->first==93 ) cout<<"Has-------->  bx = "<<bx<<endl;
0220
0221         }
0222
0223       if((RemainingLSFit) && (LumiCounter != 0))
0224         {    RemainingLSFit=false;
0225               //cout<<"Run Number ="<<RunNumber<<",    Lumi Range = "<<Lumi_lo<<" - "<<Lumi_up<<endl;
0226               int tmpLumiCounter=0;
0227               tmpLumiCounter = LumiCounter;
0228
0229               LumiCounter=0;
0230               Fit_Done++;
0231               if(runPVFitter(bxMap_,Fit_Done,tmpLumiCounter)){
0232
0233                 //store the run : LS range as Tstring
0234                 Char_t RunLSRange[300];
0235                 sprintf(RunLSRange,"%d%s%d%s%d%s",RunNumber,":",Lumi_lo," - ",Lumi_up,"\0");
0236                 TString RunLSRange_(RunLSRange);
0237                 FitDone_RunLumiRange_[Fit_Done]=RunLSRange_;
0238
0239           }
0240
0241           bxMap_.clear();
0242
0243         }//check if Desired LS are collected to perform the fit
0244
0245
0246     }//loop over Lumi/pvdata map
0247
0248
0249     } //loop over Run map
0250
0251
0252   cout<<"-----------Step - 3 : Now Filling the histograms for each Bunch Xrossing-----------------------"<<endl;
0253   outdata<<"-----------Step - 3 : Now Filling the histograms for each Bunch Xrossing-----------------------"<<endl;
0254
0255   //-------------------------------Plot the histograms and store them in a root file----------------
0256
0257   PlotHistoBX(OutPutDir, OutPutFileName );
0258
0259   outdata.close();
0260 }//NtupleChecker ends here
0261
0262
0263
0264
0265 //--------------------Plot histo after the fit results-----------------------------
0266 void  PlotHistoBX(TString ODN, TString OFN){
0267
0268   TFile f1(ODN+"/"+OFN,"recreate");
0269
0270
0271   // std::map<int, TString>::iterator RLS;
0272   // RLS=FitDone_RunLumiRange_.begin();
0273
0274   //     fit#, bx, bsfit
0275
0276   // std::map<int, map<int, std::vector<BSFitData> > >::iterator tmpIt=FitDone_Results_.begin();
0277   //  std::map<int, std::vector<BSFitData> >::iterator tmpBxIt;
0278   //  tmpBxIt = (tmpIt->second).begin();
0279
0280
0281   Int_t PointsToPlot, bunchN;
0282   bunchN=0;
0283   PointsToPlot = FitDone_Results_.size();
0284   //cout<<"Total Fit Points ="<<PointsToPlot<<endl;
0285
0286   //get the number of bunches
0287
0288   // lm you can have different set of bunches in the different fit (i.e. run/LS)
0289   // for example in fill 1815 you get that fit#1 has 1044 bunches, while fit #16 has 1043 and some of them are not in common!
0290
0291
0292   //create a map to store index of histo and bx N
0293
0294   std::map<int, int> bxNIndexMap_;
0295   bxNIndexMap_.clear();
0296
0297   for(std::map<int, map<int, std::vector<BSFitData> > >::iterator tmpIt=FitDone_Results_.begin();tmpIt!=FitDone_Results_.end();++tmpIt){
0298
0299     if((tmpIt->second).size()==0) continue;
0300
0301     for( std::map<int, std::vector<BSFitData> >::iterator BxIt=tmpIt->second.begin(); BxIt!=tmpIt->second.end(); ++BxIt){
0302       int this_bx=BxIt->first;
0303       if (bxNIndexMap_.count(this_bx)==0)
0304     {
0305       // new bx found! add id to the list
0306       bxNIndexMap_[this_bx]=bunchN;
0307       bunchN++;
0308     }
0309
0310     } // close loop on bx
0311   } // close loop on fits
0312
0313     cout << "Found " << bunchN << " unique bx numbers, creating all histos" << endl;
0314     outdata << "Found "<< bunchN << " unique bx numbers, creating all histos" << endl;
0315
0316
0317
0318
0319   TH1F* h_X_bx_[bunchN];
0320   TH1F* h_Y_bx_[bunchN];
0321   TH1F* h_Z_bx_[bunchN];
0322   TH1F* h_widthX_bx_[bunchN];
0323   TH1F* h_widthY_bx_[bunchN];
0324   TH1F* h_widthZ_bx_[bunchN];
0325   //
0326   TH1F* h_dxdz_bx_[bunchN];
0327   TH1F* h_dydz_bx_[bunchN];
0328   //
0329   TH1F* h_nPV_bx_[bunchN];
0330
0331
0332   //define the histo
0333   for (std::map<int, int>::iterator bxit=bxNIndexMap_.begin() ; bxit != bxNIndexMap_.end(); bxit++ ){
0334
0335     int x= bxit->second;
0336
0337     Char_t XName[254];
0338     sprintf(XName,"%s%d%s","h_X_bx_",bxit->first,"\0");
0339     TString XName_(XName);
0340     h_X_bx_[x] =new TH1F(XName_,XName_,PointsToPlot,0.,PointsToPlot);
0341     h_X_bx_[x]->GetYaxis()->SetTitle("BS Fit X  (cm)");
0342     DefineHistStyle(h_X_bx_[x],x);
0343
0344     Char_t YName[254];
0345     sprintf(YName,"%s%d%s","h_Y_bx_",bxit->first,"\0");
0346     TString YName_(YName);
0347     h_Y_bx_[x] =new TH1F(YName_,YName_,PointsToPlot,0.,PointsToPlot);
0348     h_Y_bx_[x]->GetYaxis()->SetTitle("BS Fit Y  (cm)");
0349     DefineHistStyle(h_Y_bx_[x],x);
0350
0351     Char_t ZName[254];
0352     sprintf(ZName,"%s%d%s","h_Z_bx_",bxit->first,"\0");
0353     TString ZName_(ZName);
0354     h_Z_bx_[x] =new TH1F(ZName_,ZName_,PointsToPlot,0.,PointsToPlot);
0355     h_Z_bx_[x]->GetYaxis()->SetTitle("BS Fit Z  (cm)");
0356     DefineHistStyle(h_Z_bx_[x],x);
0357
0358     Char_t WXName[254];
0359     sprintf(WXName,"%s%d%s","h_widthX_bx_",bxit->first,"\0");
0360     TString WXName_(WXName);
0361     h_widthX_bx_[x] =new TH1F(WXName_,WXName_,PointsToPlot,0.,PointsToPlot);
0362     h_widthX_bx_[x]->GetYaxis()->SetTitle("BS Fit #sigma X  (cm)");
0363     DefineHistStyle(h_widthX_bx_[x],x);
0364
0365     Char_t WYName[254];
0366     sprintf(WYName,"%s%d%s","h_widthY_bx_",bxit->first,"\0");
0367     TString WYName_(WYName);
0368     h_widthY_bx_[x] =new TH1F(WYName_,WYName_,PointsToPlot,0.,PointsToPlot);
0369     h_widthY_bx_[x]->GetYaxis()->SetTitle("BS Fit #sigma Y  (cm)");
0370     DefineHistStyle(h_widthY_bx_[x],x);
0371
0372
0373     Char_t WZName[254];
0374     sprintf(WZName,"%s%d%s","h_widthZ_bx_",bxit->first,"\0");
0375     TString WZName_(WZName);
0376     h_widthZ_bx_[x] =new TH1F(WZName_,WZName_,PointsToPlot,0.,PointsToPlot);
0377     h_widthZ_bx_[x]->GetYaxis()->SetTitle("BS Fit #sigma Z  (cm)");
0378     DefineHistStyle(h_widthZ_bx_[x],x);
0379
0380
0381     Char_t dxdzName[254];
0382     sprintf(dxdzName,"%s%d%s","h_dxdz_bx_",bxit->first,"\0");
0383     TString dxdzName_(dxdzName);
0384     h_dxdz_bx_[x] =new TH1F(dxdzName_,dxdzName_,PointsToPlot,0.,PointsToPlot);
0385     h_dxdz_bx_[x]->GetYaxis()->SetTitle("BS Fit dxdz  (cm)");
0386     DefineHistStyle(h_dxdz_bx_[x],x);
0387
0388     Char_t dydzName[254];
0389     sprintf(dydzName,"%s%d%s","h_dydz_bx_",bxit->first,"\0");
0390     TString dydzName_(dydzName);
0391     h_dydz_bx_[x] =new TH1F(dydzName_,dydzName_,PointsToPlot,0.,PointsToPlot);
0392     h_dydz_bx_[x]->GetYaxis()->SetTitle("BS Fit dydz  (cm)");
0393     DefineHistStyle(h_dydz_bx_[x],x);
0394
0395     //histo for PV # for each bx
0396
0397     Char_t PVSize[254];
0398     sprintf(PVSize,"%s%d%s","h_nPV_bx_",bxit->first,"\0");
0399     TString PVSize_(PVSize);
0400     h_nPV_bx_[x] =new TH1F(PVSize_,PVSize_,PointsToPlot,0.,PointsToPlot);
0401     h_nPV_bx_[x]->GetYaxis()->SetTitle("# of Primary Vertices / LS");
0402     DefineHistStyle(h_nPV_bx_[x],x);
0403
0404   }
0405
0406
0407
0408
0409
0410   for( std::map<int, map<int, std::vector<BSFitData> > >::iterator FitIt=FitDone_Results_.begin(); FitIt!=FitDone_Results_.end(); ++FitIt){
0411
0412     if(FitDone_Results_.size()==0) continue;
0413
0414     std::map<int,float> bxnpv=bxMapPVSize_[FitIt->first];
0415
0416     for( std::map<int, std::vector<BSFitData> >::iterator BxIt=FitIt->second.begin(); BxIt!=FitIt->second.end(); ++BxIt){
0417
0418       for(size_t bsfit=0; bsfit < (BxIt->second).size(); bsfit++){
0419
0420     TString RLSstring=FitDone_RunLumiRange_[FitIt->first];
0421
0422     cout<<"Fit #  = "<<FitIt->first<<"  Run:LS1-LS2  ="<<(RLSstring)<<"       Bx Number ="<< BxIt->first<<"     X0 = "<<((BxIt->second)[bsfit].xyz[0])<<endl;
0423     outdata<<"Fit #  = "<<FitIt->first<<"  Run:LS1-LS2  ="<<(RLSstring)<<"       Bx Number ="<< BxIt->first<<"     X0 = "<<((BxIt->second)[bsfit].xyz[0])<<endl;
0424
0425     int bxindex=bxNIndexMap_[BxIt->first];
0426     FillHist(h_X_bx_[bxindex],FitIt->first,(BxIt->second)[bsfit].xyz[0],(BxIt->second)[bsfit].xyzErr[0],RLSstring);
0427     FillHist(h_Y_bx_[bxindex],FitIt->first,(BxIt->second)[bsfit].xyz[1],(BxIt->second)[bsfit].xyzErr[1],RLSstring);
0428     FillHist(h_Z_bx_[bxindex],FitIt->first,(BxIt->second)[bsfit].xyz[2],(BxIt->second)[bsfit].xyzErr[2],RLSstring);
0429
0430     FillHist(h_widthX_bx_[bxindex],FitIt->first,(BxIt->second)[bsfit].xyzwidth[0],(BxIt->second)[bsfit].xyzwidthErr[0],RLSstring);
0431     FillHist(h_widthY_bx_[bxindex],FitIt->first,(BxIt->second)[bsfit].xyzwidth[1],(BxIt->second)[bsfit].xyzwidthErr[1],RLSstring);
0432     FillHist(h_widthZ_bx_[bxindex],FitIt->first,(BxIt->second)[bsfit].xyzwidth[2],(BxIt->second)[bsfit].xyzwidthErr[2],RLSstring);
0433
0434     FillHist(h_dxdz_bx_[bxindex],FitIt->first,(BxIt->second)[bsfit].dxdz,(BxIt->second)[bsfit].dxdzErr,RLSstring);
0435     FillHist(h_dydz_bx_[bxindex],FitIt->first,(BxIt->second)[bsfit].dydz,(BxIt->second)[bsfit].dydzErr,RLSstring);
0436
0437     //Fill PV Info
0438     FillnPVInfo(h_nPV_bx_[bxindex],FitIt->first, bxnpv[BxIt->first],RLSstring);
0439
0440       } //loop over position errors
0441     }//Loop over bx
0442   }//Loop over each fit
0443
0444
0445   TDirectory *X0 = f1.mkdir("X0");
0446   TDirectory *Y0 = f1.mkdir("Y0");
0447   TDirectory *Z0 = f1.mkdir("Z0");
0448   TDirectory *SigmaZ0 = f1.mkdir("Sigma_Z0");
0449   TDirectory *SigmaX0 = f1.mkdir("width_X0");
0450   TDirectory *SigmaY0 = f1.mkdir("Width_Y0");
0451
0452   TDirectory *dxdz = f1.mkdir("dxdz");
0453   TDirectory *dydz = f1.mkdir("dydz");
0454
0455   TDirectory *nPV     = f1.mkdir("bx_nPV");
0456
0457
0458
0459   for(int t=0; t< bunchN; t++){
0460     X0->cd();
0461     h_X_bx_[t]->Write();
0462     Y0->cd();
0463     h_Y_bx_[t]->Write();
0464     Z0->cd();
0465     h_Z_bx_[t]->Write();
0466     SigmaX0->cd();
0467     h_widthX_bx_[t]->Write();
0468     SigmaY0->cd();
0469     h_widthY_bx_[t]->Write();
0470     SigmaZ0->cd();
0471     h_widthZ_bx_[t]->Write();
0472
0473     dxdz->cd();
0474     h_dxdz_bx_[t]->Write();
0475
0476     dydz->cd();
0477     h_dydz_bx_[t]->Write();
0478
0479     nPV->cd();
0480     h_nPV_bx_[t]->Write();
0481
0482   }
0483
0484   //Plot all bx on the same canvas and save in root file
0485   TCanvas *All_X= new TCanvas("All_X","",7,8,699,499);
0486   PlotAllBunches(All_X, h_X_bx_, bunchN, bxNIndexMap_);
0487   X0->cd();
0488   All_X->Write();
0489   TCanvas *All_Y= new TCanvas("All_Y","",7,8,699,499);
0490   PlotAllBunches(All_Y, h_Y_bx_, bunchN, bxNIndexMap_);
0491   Y0->cd();
0492   All_Y->Write();
0493   TCanvas *All_Z= new TCanvas("All_Z","",7,8,699,499);
0494   PlotAllBunches(All_Z, h_Z_bx_, bunchN, bxNIndexMap_);
0495   Z0->cd();
0496   All_Z->Write();
0497   TCanvas *All_widthX= new TCanvas("All_widthX","",7,8,699,499);
0498   PlotAllBunches(All_widthX, h_widthX_bx_, bunchN, bxNIndexMap_);
0499   SigmaX0->cd();
0500   All_widthX->Write();
0501   TCanvas *All_widthY= new TCanvas("All_widthY","",7,8,699,499);
0502   PlotAllBunches(All_widthY, h_widthY_bx_, bunchN, bxNIndexMap_);
0503   SigmaY0->cd();
0504   All_widthY->Write();
0505   TCanvas *All_widthZ= new TCanvas("All_widthZ","",7,8,699,499);
0506   PlotAllBunches(All_widthZ, h_widthZ_bx_, bunchN, bxNIndexMap_);
0507   SigmaZ0->cd();
0508   All_widthZ->Write();
0509   TCanvas *All_dxdz= new TCanvas("All_dxdz","",7,8,699,499);
0510   PlotAllBunches(All_dxdz, h_dxdz_bx_, bunchN, bxNIndexMap_);
0511   dxdz->cd();
0512   All_dxdz->Write();
0513   TCanvas *All_dydz= new TCanvas("All_dydz","",7,8,699,499);
0514   PlotAllBunches(All_dydz, h_dydz_bx_, bunchN, bxNIndexMap_);
0515   dydz->cd();
0516   All_dydz->Write();
0517
0518   TCanvas *All_nPV= new TCanvas("All_nPV","",7,8,699,499);
0519   PlotAllBunches(All_nPV,h_nPV_bx_, bunchN, bxNIndexMap_);
0520   nPV->cd();
0521   All_nPV->Write();
0522
0523
0524   cout<<"The PV fit is performed for all the  "<<bunchN<<"  bunches"<<endl;
0525   outdata<<"The PV fit is performed for all the  "<<bunchN<<"  bunches"<<endl;
0526
0527   //  if(bunchN> 1000 || bunchN < 0)cout<<"Something Went Wrong OR there is no input to the fit!!!! "<<endl;
0528   //  if(bunchN> 1000 || bunchN < 0)outdata<<"Something Went Wrong OR there is no input to the fit!!!! "<<endl;
0529   if(bunchN> 4000 || bunchN < 0)cout<<"Something Went Wrong OR there is no input to the fit!!!! "<<endl;
0530   if(bunchN> 4000 || bunchN < 0)outdata<<"Something Went Wrong OR there is no input to the fit!!!! "<<endl;
0531
0532   f1.cd();
0533   //write the root file
0534   f1.Close();
0535
0536   //clear these vectors
0537   bxNIndexMap_.clear();
0538   FitDone_Results_.clear();
0539   FitDone_RunLumiRange_.clear();
0540   bxMapPVSize_.clear();
0541 }//PlotHistoBX ends here
0542
0543
0544 //---------Define Hstogram Styles-----
0545 void DefineHistStyle(TH1F *h1, int bx){
0546   bx=bx+1;
0547   // if(bx>8) bx=bx-8;
0548
0549   h1->SetFillColor(0);
0550   h1->SetStats(0);
0551   h1->GetYaxis()->SetTitleOffset(1.20);
0552   h1->SetOption("e1");
0553   h1->GetYaxis()->CenterTitle();
0554   h1->SetTitleOffset(1.2);
0555   h1->SetLineWidth(2);
0556
0557   h1->SetMarkerStyle((20+bx));
0558   if(bx > 10)h1->SetMarkerStyle((20+bx-10));
0559   h1->SetMarkerSize(0.9);
0560   h1->SetMarkerColor(bx);
0561   if(bx > 9)h1->SetMarkerColor(bx-9);
0562   h1->SetLineColor(bx);
0563   if(bx > 9)h1->SetLineColor(bx-9);
0564   //remove yellow
0565   if(bx==5){h1->SetMarkerColor(bx*9);
0566     h1->SetLineColor(bx*9);}
0567
0568   h1->SetLineStyle(1);
0569   if(bx>10)h1->SetLineStyle(bx-9);
0570
0571
0572   bx=0;
0573
0574 }
0575 //------------------------
0576
0577
0578 void FillHist(TH1F* h, int fitN, double pos, double posError, TString lab){
0579
0580   h->SetBinContent(fitN,pos);
0581   h->SetBinError(fitN,posError);
0582   h->GetXaxis()->SetBinLabel(fitN,lab);
0583
0584 }
0585
0586
0587 void FillnPVInfo(TH1F* h2,int fitN, float npv, TString lab ){
0588
0589   h2->SetBinContent(fitN,npv);
0590   //h->SetBinError(fitN,posError);
0591   h2->GetXaxis()->SetBinLabel(fitN,lab);
0592
0593   npv=0.0;
0594
0595 }
0596
0597
0598
0599
0600 //--------------put all bx on same canvas
0601
0602 void PlotAllBunches(TCanvas* myCanvas,TH1F* hist[], int bunchN, std::map<int,int> tmpMap_){
0603
0604
0605   TString legT;
0606   legT = "";
0607   TLegend *leg = new TLegend(0.72,0.72,0.90,0.90, legT);
0608   leg->SetFillColor(0);
0609   myCanvas->SetFillColor(0);
0610
0611   std::map<int,int>::const_iterator iter=tmpMap_.begin();
0612
0613   for(int i=0;i<bunchN;i++){
0614
0615     Char_t nbx[50];
0616     sprintf(nbx,"%s%d%s"," bx # ",(iter->first),"\0");
0617     TString nbx_(nbx);
0618
0619     leg->AddEntry(hist[i],nbx_,"PL");
0620
0621     if(i==0)hist[i]->Draw();
0622     else{hist[i]->Draw("SAMES");}
0623     iter++;
0624
0625   }
0626
0627   leg->Draw();
0628
0629 }
0630
0631
0632
0633
0634 //---------------------------------------------------------------------
0635 //------------------Here we define the Fitting module------------------
0636 //---------------------------------------------------------------------
0637 bool runPVFitter(std::map< int, std::vector<BeamSpotFitPVData> > bxMap_, int NFits, int tmpLumiCounter){
0638
0639   float errorScale_    = 0.9;
0640   float sigmaCut_      = 5.0;
0641   bool fit_ok          = true;
0642   Int_t minNrVertices_ = 10;
0643
0644
0645   for ( std::map<int,std::vector<BeamSpotFitPVData> >::const_iterator pvStore = bxMap_.begin();
0646     pvStore!=bxMap_.end(); ++pvStore) {
0647
0648
0649     //fill number of pv for each bx crossing:
0650     bxMapPVSize_[NFits][pvStore->first]=(((Float_t)(pvStore->second).size()/(Float_t)tmpLumiCounter));
0651
0652     //cout<<"  For bx ="<<pvStore->first<<"     PV # =  "<<pvStore->second.size()<<endl;
0653     if ( (pvStore->second).size() <= minNrVertices_) {
0654       cout << " Not enough PVs, Setting to zero ->"<<(pvStore->second).size() << std::endl;
0655       fit_ok = false;
0656       //continue;
0657     }
0658
0659
0660     //LL function and fitter
0661     FcnBeamSpotFitPV* fcn = new FcnBeamSpotFitPV(pvStore->second);
0662     TFitterMinuit minuitx;
0663     minuitx.SetMinuitFCN(fcn);
0664     // fit parameters: positions, widths, x-y correlations, tilts in xz and yz
0665     minuitx.SetParameter(0,"x",0.,0.02,-10.,10.);
0666     minuitx.SetParameter(1,"y",0.,0.02,-10.,10.);
0667     minuitx.SetParameter(2,"z",0.,0.20,-30.,30.);
0668     minuitx.SetParameter(3,"ex",0.015,0.01,0.,10.);
0669     minuitx.SetParameter(4,"corrxy",0.,0.02,-1.,1.);
0670     minuitx.SetParameter(5,"ey",0.015,0.01,0.,10.);
0671     minuitx.SetParameter(6,"dxdz",0.,0.0002,-0.1,0.1);
0672     minuitx.SetParameter(7,"dydz",0.,0.0002,-0.1,0.1);
0673     minuitx.SetParameter(8,"ez",1.,0.1,0.,30.);
0674     minuitx.SetParameter(9,"scale",errorScale_,errorScale_/10.,errorScale_/2.,errorScale_*2.);
0675
0676     // first iteration without correlations
0677
0678     int ierr=0;
0679     minuitx.FixParameter(4);
0680     minuitx.FixParameter(6);
0681     minuitx.FixParameter(7);
0682     minuitx.FixParameter(9);
0683     minuitx.SetMaxIterations(100);
0684
0685     minuitx.SetPrintLevel(0);
0686     minuitx.CreateMinimizer();
0687     ierr = minuitx.Minimize();
0688
0689     if ( ierr ) {
0690       cout << "3D beam spot fit failed in 1st iteration" <<endl;
0691       fit_ok =false;
0692       //continue;
0693     }
0694
0695
0696     // refit with harder selection on vertices
0697
0698     fcn->setLimits(minuitx.GetParameter(0)-sigmaCut_*minuitx.GetParameter(3),
0699            minuitx.GetParameter(0)+sigmaCut_*minuitx.GetParameter(3),
0700            minuitx.GetParameter(1)-sigmaCut_*minuitx.GetParameter(5),
0701            minuitx.GetParameter(1)+sigmaCut_*minuitx.GetParameter(5),
0702            minuitx.GetParameter(2)-sigmaCut_*minuitx.GetParameter(8),
0703            minuitx.GetParameter(2)+sigmaCut_*minuitx.GetParameter(8));
0704     ierr = minuitx.Minimize();
0705     if (ierr) {
0706       cout << "3D beam spot fit failed in 2nd iteration" << endl;
0707       fit_ok = false;
0708       //continue;
0709     }
0710
0711     // refit with correlations
0712     minuitx.ReleaseParameter(4);
0713     minuitx.ReleaseParameter(6);
0714     minuitx.ReleaseParameter(7);
0715     ierr = minuitx.Minimize();
0716
0717     if ( ierr ) {
0718       cout << "3D beam spot fit failed in 3rd iteration: Setting to zero" <<endl;
0719       fit_ok = false;
0720       //continue;
0721     }
0722
0723     BSFitData PVFitDataForNLumi;
0724
0725     PVFitDataForNLumi.xyz[0] = minuitx.GetParameter(0);
0726     PVFitDataForNLumi.xyz[1] = minuitx.GetParameter(1);
0727     PVFitDataForNLumi.xyz[2] = minuitx.GetParameter(2);
0728     PVFitDataForNLumi.xyzErr[0] = minuitx.GetParError(0);
0729     PVFitDataForNLumi.xyzErr[1] = minuitx.GetParError(1);
0730     PVFitDataForNLumi.xyzErr[2] = minuitx.GetParError(2);
0731
0732     PVFitDataForNLumi.xyzwidth[0] = minuitx.GetParameter(3);
0733     PVFitDataForNLumi.xyzwidth[1] = minuitx.GetParameter(5);
0734     PVFitDataForNLumi.xyzwidth[2] = minuitx.GetParameter(8);
0735
0736     PVFitDataForNLumi.xyzwidthErr[0] = minuitx.GetParError(3);
0737     PVFitDataForNLumi.xyzwidthErr[1]= minuitx.GetParError(5);
0738     PVFitDataForNLumi.xyzwidthErr[2] = minuitx.GetParError(8);
0739
0740     PVFitDataForNLumi.dxdz = minuitx.GetParameter(6);
0741     PVFitDataForNLumi.dydz = minuitx.GetParameter(7);
0742     PVFitDataForNLumi.dxdzErr = minuitx.GetParError(6);
0743     PVFitDataForNLumi.dydzErr = minuitx.GetParError(7);
0744
0745
0746
0747
0748     if(!fit_ok){//set zero if fit fails//or not enough PV
0749       PVFitDataForNLumi.xyz[0] = 0.0;
0750       PVFitDataForNLumi.xyz[1] = 0.0;
0751       PVFitDataForNLumi.xyz[2] = 0.0;
0752       PVFitDataForNLumi.xyzErr[0] = 0.0;
0753       PVFitDataForNLumi.xyzErr[1] = 0.0;
0754       PVFitDataForNLumi.xyzErr[2] = 0.0;
0755
0756       PVFitDataForNLumi.xyzwidth[0] = 0.0;
0757       PVFitDataForNLumi.xyzwidth[1] = 0.0;
0758       PVFitDataForNLumi.xyzwidth[2] = 0.0;
0759
0760       PVFitDataForNLumi.xyzwidthErr[0] = 0.0;
0761       PVFitDataForNLumi.xyzwidthErr[1] = 0.0;
0762       PVFitDataForNLumi.xyzwidthErr[2] = 0.0;
0763
0764       PVFitDataForNLumi.dxdz = 0.0;
0765       PVFitDataForNLumi.dydz = 0.0;
0766       PVFitDataForNLumi.dxdzErr = 0.0;
0767       PVFitDataForNLumi.dydzErr = 0.0;
0768     }
0769
0770     fit_ok=true; //now safe to put it ture
0771     /*
0772       cout<<" Fitted X="<< minuitx.GetParameter(0)<<"    Fitted WidthX ="<<minuitx.GetParameter(3)<<endl;
0773       cout<<" Fitted Y="<< minuitx.GetParameter(1)<<"    Fitted WidthY ="<<minuitx.GetParameter(5)<<endl;
0774       cout<<" Fitted Z="<< minuitx.GetParameter(2)<<"    Fitted WidthZ ="<<minuitx.GetParameter(8)<<endl;
0775     */
0776
0777
0778     //store beam spot fit result as a function of bx
0779     FitDone_Results_[NFits][pvStore->first].push_back(PVFitDataForNLumi);
0780
0781   }//loop over map
0782
0783 }//endof runBXFiiter