text: // -*- C++ -*-
  //
  // Package:    PatBTag
  // Class:      PatBTagCommonHistos
  // 
  /**\class PatBTagCommonHistos PatBTagCommonHistos.cc
  
   Description: <Define and Fill common set of histograms depending on flavor>
  
   Implementation:
   
   Create a container of histograms. 
  */
  //
  // Original Author:  J. E. Ramirez
  //
  //
  
  
  // system include files
  #include <memory>
  
  // user include files
  #include "FWCore/ParameterSet/interface/ParameterSet.h"
  #include "FWCore/Utilities/interface/InputTag.h"
  #include "DataFormats/PatCandidates/interface/Jet.h"
  #include "PhysicsTools/PatExamples/interface/PatBTagCommonHistos.h"
  
  #include "FWCore/ServiceRegistry/interface/Service.h"
  #include "CommonTools/UtilAlgos/interface/TFileService.h"
  
  //
  // constructors and destructor
  //
  PatBTagCommonHistos::PatBTagCommonHistos(const edm::ParameterSet& iConfig):
    histocontainer_()
  ,BTagger(iConfig.getParameter< edm::ParameterSet >("BJetOperatingPoints"))
     //now do what ever initialization is needed
  {
    edm::ParameterSet PatBjet_(iConfig.getParameter< edm::ParameterSet >("BjetTag"));
    BTagdisccut_      = PatBjet_.getUntrackedParameter<double>("mindiscriminatorcut",5.0);
    BTagdiscriminator_= PatBjet_.getParameter<std::string>("discriminator");
    BTagpurity_       = PatBjet_.getParameter<std::string>("purity");
  }
  
  
  PatBTagCommonHistos::~PatBTagCommonHistos()
  {
   
     // do anything here that needs to be done at desctruction time
     // (e.g. close files, deallocate resources etc.)
  
  }
  
  
  //
  // member functions
  //
  
  // ------------ method called to for each event  ------------
  void
  PatBTagCommonHistos::Fill( edm::View<pat::Jet>::const_iterator& jet_iter, std::string flavor)
  {
  
  float isb    =jet_iter->bDiscriminator(BTagdiscriminator_);
  //
  // no pt cut histograms
  //
  if(
  	fabs(jet_iter->eta()) < 2.4
  	&& fabs(jet_iter->eta()) > 0
  ){
  	//
  	//Fill histos for Loose(defined in configuration file) cut
  	//tagged using auxiliar function (Loose)
  	//
  	if(BTagger.IsbTag(*jet_iter,BTagpurity_,BTagdiscriminator_)){
  		histocontainer_["jet_pt_uncorr_"+flavor+"_tagged"]->Fill(jet_iter->correctedJet("raw").pt());
  		histocontainer_["jet_pt_"+flavor+"_tagged"]->Fill(jet_iter->pt());
  		histocontainer_["jet_eta_"+flavor+"_tagged"]->Fill(jet_iter->eta());
  	}
  	//
  	// Fill histos 
  	//    tagged minimum discriminant cut criteria
  	//
  	if( isb > float(BTagdisccut_) ){
  		histocontainer_["jet_pt_uncorr_"+flavor+"_taggedmincut"]->Fill(jet_iter->correctedJet("raw").pt());
  		histocontainer_["jet_pt_"+flavor+"_taggedmincut"]->Fill(jet_iter->pt());
  		histocontainer_["jet_eta_"+flavor+"_taggedmincut"]->Fill(jet_iter->eta());
  	}
  	//
  	//fill taggable jets histograms (tracks in jet > 0,1,2) 
          std::map<int,std::string> tagbl;
          tagbl[0]="0";tagbl[1]="1";tagbl[2]="2";
  	for (size_t i=0; i < tagbl.size(); i++)
  		if( jet_iter->associatedTracks().size() > i ){
  			histocontainer_["jet_pt_uncorr_"+flavor+"_taggable"+tagbl[i]]->Fill(jet_iter->correctedJet("raw").pt());
  			histocontainer_["jet_pt_"+flavor+"_taggable"+tagbl[i]]->Fill(jet_iter->pt());
  			histocontainer_["jet_eta_"+flavor+"_taggable"+tagbl[i]]->Fill(jet_iter->eta());
  			if (jet_iter->pt() <  30 ){
  				histocontainer_["jet_eta_"+flavor+"_taggable"+tagbl[i]+"_030"]->Fill(jet_iter->eta());
  			}
  			if (jet_iter->pt() >  30 && jet_iter->pt() < 50){
  				histocontainer_["jet_eta_"+flavor+"_taggable"+tagbl[i]+"_3050"]->Fill(jet_iter->eta());
  			}
  			if (jet_iter->pt() >  50 ){
  				histocontainer_["jet_eta_"+flavor+"_taggable"+tagbl[i]+"_50"]->Fill(jet_iter->eta());
  			}
  		}
  	//
  	// Fill histos needed for normalization
  	// uncorrected pt distributions
  	// corrected pt distributions
  	// eta distributions
  	// scatter plots for pts 
  	histocontainer_["jet_pt_uncorr_"+flavor]->Fill(jet_iter->correctedJet("raw").pt());
  	histocontainer_["jet_pt_"+flavor]->Fill(jet_iter->pt());
  	histocontainer_["jet_eta_"+flavor]->Fill(jet_iter->eta());
  	histocontainer_["tracks_in_jet_"+flavor]->Fill(jet_iter->associatedTracks().size());
  	h2_["jet_scatter_pt_"+flavor]->Fill(jet_iter->pt(),jet_iter->correctedJet("raw").pt());
  	for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  		histocontainer_["pt_tracks_in_jet_"+flavor]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  	}
  	if (jet_iter->pt() <  30 ){
  		histocontainer_["jet_eta_"+flavor+"030"]->Fill(jet_iter->eta());
  		histocontainer_["tracks_in_jet_"+flavor+"030"]->Fill(jet_iter->associatedTracks().size());
  		for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  			histocontainer_["pt_tracks_in_jet_"+flavor+"030"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  		}
  		if (fabs(jet_iter->eta()) <  1.5 ){
  			for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  				histocontainer_["pt_tracks_in_jet_"+flavor+"_center030"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  			}
  		}
  		if (fabs(jet_iter->eta()) >  1.5 && fabs(jet_iter->eta()) <  2.4){
  			for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  				histocontainer_["pt_tracks_in_jet_"+flavor+"_sides030"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  			}
  		}
  	}
  	if (jet_iter->pt() >  30 && jet_iter->pt() < 50){
  		histocontainer_["jet_eta_"+flavor+"3050"]->Fill(jet_iter->eta());
  		histocontainer_["tracks_in_jet_"+flavor+"3050"]->Fill(jet_iter->associatedTracks().size());
  		for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  			histocontainer_["pt_tracks_in_jet_"+flavor+"3050"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  		}
  		if (fabs(jet_iter->eta()) <  1.5 ){
  			for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  				histocontainer_["pt_tracks_in_jet_"+flavor+"_center3050"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  			}
  		}
  		if (fabs(jet_iter->eta()) >  1.5 && fabs(jet_iter->eta()) <  2.4){
  			for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  				histocontainer_["pt_tracks_in_jet_"+flavor+"_sides3050"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  			}
  		}
  	}
  	if (jet_iter->pt() >  50 ){
  		histocontainer_["jet_eta_"+flavor+"50"]->Fill(jet_iter->eta());
  		histocontainer_["tracks_in_jet_"+flavor+"50"]->Fill(jet_iter->associatedTracks().size());
  		for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  			histocontainer_["pt_tracks_in_jet_"+flavor+"50"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  		}
  		if (fabs(jet_iter->eta()) <  1.5 ){
  			for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  				histocontainer_["pt_tracks_in_jet_"+flavor+"_center50"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  			}
  		}
  		if (fabs(jet_iter->eta()) >  1.5 && fabs(jet_iter->eta()) <  2.4){
  			for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  				histocontainer_["pt_tracks_in_jet_"+flavor+"_sides50"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  			}
  		}
  	}
  	if (fabs(jet_iter->eta()) <  1.5 ){
  		for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  			histocontainer_["pt_tracks_in_jet_"+flavor+"_center"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  		}
  	}
  	if (fabs(jet_iter->eta()) >  1.5 && fabs(jet_iter->eta()) <  2.4){
  		for(size_t itrack=0;itrack < jet_iter->associatedTracks().size() ;++itrack){
  			histocontainer_["pt_tracks_in_jet_"+flavor+"_sides"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
  		}
  	}
  
   }//endif
  }//end function
  // ------------ method called once each job just before starting event loop  ------------
  // ------------  This function is needed to set a group of histogram  -------------------
  void 
  PatBTagCommonHistos::Set(std::string flavor)
  {
  
    const int ntptarray = 23;
    const int njptarray = 14;
    const int netaarray = 19;
    Double_t jetetabins[netaarray] = {-2.5,-2.0,-1.75,-1.5,-1.25,-1.0,-0.75,-0.5,-0.25,0.0,0.25,0.5,0.75,1.0,1.25,1.5,1.75,2.0,2.5};
    Double_t jetptxbins[njptarray] = {0.,10.,20.,30., 40., 50., 60., 70., 80, 90., 100., 120., 140., 230.};
    Double_t jetpttbins[ntptarray] = {0.,1.,2.,4.,6.,8.,10.,15.,20.,25.,30., 35.,40., 45., 50., 60., 70., 80, 90., 100., 120., 140., 230.};
    edm::Service<TFileService> fs;
    std::string histoid;
    std::string histotitle;
    
    //Define histograms
    histoid = "jet_pt_"+flavor;        histotitle = flavor+" jet p_{T} [GeV/c]";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),njptarray-1,jetptxbins);
    histoid = "jet_pt_uncorr_"+flavor; histotitle = flavor+" jet uncorrected p_{T} [GeV/c]";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),njptarray-1,jetptxbins);
    histoid = "jet_eta_"+flavor;       
    histocontainer_["jet_eta_"+flavor]=fs->make<TH1D>(histoid.c_str(),"jet #eta",netaarray-1,jetetabins);
    histoid = "jet_eta_"+flavor+"030";
    histocontainer_["jet_eta_"+flavor+"030"]=fs->make<TH1D>(histoid.c_str(),"jet #eta",netaarray-1,jetetabins);
    histoid = "jet_eta_"+flavor+"3050";
    histocontainer_["jet_eta_"+flavor+"3050"]=fs->make<TH1D>(histoid.c_str(),"jet #eta",netaarray-1,jetetabins);
    histoid = "jet_eta_"+flavor+"50";
    histocontainer_["jet_eta_"+flavor+"50"]=fs->make<TH1D>(histoid.c_str(),"jet #eta",netaarray-1,jetetabins);
    histoid = "jet_pt_"+flavor+"_tagged";
    histocontainer_["jet_pt_"+flavor+"_tagged"]=fs->make<TH1D>(histoid.c_str(),"jet_tagged p_{T} [GeV/c]",njptarray-1,jetptxbins);
    histoid = "jet_pt_uncorr_"+flavor+"_tagged";
    histocontainer_["jet_pt_uncorr_"+flavor+"_tagged"]=fs->make<TH1D>(histoid.c_str(),"jet_tagged p_{T} [GeV/c]",njptarray-1,jetptxbins);
    histoid = "jet_eta_"+flavor+"_tagged";
    histocontainer_["jet_eta_"+flavor+"_tagged"]=fs->make<TH1D>(histoid.c_str(),"jet_tagged #eta",netaarray-1,jetetabins);
  
    //tagged min cut
    histoid = "jet_pt_"+flavor+"_taggedmincut";
    histocontainer_["jet_pt_"+flavor+"_taggedmincut"]=fs->make<TH1D>(histoid.c_str(),"jet_tagged p_{T} [GeV/c]",njptarray-1,jetptxbins);
    histoid = "jet_pt_uncorr_"+flavor+"_taggedmincut";
    histocontainer_["jet_pt_uncorr_"+flavor+"_taggedmincut"]=fs->make<TH1D>(histoid.c_str(),"jet_tagged p_{T} [GeV/c]",njptarray-1,jetptxbins);
    histoid = "jet_eta_"+flavor+"_taggedmincut";
    histocontainer_["jet_eta_"+flavor+"_taggedmincut"]=fs->make<TH1D>(histoid.c_str(),"jet_tagged #eta",netaarray-1,jetetabins);
  
    //#tracks per jet
    histoid = "tracks_in_jet_"+flavor;  histotitle = "traks per jet "+flavor;
    histocontainer_["tracks_in_jet_"+flavor]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),31,-0.5,30.5);
    histoid = "tracks_in_jet_"+flavor+"030"; histotitle = "traks per jet "+flavor+ "pt_{T} < 30[GeV/c]";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),31,-0.5,30.5);
    histoid = "tracks_in_jet_"+flavor+"3050"; histotitle = "traks per jet "+flavor+ "30 < pt_{T} < 50[GeV/c]";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),31,-0.5,30.5);
    histoid = "tracks_in_jet_"+flavor+"50"; histotitle = "traks per jet "+flavor+ "pt_{T} > 50[GeV/c]";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),31,-0.5,30.5);
  
    // pt of tracks in bins of jet pt 0-30,30-50,50
    histoid= "pt_tracks_in_jet_"+flavor; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
    histoid= "pt_tracks_in_jet_"+flavor+"030"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
    histoid= "pt_tracks_in_jet_"+flavor+"3050"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
    histoid= "pt_tracks_in_jet_"+flavor+"50"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
  
    // pt of tracks in bins of jet eta abs(eta)<1.5, 1.5<abs(eta)<2.4
    // combined with bins of jet pt
    histoid= "pt_tracks_in_jet_"+flavor+"_center"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
    histoid= "pt_tracks_in_jet_"+flavor+"_center030"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
    histoid= "pt_tracks_in_jet_"+flavor+"_center3050"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
    histoid= "pt_tracks_in_jet_"+flavor+"_center50"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
    histoid= "pt_tracks_in_jet_"+flavor+"_sides"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
    histoid= "pt_tracks_in_jet_"+flavor+"_sides030"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
    histoid= "pt_tracks_in_jet_"+flavor+"_sides3050"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
    histoid= "pt_tracks_in_jet_"+flavor+"_sides50"; histotitle = "track p_{T} [GeV/c] "+ flavor+" jets";
    histocontainer_[histoid]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),ntptarray-1,jetpttbins);
  
    //taggable 0,1,2
    std::map<int,std::string> tagbl;
    tagbl[0]="0";tagbl[1]="1";tagbl[2]="2";
    for (size_t i=0; i < tagbl.size(); i++){
  	histoid = "jet_pt_"+flavor+"_taggable"+tagbl[i];        histotitle = flavor+" jet_taggable p_{T} [GeV/c]";
  	histocontainer_["jet_pt_"+flavor+"_taggable"+tagbl[i]]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),njptarray-1,jetptxbins);
  	histoid = "jet_pt_uncorr_"+flavor+"_taggable"+tagbl[i]; histotitle = flavor+" jet_taggable p_{T} uncorr [GeV/c]";
  	histocontainer_["jet_pt_uncorr_"+flavor+"_taggable"+tagbl[i]]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),njptarray-1,jetptxbins);
  	histoid = "jet_eta_"+flavor+"_taggable"+tagbl[i];       histotitle = flavor+" jet_taggable #eta";
  	histocontainer_["jet_eta_"+flavor+"_taggable"+tagbl[i]]=fs->make<TH1D>(histoid.c_str(),histotitle.c_str(),netaarray-1,jetetabins);
  	histoid = "jet_eta_"+flavor+"_taggable"+tagbl[i]+"_030";
  	histocontainer_["jet_eta_"+flavor+"_taggable"+tagbl[i]+"_030"]=fs->make<TH1D>(histoid.c_str(),"jet #eta",netaarray-1,jetetabins);
  	histoid = "jet_eta_"+flavor+"_taggable"+tagbl[i]+"_3050";
  	histocontainer_["jet_eta_"+flavor+"_taggable"+tagbl[i]+"_3050"]=fs->make<TH1D>(histoid.c_str(),"jet #eta",netaarray-1,jetetabins);
  	histoid = "jet_eta_"+flavor+"_taggable"+tagbl[i]+"_50";
  	histocontainer_["jet_eta_"+flavor+"_taggable"+tagbl[i]+"_50"]=fs->make<TH1D>(histoid.c_str(),"jet #eta",netaarray-1,jetetabins);
    }
  
    histoid = "jet_scatter_pt_"+flavor;
    h2_["jet_scatter_pt_"+flavor]=fs->make<TH2D>(histoid.c_str(),"jet p_{T} uncorr(y) vs corr(x) [GeV/c]",
                                                  njptarray-1,jetptxbins,njptarray-1,jetptxbins);
  
  }//end function Set
  
  // ------------ method called once each job just before starting event loop  ------------
  // ------------              after setting histograms                --------------------
  // ------------  This function is needed to save histogram errors -----------------------
  void 
  PatBTagCommonHistos::Sumw2()
  {
    for (std::map<std::string,TH1D*>::const_iterator ih=histocontainer_.begin();
           ih!= histocontainer_.end(); ++ih) {
        TH1D *htemp = (TH1D*) ih->second;
        htemp->Sumw2();
    }
  }//end function Sumw2
  
  