FirstData/Macros/runProjectQQ.C

#include "ProjectQQ.C"

void BookHistos();
void WriteHistos(Bool_t normalize, Double_t intLumi, Char_t *fNameOut);
//==========================================
void runProjectQQ(Char_t *fNameOut = "histos_jpsi_MC_900GeV_STARTUP.root", //fileName where the histos will be stored
                  Bool_t matchMC = kFALSE,//set to kTRUE for J/psi MC sample and to kFALSE for BG MC and real data
                  Bool_t removeQQ = kFALSE, //for MB/ppMuXLoose samples set to kTRUE; for J/psi MC and real data: set to kFALSE
                  Bool_t normalize = kFALSE, //set to true *only* if you want to normalise the histos; in standard mode: should be kFALSE
                  Double_t intLumi = 16.675//only relevant if "normalise"=kTRUE; 2.36 TeV: 13557/813=16.675; 900 GeV: 6724./404=16.644
                  ){

  ProjectQQ tree;

  BookHistos();

  tree.Loop(removeQQ, matchMC);
  WriteHistos(normalize, intLumi, fNameOut);

}
//==========================================
void BookHistos(){

  //mass
//   Int_t nBinsMass = 20;
  Int_t nBinsMass = 400;
  Double_t massMin = 0.0, massMax = 4.0;
  //pt
  Int_t nBinsPT = 20;
  Double_t pTMin = 0., pTMax = 10.;
  //rap
  Int_t nBinsRap = 60;
  Double_t rapMin = -3.0, rapMax = 3.0;

  double ptbins[] = {0.0, 0.5, 1.0, 1.5,
                     2.0, 2.5, 3.0, 3.5,
                     4.0, 4.5, 5.0, 5.5,
                     6.0, 6.5, 7.0, 7.5,
                     8.0, 8.5, 9.0, 9.5,
                     10.0, 10.5, 11.0, 11.5,
                     12.0, 12.5, 13.0, 13.5,
                     14.0, 14.5, 15.0, 15.5,
                     16.0, 17.0, 18.0, 19.0,
                     20.0, 22.0, 24.0, 26.0,
                     28.0, 30.0, 35.0, 40.0,
                     50.0, 70.0, 100.};
  int nbinspt = sizeof(ptbins)/sizeof(double)-1;


  Char_t name[100], title[300];
  for(int iCharge = 0; iCharge < kNbCharge; iCharge++){
    for(int iCat = 0; iCat < kNbCath; iCat++){
      for (int iCut = 0; iCut < kNbCuts; iCut++) {

        //mass histos
        sprintf(name, "hReco_mass_%s_%s_%s", chargeName[iCharge], oniaCatName[iCat], muCutName[iCut]);
        sprintf(title, ";M [GeV/c^{2}];dN/dM [per %1.0f MeV/c^{2}]", 
                1000.* (massMax-massMin)/nBinsMass);
        hReco_mass[iCharge][iCat][iCut] = new TH1F(name, title, nBinsMass, massMin, massMax);
        hReco_mass[iCharge][iCat][iCut]->Sumw2();
        
        for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
          //pT histos
          sprintf(name, "hReco_pT_%d_%s_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat], muCutName[iCut]);
          if(iSet == 0)
            sprintf(title, ";p_{T}(J/#psi) [GeV/c];dN/dp_{T} [per %1.0f MeV/c]", 1000.* (pTMax-pTMin)/nBinsPT);
          else
            sprintf(title, ";p_{T} [GeV/c];dN/dp_{T} [per %1.0f MeV/c]", 1000.* (pTMax-pTMin)/nBinsPT);
          hReco_pT[iSet][iCharge][iCat][iCut] = new TH1F(name, title, nBinsPT, pTMin, pTMax);
          hReco_pT[iSet][iCharge][iCat][iCut]->Sumw2();
          //rap histos
          sprintf(name, "hReco_rap_%d_%s_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat], muCutName[iCut]);
          if(iSet == 0)
            sprintf(title, ";y(J/#psi);dN/dy");
          else
            sprintf(title, ";y;dN/dy");
          hReco_rap[iSet][iCharge][iCat][iCut] = new TH1F(name, title, nBinsRap, rapMin, rapMax);
          hReco_rap[iSet][iCharge][iCat][iCut]->Sumw2();

          //eta(mu1) vs eta(mu2)
          sprintf(name, "hReco_eta1_eta2_%d_%s_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat], muCutName[iCut]);
          hReco_eta1_eta2[iSet][iCharge][iCat][iCut] = new TH2F(name, ";#eta(#mu^{+});#eta(#mu^{-})", nBinsRap, rapMin, rapMax, nBinsRap, rapMin, rapMax);

          if(iSet == 0){
            //multiplicity
            sprintf(name, "hMult_QQ_%s_%s_%s", chargeName[iCharge], oniaCatName[iCat], muCutName[iCut]);
            hMult_QQ[iCharge][iCat][iCut] = new TH1F(name, ";multiplicity", 10, 0., 10.);
          }
          if(iCut == 0){

            //dimuon vertexing probability
            sprintf(name, "hDimuVtxProb_%d_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat]);
            hDimuVtxProb[iSet][iCharge][iCat] = new TH1F(name, ";#chi^{2} prob. of dimuon vertex", 2000, 0., 1.);
          }
        } //iSet
      } //iCut

        //mass histos: best QQbar
      sprintf(name, "hRecoBest_mass_%s_%s", chargeName[iCharge], oniaCatName[iCat]);
      sprintf(title, ";M [GeV/c^{2}];dN/dM [per %1.0f MeV/c^{2}]", 
              1000.* (massMax-massMin)/nBinsMass);
      hRecoBest_mass[iCharge][iCat] = new TH1F(name, title, nBinsMass, massMin, massMax);
      hRecoBest_mass[iCharge][iCat]->Sumw2();

      for(int iSet = 0; iSet < kNbDimuonSet; iSet++){

        //pT histos: best QQbar
        sprintf(name, "hRecoBest_pT_%d_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat]);
        sprintf(title, ";p_{T}(J/#psi) [GeV/c];dN/dp_{T} [per %1.0f MeV/c]", 
                1000.* (pTMax-pTMin)/nBinsPT);
        hRecoBest_pT[iSet][iCharge][iCat] = new TH1F(name, title, nBinsPT, pTMin, pTMax);
        hRecoBest_pT[iSet][iCharge][iCat]->Sumw2();
        //rap histos: best QQbar
        sprintf(name, "hRecoBest_rap_%d_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat]);
        sprintf(title, ";y(J/#psi);dN/dy");
        hRecoBest_rap[iSet][iCharge][iCat] = new TH1F(name, title, nBinsRap, rapMin, rapMax);
        hRecoBest_rap[iSet][iCharge][iCat]->Sumw2();
      }
    } // iCat
  } // iCharge

  //counter for events where no "bestQQ" could be found:
  sprintf(name, "hRecoNoBest");
  hRecoNoBest = new TH1F(name, ";J/#psi category", kNbCath, 0., kNbCath);

  //single muon histograms:
  for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
    for(int iCharge = 0; iCharge < kNbCharge; iCharge++){
      for(int iCat = 0; iCat < kNbMuCath; iCat++){

        //eta
        sprintf(name, "hMuon_eta_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
        hMuon_eta[iSet][iCharge][iCat] = new TH1F(name, ";#eta(#mu)", nBinsRap, rapMin, rapMax);
        hMuon_eta[iSet][iCharge][iCat]->Sumw2();

        //pT
        sprintf(name, "hMuon_pT_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
        hMuon_pT[iSet][iCharge][iCat] = new TH1F(name, ";p_{T}(#mu) [GeV/c]", nBinsPT, pTMin, pTMax);
        hMuon_pT[iSet][iCharge][iCat]->Sumw2();

        //pT vs eta
        sprintf(name, "hMuon_pTvsEta_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
        hMuon_pT_eta[iSet][iCharge][iCat] = new TH2F(name, ";#eta(#mu);p_{T}(#mu) [GeV/c]", nBinsRap, rapMin, rapMax, nBinsPT, pTMin, pTMax);
        hMuon_pT_eta[iSet][iCharge][iCat]->Sumw2();

        //chi2 of global fit
        if(iCat == 0){//global muons
          sprintf(name, "hChi2GlobalFit_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
          hChi2GlobalFit[iSet][iCharge] = new TH1F(name, ";#chi^{2}/ndf (global fit)", 200, 0, 20.);
        }
        //chi2 of tracker fit
        sprintf(name, "hChi2TrackerFit_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
        hChi2TrackerFit[iSet][iCharge][iCat] = new TH1F(name, ";#chi^{2}/ndf (tracker fit)", 200, 0, 20.);

        //#hits in silicon tracker
        sprintf(name, "hNHitsSilicon_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
        hNHitsSilicon[iSet][iCharge][iCat] = new TH1F(name, ";#hits in silicon tracker", 40, 0, 40.);

        if(iCat == 1){//tracker muons
          //hTM2DCompatibilityTight
          sprintf(name, "hTM2DCompatibilityTight_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
          hTM2DCompatibilityTight[iSet][iCharge] = new TH1F(name, ";TM2DCompatibilityTight", 2, 0, 2);
          //hTMLastStationOptimizedLowPtLoose
          sprintf(name, "hTMLastStationOptimizedLowPtLoose_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
          hTMLastStationOptimizedLowPtLoose[iSet][iCharge] = new TH1F(name, ";TMLastStationOptimizedLowPtLoose", 2, 0, 2);
          //hTrackerMuonArbitrated
          sprintf(name, "hTrackerMuonArbitrated_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
          hTrackerMuonArbitrated[iSet][iCharge] = new TH1F(name, ";TrackerMuonArbitrated", 2, 0, 2);
        }
        //calo compatibility
        sprintf(name, "hCaloCompatibility_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
        hCaloCompatibility[iSet][iCharge][iCat] = new TH1F(name, ";CaloCompatibility", 100, 0., 1.);
      }//cat
    }//charge
  }//set
}

//==========================================
void WriteHistos(Bool_t normalize, Double_t intLumi, Char_t *fNameOut){

  TFile *fOut = new TFile(fNameOut, "RECREATE");
  for(int iCharge = 0; iCharge < kNbCharge; iCharge++){
    for(int iCat = 0; iCat < kNbCath; iCat++){
      for(int iCut = 0; iCut < kNbCuts; iCut++){
        //normalise the histo to 1 nb-1:
        if(normalize){
          hReco_mass[iCharge][iCat][iCut]->Scale(1./intLumi);
          for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
            hReco_pT[iSet][iCharge][iCat][iCut]->Scale(1./intLumi);
            hReco_rap[iSet][iCharge][iCat][iCut]->Scale(1./intLumi);
          }
        }
        hReco_mass[iCharge][iCat][iCut]->Write();
        for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
          hReco_pT[iSet][iCharge][iCat][iCut]->Write();
          hReco_rap[iSet][iCharge][iCat][iCut]->Write();
          if(iCat == 0 || iCat == 2)
            hReco_eta1_eta2[iSet][iCharge][iCat][iCut]->Write();
        }
        hMult_QQ[iCharge][iCat][iCut]->Write();
      }//cut
      
      //normalise the histo to 1 nb-1:
      if(normalize){
        hRecoBest_mass[iCharge][iCat]->Scale(1./intLumi);
        for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
          hRecoBest_pT[iSet][iCharge][iCat]->Scale(1./intLumi);
          hRecoBest_rap[iSet][iCharge][iCat]->Scale(1./intLumi);
        }
      }
      //mass
      hRecoBest_mass[iCharge][iCat]->Write();
      for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
        hRecoBest_pT[iSet][iCharge][iCat]->Write();
        hRecoBest_rap[iSet][iCharge][iCat]->Write();
      }

      for(int iSet = 0; iSet < kNbDimuonSet; iSet++)
        if(iCharge == 0)
          hDimuVtxProb[iSet][iCharge][iCat]->Write();
    }
  }

  hRecoNoBest->Write();


  for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
//     for(int iCharge = 0; iCharge < kNbCharge; iCharge++){
    for(int iCharge = 0; iCharge < 1; iCharge++){//only muons from OS pairs are filled 
      for(int iCat = 0; iCat < kNbMuCath; iCat++){

        hMuon_eta[iSet][iCharge][iCat]->Write();
        hMuon_pT[iSet][iCharge][iCat]->Write();
        hMuon_pT_eta[iSet][iCharge][iCat]->Write();

        if(iCat == 0)
          hChi2GlobalFit[iSet][iCharge]->Write();
        if(iCat == 1){
          hTM2DCompatibilityTight[iSet][iCharge]->Write();
          hTMLastStationOptimizedLowPtLoose[iSet][iCharge]->Write();
          hTrackerMuonArbitrated[iSet][iCharge]->Write();
        }
        if(iCat > 0)
          hChi2TrackerFit[iSet][iCharge][iCat]->Write();

        hNHitsSilicon[iSet][iCharge][iCat]->Write();
        hCaloCompatibility[iSet][iCharge][iCat]->Write();
      }
    }
  }

  fOut->Close();
}
Revision:	1.1
Committed:	Thu Dec 17 18:09:59 2009 UTC (15 years, 4 months ago) by hwoehri
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Log Message:	* empty log message *
#	User	Rev	Content
1	hwoehri	1.1	#include "ProjectQQ.C"
2
3			void BookHistos();
4			void WriteHistos(Bool_t normalize, Double_t intLumi, Char_t *fNameOut);
5			//==========================================
6			void runProjectQQ(Char_t *fNameOut = "histos_jpsi_MC_900GeV_STARTUP.root", //fileName where the histos will be stored
7			Bool_t matchMC = kFALSE,//set to kTRUE for J/psi MC sample and to kFALSE for BG MC and real data
8			Bool_t removeQQ = kFALSE, //for MB/ppMuXLoose samples set to kTRUE; for J/psi MC and real data: set to kFALSE
9			Bool_t normalize = kFALSE, //set to true only if you want to normalise the histos; in standard mode: should be kFALSE
10			Double_t intLumi = 16.675//only relevant if "normalise"=kTRUE; 2.36 TeV: 13557/813=16.675; 900 GeV: 6724./404=16.644
11			){
12
13			ProjectQQ tree;
14
15			BookHistos();
16
17			tree.Loop(removeQQ, matchMC);
18			WriteHistos(normalize, intLumi, fNameOut);
19
20			}
21			//==========================================
22			void BookHistos(){
23
24			//mass
25			// Int_t nBinsMass = 20;
26			Int_t nBinsMass = 400;
27			Double_t massMin = 0.0, massMax = 4.0;
28			//pt
29			Int_t nBinsPT = 20;
30			Double_t pTMin = 0., pTMax = 10.;
31			//rap
32			Int_t nBinsRap = 60;
33			Double_t rapMin = -3.0, rapMax = 3.0;
34
35			double ptbins[] = {0.0, 0.5, 1.0, 1.5,
36			2.0, 2.5, 3.0, 3.5,
37			4.0, 4.5, 5.0, 5.5,
38			6.0, 6.5, 7.0, 7.5,
39			8.0, 8.5, 9.0, 9.5,
40			10.0, 10.5, 11.0, 11.5,
41			12.0, 12.5, 13.0, 13.5,
42			14.0, 14.5, 15.0, 15.5,
43			16.0, 17.0, 18.0, 19.0,
44			20.0, 22.0, 24.0, 26.0,
45			28.0, 30.0, 35.0, 40.0,
46			50.0, 70.0, 100.};
47			int nbinspt = sizeof(ptbins)/sizeof(double)-1;
48
49
50			Char_t name[100], title[300];
51			for(int iCharge = 0; iCharge < kNbCharge; iCharge++){
52			for(int iCat = 0; iCat < kNbCath; iCat++){
53			for (int iCut = 0; iCut < kNbCuts; iCut++) {
54
55			//mass histos
56			sprintf(name, "hReco_mass_%s_%s_%s", chargeName[iCharge], oniaCatName[iCat], muCutName[iCut]);
57			sprintf(title, ";M [GeV/c^{2}];dN/dM [per %1.0f MeV/c^{2}]",
58			1000.* (massMax-massMin)/nBinsMass);
59			hReco_mass[iCharge][iCat][iCut] = new TH1F(name, title, nBinsMass, massMin, massMax);
60			hReco_mass[iCharge][iCat][iCut]->Sumw2();
61
62			for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
63			//pT histos
64			sprintf(name, "hReco_pT_%d_%s_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat], muCutName[iCut]);
65			if(iSet == 0)
66			sprintf(title, ";p_{T}(J/#psi) [GeV/c];dN/dp_{T} [per %1.0f MeV/c]", 1000.* (pTMax-pTMin)/nBinsPT);
67			else
68			sprintf(title, ";p_{T} [GeV/c];dN/dp_{T} [per %1.0f MeV/c]", 1000.* (pTMax-pTMin)/nBinsPT);
69			hReco_pT[iSet][iCharge][iCat][iCut] = new TH1F(name, title, nBinsPT, pTMin, pTMax);
70			hReco_pT[iSet][iCharge][iCat][iCut]->Sumw2();
71			//rap histos
72			sprintf(name, "hReco_rap_%d_%s_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat], muCutName[iCut]);
73			if(iSet == 0)
74			sprintf(title, ";y(J/#psi);dN/dy");
75			else
76			sprintf(title, ";y;dN/dy");
77			hReco_rap[iSet][iCharge][iCat][iCut] = new TH1F(name, title, nBinsRap, rapMin, rapMax);
78			hReco_rap[iSet][iCharge][iCat][iCut]->Sumw2();
79
80			//eta(mu1) vs eta(mu2)
81			sprintf(name, "hReco_eta1_eta2_%d_%s_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat], muCutName[iCut]);
82			hReco_eta1_eta2[iSet][iCharge][iCat][iCut] = new TH2F(name, ";#eta(#mu^{+});#eta(#mu^{-})", nBinsRap, rapMin, rapMax, nBinsRap, rapMin, rapMax);
83
84			if(iSet == 0){
85			//multiplicity
86			sprintf(name, "hMult_QQ_%s_%s_%s", chargeName[iCharge], oniaCatName[iCat], muCutName[iCut]);
87			hMult_QQ[iCharge][iCat][iCut] = new TH1F(name, ";multiplicity", 10, 0., 10.);
88			}
89			if(iCut == 0){
90
91			//dimuon vertexing probability
92			sprintf(name, "hDimuVtxProb_%d_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat]);
93			hDimuVtxProb[iSet][iCharge][iCat] = new TH1F(name, ";#chi^{2} prob. of dimuon vertex", 2000, 0., 1.);
94			}
95			} //iSet
96			} //iCut
97
98			//mass histos: best QQbar
99			sprintf(name, "hRecoBest_mass_%s_%s", chargeName[iCharge], oniaCatName[iCat]);
100			sprintf(title, ";M [GeV/c^{2}];dN/dM [per %1.0f MeV/c^{2}]",
101			1000.* (massMax-massMin)/nBinsMass);
102			hRecoBest_mass[iCharge][iCat] = new TH1F(name, title, nBinsMass, massMin, massMax);
103			hRecoBest_mass[iCharge][iCat]->Sumw2();
104
105			for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
106
107			//pT histos: best QQbar
108			sprintf(name, "hRecoBest_pT_%d_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat]);
109			sprintf(title, ";p_{T}(J/#psi) [GeV/c];dN/dp_{T} [per %1.0f MeV/c]",
110			1000.* (pTMax-pTMin)/nBinsPT);
111			hRecoBest_pT[iSet][iCharge][iCat] = new TH1F(name, title, nBinsPT, pTMin, pTMax);
112			hRecoBest_pT[iSet][iCharge][iCat]->Sumw2();
113			//rap histos: best QQbar
114			sprintf(name, "hRecoBest_rap_%d_%s_%s", iSet, chargeName[iCharge], oniaCatName[iCat]);
115			sprintf(title, ";y(J/#psi);dN/dy");
116			hRecoBest_rap[iSet][iCharge][iCat] = new TH1F(name, title, nBinsRap, rapMin, rapMax);
117			hRecoBest_rap[iSet][iCharge][iCat]->Sumw2();
118			}
119			} // iCat
120			} // iCharge
121
122			//counter for events where no "bestQQ" could be found:
123			sprintf(name, "hRecoNoBest");
124			hRecoNoBest = new TH1F(name, ";J/#psi category", kNbCath, 0., kNbCath);
125
126			//single muon histograms:
127			for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
128			for(int iCharge = 0; iCharge < kNbCharge; iCharge++){
129			for(int iCat = 0; iCat < kNbMuCath; iCat++){
130
131			//eta
132			sprintf(name, "hMuon_eta_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
133			hMuon_eta[iSet][iCharge][iCat] = new TH1F(name, ";#eta(#mu)", nBinsRap, rapMin, rapMax);
134			hMuon_eta[iSet][iCharge][iCat]->Sumw2();
135
136			//pT
137			sprintf(name, "hMuon_pT_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
138			hMuon_pT[iSet][iCharge][iCat] = new TH1F(name, ";p_{T}(#mu) [GeV/c]", nBinsPT, pTMin, pTMax);
139			hMuon_pT[iSet][iCharge][iCat]->Sumw2();
140
141			//pT vs eta
142			sprintf(name, "hMuon_pTvsEta_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
143			hMuon_pT_eta[iSet][iCharge][iCat] = new TH2F(name, ";#eta(#mu);p_{T}(#mu) [GeV/c]", nBinsRap, rapMin, rapMax, nBinsPT, pTMin, pTMax);
144			hMuon_pT_eta[iSet][iCharge][iCat]->Sumw2();
145
146			//chi2 of global fit
147			if(iCat == 0){//global muons
148			sprintf(name, "hChi2GlobalFit_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
149			hChi2GlobalFit[iSet][iCharge] = new TH1F(name, ";#chi^{2}/ndf (global fit)", 200, 0, 20.);
150			}
151			//chi2 of tracker fit
152			sprintf(name, "hChi2TrackerFit_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
153			hChi2TrackerFit[iSet][iCharge][iCat] = new TH1F(name, ";#chi^{2}/ndf (tracker fit)", 200, 0, 20.);
154
155			//#hits in silicon tracker
156			sprintf(name, "hNHitsSilicon_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
157			hNHitsSilicon[iSet][iCharge][iCat] = new TH1F(name, ";#hits in silicon tracker", 40, 0, 40.);
158
159			if(iCat == 1){//tracker muons
160			//hTM2DCompatibilityTight
161			sprintf(name, "hTM2DCompatibilityTight_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
162			hTM2DCompatibilityTight[iSet][iCharge] = new TH1F(name, ";TM2DCompatibilityTight", 2, 0, 2);
163			//hTMLastStationOptimizedLowPtLoose
164			sprintf(name, "hTMLastStationOptimizedLowPtLoose_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
165			hTMLastStationOptimizedLowPtLoose[iSet][iCharge] = new TH1F(name, ";TMLastStationOptimizedLowPtLoose", 2, 0, 2);
166			//hTrackerMuonArbitrated
167			sprintf(name, "hTrackerMuonArbitrated_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
168			hTrackerMuonArbitrated[iSet][iCharge] = new TH1F(name, ";TrackerMuonArbitrated", 2, 0, 2);
169			}
170			//calo compatibility
171			sprintf(name, "hCaloCompatibility_%d_%s_%s", iSet, chargeName[iCharge], muCatName[iCat]);
172			hCaloCompatibility[iSet][iCharge][iCat] = new TH1F(name, ";CaloCompatibility", 100, 0., 1.);
173			}//cat
174			}//charge
175			}//set
176			}
177
178			//==========================================
179			void WriteHistos(Bool_t normalize, Double_t intLumi, Char_t *fNameOut){
180
181			TFile *fOut = new TFile(fNameOut, "RECREATE");
182			for(int iCharge = 0; iCharge < kNbCharge; iCharge++){
183			for(int iCat = 0; iCat < kNbCath; iCat++){
184			for(int iCut = 0; iCut < kNbCuts; iCut++){
185			//normalise the histo to 1 nb-1:
186			if(normalize){
187			hReco_mass[iCharge][iCat][iCut]->Scale(1./intLumi);
188			for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
189			hReco_pT[iSet][iCharge][iCat][iCut]->Scale(1./intLumi);
190			hReco_rap[iSet][iCharge][iCat][iCut]->Scale(1./intLumi);
191			}
192			}
193			hReco_mass[iCharge][iCat][iCut]->Write();
194			for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
195			hReco_pT[iSet][iCharge][iCat][iCut]->Write();
196			hReco_rap[iSet][iCharge][iCat][iCut]->Write();
197			if(iCat == 0 \|\| iCat == 2)
198			hReco_eta1_eta2[iSet][iCharge][iCat][iCut]->Write();
199			}
200			hMult_QQ[iCharge][iCat][iCut]->Write();
201			}//cut
202
203			//normalise the histo to 1 nb-1:
204			if(normalize){
205			hRecoBest_mass[iCharge][iCat]->Scale(1./intLumi);
206			for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
207			hRecoBest_pT[iSet][iCharge][iCat]->Scale(1./intLumi);
208			hRecoBest_rap[iSet][iCharge][iCat]->Scale(1./intLumi);
209			}
210			}
211			//mass
212			hRecoBest_mass[iCharge][iCat]->Write();
213			for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
214			hRecoBest_pT[iSet][iCharge][iCat]->Write();
215			hRecoBest_rap[iSet][iCharge][iCat]->Write();
216			}
217
218			for(int iSet = 0; iSet < kNbDimuonSet; iSet++)
219			if(iCharge == 0)
220			hDimuVtxProb[iSet][iCharge][iCat]->Write();
221			}
222			}
223
224			hRecoNoBest->Write();
225
226
227			for(int iSet = 0; iSet < kNbDimuonSet; iSet++){
228			// for(int iCharge = 0; iCharge < kNbCharge; iCharge++){
229			for(int iCharge = 0; iCharge < 1; iCharge++){//only muons from OS pairs are filled
230			for(int iCat = 0; iCat < kNbMuCath; iCat++){
231
232			hMuon_eta[iSet][iCharge][iCat]->Write();
233			hMuon_pT[iSet][iCharge][iCat]->Write();
234			hMuon_pT_eta[iSet][iCharge][iCat]->Write();
235
236			if(iCat == 0)
237			hChi2GlobalFit[iSet][iCharge]->Write();
238			if(iCat == 1){
239			hTM2DCompatibilityTight[iSet][iCharge]->Write();
240			hTMLastStationOptimizedLowPtLoose[iSet][iCharge]->Write();
241			hTrackerMuonArbitrated[iSet][iCharge]->Write();
242			}
243			if(iCat > 0)
244			hChi2TrackerFit[iSet][iCharge][iCat]->Write();
245
246			hNHitsSilicon[iSet][iCharge][iCat]->Write();
247			hCaloCompatibility[iSet][iCharge][iCat]->Write();
248			}
249			}
250			}
251
252			fOut->Close();
253			}