#include #include #include "TPostScript.h" TCanvas* c1; TPostScript* myps; int const nPtHat = 15; int const rootFN = 1; TFile* f[rootFN]; TString titleArray[9]; int globalK; TH1F* Pthat0; TFile* fileCreate ; TString lam; TH2F* Sum2DHisto(TString hname){ TH2F* TwoDH[nPtHat]; for (int i =0; i < nPtHat; ++i) TwoDH[i] = new TH2F(); //Just for initailizing the bin and title for the new histo TString tStr = "rootFiles/" + lam + "/" + "Lambda" + lam + "Pt" + "018" + ".root"; TFile* tmpFile = new TFile(tStr); TH2F* h2 = (TH2F*)tmpFile->Get(hname); TString title = h2->GetTitle(); int nbinsX = h2->GetNbinsX(); //pthat int nbinsY = h2->GetNbinsY(); // mass bin int const nbinX = nbinsX; int const nbinY = nbinsY; //Initialize bins for the merge histogram double MjjBin[nbinY +1]; double lowEdge =0; //Initialize bins for merged histogram for (int i = 1; i < nbinY+2; ++i){ lowEdge = h2->GetYaxis()->GetBinLowEdge(i); MjjBin[i-1] = lowEdge; }//i TH2F* histoF = new TH2F("histoF", title, 15, 0., 15., nbinY, MjjBin); for (int i =0; i < nbinX; ++i){ for (int j =0; j < nbinY; ++j){ float value =0; float error =0; int counter = 0; //File Starts here for (int k =0; k < rootFN; ++k){ char buffer[20]; sprintf(buffer, "%d%d%d", i,j,k); TwoDH[counter] = (TH2F*)f[counter]->Get(hname)->Clone(hname + buffer); float binCt = TwoDH[counter]->GetBinContent(i+1, j+1); float binError = TwoDH[counter]->GetBinError(i+1, j+1); value = value + binCt; error = error + (binError*binError); ++counter; }//k float errorF = TMath::Sqrt(error); histoF->SetBinContent(i+1,j+1,value); histoF->SetBinError(i+1,j+1,errorF); }//j }//i delete [] TwoDH; tmpFile = 0; delete tmpFile; return histoF; }//Sum2DHisto TH1F* Sum1DHisto(TString hname){ TH1F* OneDH[nPtHat]; for (int m =0; m < nPtHat; ++m) OneDH[m] = new TH1F(); //Just for initailizing the bin and title for the new histo TString tStr = "rootFiles/" + lam + "/" + "Lambda" + lam + "Pt" + "018" + ".root"; TFile* tmpFile = new TFile(tStr); TH1F* h1 = (TH1F*)tmpFile->Get(hname); TString title = h1->GetTitle(); int nbinsX = h1->GetNbinsX(); //pthat int const nbinX = nbinsX; TH1F* histoF = new TH1F("histoF", title, 2000, 0., 2000); for (int i =0; i < nbinX; ++i){ float value =0.; float error =0.; int counter = 0; //File Starts here for (int k =0; k < rootFN; ++k){ char buffer[20]; sprintf(buffer, "%d%d", i,k); OneDH[counter] = (TH1F*)f[counter]->Get(hname)->Clone(hname + buffer); float binCt = OneDH[counter]->GetBinContent(i+1); float binError = OneDH[counter]->GetBinError(i+1); value = value + binCt; error = error + (binError * binError); ++counter; }//k histoF->SetBinContent(i+1,value); float errorF = TMath::Sqrt(error); histoF->SetBinError(i+1,errorF); if (i % 500 == 0) cout << "Done bin no. " << i << endl; }//i for (int m =0; m < nPtHat; ++m){ OneDH[m] = 0; delete OneDH[m] ; } delete h1; return histoF; }//Sum1DHisto TString getTitle(TString hName){ TString tit1; int len = hName.Length(); hName.Remove(len - 1); char* chrStr = hName.Data(); if ( (!strcmp(chrStr,"MjjL5VPthat2")) ) tit1 = titleArray[0]; else if ( (!strcmp(chrStr,"MjjL5VPthat0")) ) tit1 = titleArray[6]; else if ( (!strcmp(chrStr,"MjjL5VPthat1")) ) tit1 = titleArray[7]; else if ( (!strcmp(chrStr,"MjjL5VPthat2")) ) tit1 = titleArray[8]; else if ( (!strcmp(chrStr,"MjjL5VPthat3")) ) tit1 = titleArray[0] + titleArray[1]; else if ( (!strcmp(chrStr,"MjjL5VPthat4")) ) tit1 = titleArray[0] + titleArray[1] + titleArray[2]; else if ( (!strcmp(chrStr,"MjjL5VPthat5")) ) tit1 = titleArray[0] + titleArray[1] + titleArray[3]; else if ( (!strcmp(chrStr,"MjjL5VPthat6")) ) tit1 = titleArray[0] + titleArray[4]; else if ( (!strcmp(chrStr,"MjjL5VPthat7")) ) tit1 = titleArray[0] + titleArray[4] + titleArray[2]; else if ( (!strcmp(chrStr,"MjjL5VPthat8")) ) tit1 = titleArray[0] + titleArray[4] + titleArray[3]; else if ( (!strcmp(chrStr,"MjjL5VPthat9")) ) tit1 = titleArray[0] + titleArray[5]; else if ( (!strcmp(chrStr,"MjjL5VPthat10")) ) tit1 = titleArray[0] + titleArray[5] + titleArray[2]; else if ( (!strcmp(chrStr,"MjjL5VPthat11")) ) tit1 = titleArray[0] + titleArray[5] + titleArray[3]; else tit1 = ""; return tit1 ; } TH1F* Addition(TString hName){ TH2F* htemp = Sum2DHisto(hName); int nbinsX = htemp->GetNbinsX(); //pthat int nbinsY = htemp->GetNbinsY(); // mass bin int const nbinX = nbinsX; int const nbinY = nbinsY; double xsec[nbinX]; //cross section from pythia in each pthat bin int nev[nbinX]; // No. of events in each pthat for (int i =0; i < nbinX; ++i){ xsec[i] = 0. ; nev[i] = 0.; } //Initialize bins for the merge histogram double MjjBin[nbinY +1]; double lowEdge =0; //Initialize bins for merged histogram for (int i = 1; i < nbinY+2; ++i){ lowEdge = htemp->GetYaxis()->GetBinLowEdge(i); MjjBin[i-1] = lowEdge; }//i //Initialize histogram for each pthat TH1F* ptHatHisto[nbinX]; //Make title char *p = hName.Data(); TString tit0; if ( (!strcmp(p,"MjjL5VPthat2h")) || (!strcmp(p,"MjjL5VPthat0h")) || (!strcmp(p,"MjjL5VPthat1h")) || (!strcmp(p,"MjjL5VPthat2h")) || (!strcmp(p,"MjjL5VPthat3h")) || (!strcmp(p,"MjjL5VPthat4h")) || (!strcmp(p,"MjjL5VPthat5h")) || (!strcmp(p,"MjjL5VPthat6h")) || (!strcmp(p,"MjjL5VPthat7h")) || (!strcmp(p,"MjjL5VPthat8h")) || (!strcmp(p,"MjjL5VPthat9h")) || (!strcmp(p,"MjjL5VPthat10h")) || (!strcmp(p,"MjjL5VPthat11h")) ) tit0 = "Had Level: "; else tit0 = "Parton Level: "; TString tit1 = getTitle(hName); TString title = tit0 + tit1; cout << " hName, Title = " << hName << " " << title << endl; for (int i = 0; i < nbinX; ++i){ char buff[20]; sprintf(buff, "%d", i); TString pthat = TString(buff); TString hTitle = title + ",_Pthat" + pthat; TString histoName = "_Pthat" + pthat; ptHatHisto[i] = new TH1F(histoName, hTitle, nbinY, MjjBin); ptHatHisto[i]->GetXaxis()->SetTitle("Dijet Mass (in GeV)"); }//i //SetBinContent for ptHatHisto from 2D Histo for (int i =0; i < nbinX; ++i){ for (int j =0; j < nbinY; ++j){ double binContent = htemp->GetBinContent(i+1, j+1); double binError = htemp->GetBinError(i+1,j+1); // cout << "binContent = " << binContent << endl; ptHatHisto[i]->SetBinContent(j+1, binContent); ptHatHisto[i]->SetBinError(j+1, binError); }//j }//i myps->NewPage(); c1->Clear(); c1->Divide(3,4); for (int i =0; i < nbinX; ++i){ int canvasCnt = i%12 + 1; c1->cd(canvasCnt); ptHatHisto[i]->SetMarkerStyle(29); ptHatHisto[i]->SetMarkerColor(2); ptHatHisto[i]->SetStats(kFALSE); ptHatHisto[i]->Draw(); if (canvasCnt%12 == 0){ c1->Update(); myps->NewPage(); } }//i c1->Update(); //Draw only once if (globalK == 0){ Pthat0 = Sum1DHisto("Pthat0"); c1->cd(0); c1->SetLogy(); Pthat0->GetXaxis()->SetTitle("Dijet Mass (in GeV)"); Pthat0->GetYaxis()->SetTitle("No. of events"); Pthat0->SetMarkerStyle(29); Pthat0->SetMarkerColor(2); Pthat0->SetStats(kFALSE); Pthat0->Draw("P"); c1->Update(); } //Calculate no. of events in each pthat for (int i=0;i<2000;i++){ int ipthat; int idata; if (i>= 0 && i< 10) ipthat=0; if (i>= 10 && i< 18) ipthat=1; if (i>= 18 && i< 40) ipthat=2; if (i>= 40 && i< 60) ipthat=3; if (i>= 60 && i< 90) ipthat=4; if (i>= 90 && i< 120) ipthat=5; if (i>=120 && i< 150) ipthat=6; if (i>=150 && i< 200) ipthat=7; if (i>=200 && i< 300) ipthat=8; if (i>=300 && i< 400) ipthat=9; if (i>=400 && i< 500) ipthat=10; if (i>=500 && i< 600) ipthat=11; if (i>=600 && i<1000) ipthat=12; nev[ipthat] =nev[ipthat] + Pthat0->GetBinContent(i+1); }//i //Cross section from pythia pthat bin if (!strcmp(lam.Data(), "1600")){ xsec[0] = 0. ; xsec[1] = 0. ; xsec[2] = 4.961E-02 - 1.317E-03 ; // 18 -40 xsec[3] = 1.317E-03 - 1.757E-04 ; // 40 -60 xsec[4] = 1.757E-04 - 2.027E-05 ; //60 -90 xsec[5] = 2.027E-05 - 3.910E-06 ; //90 - 120 xsec[6] = 3.910E-06 - 1.015E-06 ; //120 - 150 xsec[7] = 1.015E-06 - 1.605E-07 ; //150 - 200 xsec[8] = 1.605E-07 - 1.014E-08 ; // 200 -300 xsec[9] = 1.014E-08 - 1.145E-09 ; // 300 -400 xsec[10] = 1.145E-09 - 1.279E-10 ; // 400 - 500 xsec[11] = 1.279E-10 - 9.534E-12 ; // 500 - 600 xsec[12] = 9.534E-12 ; // > 600 } else if (!strcmp(lam.Data(), "2000")){ xsec[0] = 0. ; xsec[1] = 0. ; xsec[2] = 4.943E-02 - 1.319E-03 ; xsec[3] = 1.319E-03 - 1.756E-04 ; xsec[4] = 1.756E-04 - 2.026E-05 ; xsec[5] = 2.026E-05 - 3.903E-06 ; xsec[6] = 3.903E-06 - 9.977E-07 ; xsec[7] = 9.977E-07 - 1.512E-07 ; xsec[8] = 1.512E-07 - 7.700E-09 ; xsec[9] = 7.700E-09 - 6.664E-10 ; xsec[10] = 6.664E-10 - 6.225E-11 ; xsec[11] = 6.225E-11 - 4.273E-12 ; xsec[12] = 4.273E-12 ; } else if (!strcmp(lam.Data(), "2400")){ xsec[0] = 0. ; xsec[1] = 0. ; xsec[2] = 4.957E-02 - 1.319E-03 ; xsec[3] = 1.319E-03 - 1.756E-04 ; xsec[4] = 1.756E-04 - 2.029E-05 ; xsec[5] = 2.029E-05 - 3.884E-06 ; xsec[6] = 3.884E-06 - 9.929E-07 ; xsec[7] = 9.929E-07 - 1.480E-07 ; xsec[8] = 1.480E-07 - 6.859E-09 ; xsec[9] = 6.859E-09 - 4.960E-10 ; xsec[10] = 4.960E-10 - 3.903E-11 ; xsec[11] = 3.903E-11 - 2.387E-12 ; xsec[12] = 2.387E-12 ; } else if (!strcmp(lam.Data(), "2800")){ xsec[0] = 0. ; xsec[1] = 0. ; xsec[2] = 4.957E-02 - 1.315E-03 ; xsec[3] = 1.315E-03 - 1.760E-04 ; xsec[4] = 1.760E-04 - 2.020E-05 ; xsec[5] = 2.020E-05 - 3.882E-06 ; xsec[6] = 3.882E-06 - 9.926E-07 ; xsec[7] = 9.926E-07 - 1.476E-07 ; xsec[8] = 1.476E-07 - 6.520E-09 ; xsec[9] = 6.520E-09 - 4.233E-10 ; xsec[10] = 4.233E-10 - 2.915E-11 ; xsec[11] = 2.915E-11 - 1.583E-12 ; xsec[12] = 1.583E-12 ; } else if (!strcmp(lam.Data(), "3200")){ xsec[0] = 0. ; xsec[1] = 0. ; xsec[2] = 4.955E-02 - 1.315E-03 ; xsec[3] = 1.315E-03 - 1.760E-04 ; xsec[4] = 1.760E-04 - 2.025E-05 ; xsec[5] = 2.025E-05 - 3.888E-06 ; xsec[6] = 3.888E-06 - 9.923E-07 ; xsec[7] = 9.923E-07 - 1.468E-07 ; xsec[8] = 1.468E-07 - 6.361E-09 ; xsec[9] = 6.361E-09 - 3.902E-10 ; xsec[10] = 3.902E-10 - 2.455E-11 ; xsec[11] = 2.455E-11 - 1.203E-12 ; xsec[12] = 1.203E-12 ; } else if (!strcmp(lam.Data(), "4000")){ xsec[0] = 0. ; xsec[1] = 0. ; xsec[2] = 4.947E-02 - 1.317E-03 ; xsec[3] = 1.317E-03 - 1.760E-04 ; xsec[4] = 1.760E-04 - 2.024E-05 ; xsec[5] = 2.024E-05 - 3.888E-06 ; xsec[6] = 3.888E-06 - 9.939E-07 ; xsec[7] = 9.939E-07 - 1.467E-07 ; xsec[8] = 1.467E-07 - 6.228E-09 ; xsec[9] = 6.228E-09 - 3.641E-10 ; xsec[10] = 3.641E-10 - 2.072E-11 ; xsec[11] = 2.072E-11 - 8.873E-13 ; xsec[12] = 8.873E-13 ; } else cout << "Lambda don't exists " << lam << endl; //Initialize other value of xsec to zero for (int k1 = 13; k1 < nbinX; ++k1){ xsec[k1] =0.; } TH1F* histoMerged = new TH1F("histoMerged", title, nbinY, MjjBin); for (int i =0; i < nbinY; ++i){ float binValue = 0.; float _error = 0.; for (int j =0; j < nbinX; ++j){ float scale; if (nev[j] != 0.) scale = (xsec[j] * 10.)/nev[j]; // in nb, divided by 10E-5 (nevents) else scale = 0.; binValue = binValue + scale*(htemp->GetBinContent(j+1, i+1)); float tmpError = htemp->GetBinError(j+1, i+1); _error = _error + (scale * scale) * ( tmpError * tmpError); }//j float _errorF = TMath::Sqrt(_error); histoMerged->SetBinContent(i+1, binValue); histoMerged->SetBinError(i+1, _errorF); } myps->NewPage(); c1->cd(0); c1->SetLogy(); scale(histoMerged); histoMerged->GetXaxis()->SetTitle("Dijet Mass (in GeV)"); histoMerged->GetXaxis()->SetRangeUser(0., 1400.); histoMerged->GetYaxis()->SetTitle("d#sigma/dm (in nb/GeV)"); histoMerged->SetMarkerStyle(29); histoMerged->SetMarkerColor(2); histoMerged->SetStats(kFALSE); histoMerged->Draw("P"); c1->Update(); // myps->NewPage(); histoMerged->Clone(hName); fileCreate->cd(); histoMerged->Write(hName); for (int i =0; i < nbinX; ++i){ ptHatHisto[i] = 0; delete ptHatHisto[i]; } ++globalK; htemp =0; p =0; // delete [] ptHatHisto; delete htemp; return histoMerged; }//Addition void scale(TH1* &h1) { Int_t nbin = h1->GetNbinsX(); TH1* h2 = (TH1F*) h1->Clone(); // h1->Sumw2(); for (Int_t i =1; i < nbin; ++i){ float binWidth = h2->GetBinWidth(i); float binContent = h2->GetBinContent(i); float value = binContent/binWidth; float error = h2->GetBinError(i); float errVal = error/binWidth; h1->SetBinContent(i, value); h1->SetBinError(i, errVal); } }//scale void JoinHisto(TString hname){ TString had = hname + "h"; TH1F* histoHad = Addition(had); } int Signal(){ //Initialize Input file TString pt[11] = {"018", "040", "060", "090", "120", "150", "200", "300", "400", "500", "600"}; lam = "1600"; TString ext = ".root"; for (int k =0; k < rootFN; ++k){ // TString kFile = "rootFiles/" + lam + "/" + "Lambda" + lam + "Pt" + pt[k] + ".root"; TString kFile = "Lambda" + lam + "Pt" + pt[k] + ".root"; f[k] = new TFile(kFile); }//k c1 = new TCanvas("c1", "c1",3,25,999,799); TString psFile = "Lamda" + lam + ".ps"; myps = new TPostScript(psFile, 112); globalK = 0; Pthat0 = new TH1F(); //Initialize titles titleArray[0] = "Jet.size() >=2 && metSigCut && zvtx < 60."; titleArray[1] = " && 0. <= |rap{jj}| < 0.5"; titleArray[2] = " && SS "; titleArray[3] = " && OS "; titleArray[4] = " && 0.5 <= |rap{jj}| < 1.1"; titleArray[5] = " && 1.1 <= |rap{jj}| < 2.1"; titleArray[6] = "Jet.size() >=2 "; titleArray[7] = titleArray[6] + "&& _metCut "; titleArray[8] = titleArray[7] + "&& _xvtx < 60."; //Outout files TString tOut = "Lam" + lam + ".root"; fileCreate = new TFile(tOut, "RECREATE"); JoinHisto("MjjL5VPthat0"); /* RatioHadOvCal("MjjL5VPthat1"); RatioHadOvCal("MjjL5VPthat2"); RatioHadOvCal("MjjL5VPthat3"); RatioHadOvCal("MjjL5VPthat4"); RatioHadOvCal("MjjL5VPthat5"); RatioHadOvCal("MjjL5VPthat6"); RatioHadOvCal("MjjL5VPthat7"); RatioHadOvCal("MjjL5VPthat8"); RatioHadOvCal("MjjL5VPthat9"); RatioHadOvCal("MjjL5VPthat10"); RatioHadOvCal("MjjL5VPthat11"); */ fileCreate->Close(); for (int k =0; k < rootFN; ++k) f[k]->Close(); myps->Close(); delete c1 ; delete Pthat0; delete [] f; TIter next(gDirectory->GetListOfKeys()); TKey *key; while ((key=(TKey*)next())) { printf("key: %s points to an object of class: %s at %d\n", key->GetName(), key->GetClassName(),key->GetSeekKey()); } gObjectTable->Print(); // gSystem->Exec("gv " + psFile); return 0; }