#include #include gRandom->SetSeed(); TCanvas* c1=NULL; TCanvas* c2=NULL; TCanvas* c3=NULL; TCanvas* c4=NULL; TCanvas* cMT=NULL; TCanvas* cHt=NULL; TCanvas* cJETS=NULL; TCanvas* cCLJ=NULL; TCanvas* cPT=NULL; TCanvas* cPTjets=NULL; TCanvas* cMass=NULL; TCanvas* cTopcharge=NULL; TCanvas* cJetcharge=NULL; TCanvas* cJetchargewl=NULL; TCanvas* cDelta=NULL; TCanvas* cDeltaB=NULL; TCanvas* cDiscriminant=NULL; TCanvas* cHighestBtag=NULL; TCanvas* cNeutrinoCheck=NULL; TCanvas* cMLP1=NULL; TCanvas* cMLP2=NULL; TCanvas* cMLP3=NULL; TH1F *hpmu3x; TH1F *hpmu3y; TH1F *hpmu3z; TH1F *hpWx; TH1F *hpWy; TH1F *hpWz1; TH1F *hpWz2; TH1F *hMT; TH1F *hMTcut; TH1F *hMET; TH1F *hDeltaMET; TH1F *hCLJ; TH1F *hPT; TH1F *hMass; TH1F *hMasscut; TH1F *hTopcharge; TH1F *hJetcharge1; TH1F *hJetcharge2; TH1F *hJetcharge3; TH1F *hJetchargewl1; TH1F *hJetchargewl2; TH1F *hJetchargewl3; TH1F *hHt; TH1F *hLogHt; TH1F *hNj; TH1F *hNjveto; TH1F *hNjBarrel; TH1F *hNjEndcap; TH1F *hNbj; TH1F *hNgj; TH1F *hNuncalib; TH1F *hNcalib; TH1F *hNl; TH1F *hNe; TH1F *hDelta; TH1F *hDeltaB; TH1F *hDiscriminant; TH1F *hEleIso; TH1F *hHighestBtag; TH1F *hBtag1; TH1F *hBtag2; TH1F *hNeutrinoCheck; TH1F *hPtJ1; TH1F *hPtJ2; TH1F *hPtTotal; TH1F *hPhiTB; TH1F *hMLP1; TH1F *hMLP2; TH1F *hMLP3; Double_t PW_x; Double_t PW_y; Double_t PW_z1; Double_t PW_z2; Double_t EW1; Double_t EW2; void top() { ////Arguments for run(): unsigned int verbosity=1, int metopt=1, unsigned int neutrino=0, unsigned int bchoice=0, unsigned int whichlepton=3, Char_t nomesample[40]="TTbar_inclusive", unsigned int numfiles=10, Int_t maxevent=-1, bool plots=true, bool savetree=false run(0,0,4,7,3,"mccard_ewtop_t_l_1_famos_140",2,-1,true,true); // run(0,0,4,7,3,"sm05_st_s_lvb",11,-1,true,true); // compareMET("mtw",5,3,"jm03b_TTbar_leptonic"); // plot("met",-1,4,7,3,"jm03b_TTbar_inclusive",-1,4,7,3,"jm03b_tt_inclusive_Id_883_famos_131"); // plot("met",-1,4,7,3,"jm03b_TTbar_inclusive",-1,4,7,3,"jm03b_tt_inclusive_Id_883_famos_1_3_1"); // plot("deltamet",-1,4,7,3,"jm03b_TTbar_inclusive",-1,4,7,3,"jm03b_tt_inclusive_Id_883_famos_1_2_0"); // plot("deltamet",-1,4,7,3,"jm03b_TTbar_inclusive",-1,4,7,3,"jm03b_tt_inclusive_Id_883_famos_1_3_1"); // plot("deltamet",-1,4,7,3,"jm03b_TTbar_inclusive",-1,4,7,3,"jm03b_tt_inclusive_Id_883_famos_1_3_2_pre3"); // plot("njetveto",0,4,7,3,"jm03b_tt_inclusive_Id_883_famos_1_3_1",0,4,7,3,"jm03b_tt_inclusive_Id_883_famos_1_3_2_pre3"); // plot("pttot",0,4,7,3,"mccard_ewtop_s_l_famos_1_3_1",0,4,7,3,"jm03b_TTbar_inclusive"); // plot("pttot",0,4,7,3,"mccard_ewtop_s_l_famos_1_3_1",0,4,7,3,"wbb_mu_famos_1_3_1"); } void plot(Char_t var[10]="mtop",int metopt1=1,unsigned int neutr1=4, unsigned int bchoice1=7, unsigned int whichlepton1=1, Char_t sample_1[40]="TTbar_leptonic",int metopt2=1,unsigned int neutr2=4, unsigned int bchoice2=7, unsigned int whichlepton2=2,Char_t sample_2[40]="TTbar_inclusive") { gStyle->SetOptStat(0000); TFile *f1; TFile *f2; Char_t temp[80]; Char_t idlep1[2]; Char_t idlep2[2]; if (whichlepton1==1) {sprintf(idlep1,"e");} else if (whichlepton1==2) {sprintf(idlep1,"m");} else {sprintf(idlep1,"em");} if (whichlepton2==1) {sprintf(idlep2,"e");} else if (whichlepton2==2) {sprintf(idlep2,"m");} else {sprintf(idlep2,"em");} sprintf(temp,"histo_%d_%d_%d_%s_%s.root",metopt1,neutr1,bchoice1,idlep1,sample_1); f1 = new TFile(temp); sprintf(temp,"histo_%d_%d_%d_%s_%s.root",metopt2,neutr2,bchoice2,idlep2,sample_2); f2 = new TFile(temp); f1->ls(); f2->ls(); TH1F *histo1 = (TH1F*)f1->Get(var); TH1F *histo2 = (TH1F*)f2->Get(var); ///// normalizzazione: Double_t scale; if (histo1->Integral()!=0) scale = 1/histo1->Integral(); histo1->Scale(scale); if (histo2->Integral()!=0) scale = 1/histo2->Integral(); histo2->Scale(scale); //// plot degli istogrammi: // c1=new TCanvas("c1",var,0,0,600,400); TCanvas c1=new TCanvas("c1",var,2); c1->cd(); histo1->SetLineWidth(2); histo2->SetLineWidth(2); histo1->SetTitle(""); histo2->SetTitle(""); Char_t titolo[80]; if (var=="mtw") {sprintf(titolo,"W transv. mass (GeV)");} if (var=="met") {sprintf(titolo,"missing ET (GeV)");} if (var=="deltamet") {sprintf(titolo,"missing ET resolution");} if (var=="mtop") {sprintf(titolo,"top mass (GeV)");} if (var=="njet") {sprintf(titolo,"N. of jets");} if (var=="ngjet") {sprintf(titolo,"N. of generated jets");} if (var=="nuncalib") {sprintf(titolo,"N. of uncalib. jets");} if (var=="ncalib") {sprintf(titolo,"N. of calib. jets");} if (var=="njetbarrel") {sprintf(titolo,"N. of jets in the barrel");} if (var=="njetendcap") {sprintf(titolo,"N. of jets in the endcap");} if (var=="nbjet") {sprintf(titolo,"N. of b-jets");} if (var=="njetveto") {sprintf(titolo,"N. of vetoing jets");} if (var=="nlep") {sprintf(titolo,"N. of leptons");} if (var=="nele") {sprintf(titolo,"N. of electrons");} if (var=="ht") {sprintf(titolo,"Ht (GeV)");} if (var=="loght") {sprintf(titolo,"log(Ht)");} if (var=="discriminant") {sprintf(titolo,"Discriminant");} if (var=="eiso") {sprintf(titolo,"Electron isolation");} if (var=="coslj") {sprintf(titolo,"cos theta_lj");} if (var=="delta") {sprintf(titolo,"Delta top true / top reco");} if (var=="deltab") {sprintf(titolo,"Delta b true / b reco");} if (var=="nucheck") {sprintf(titolo,"DeltaZ W true / W reco");} if (var=="highestbtag") {sprintf(titolo,"Highest btag jet");} if (var=="ptj1") {sprintf(titolo,"Pt jet 1, GeV");} if (var=="ptj2") {sprintf(titolo,"Pt jet 2, GeV");} if (var=="btag1") {sprintf(titolo,"b-tag jet 1, GeV");} if (var=="btag2") {sprintf(titolo,"b-tag jet 2, GeV");} if (var=="pttot") {sprintf(titolo,"Pt of top+bottom system");} if (var=="mlp1") {sprintf(titolo,"NN 1");} if (var=="mlp2") {sprintf(titolo,"NN 2");} if (var=="mlp3") {sprintf(titolo,"NN 3");} histo1->SetXTitle(titolo); histo1->DrawCopy(); histo2->SetLineColor(kRed); histo2->DrawCopy("same"); ///// legenda: Char_t met1[80]; if (metopt1==-3) {sprintf(met1,"MET from MC, smeared");} if (metopt1==-2) {sprintf(met1,"MET from MC");} if (metopt1==-1) {sprintf(met1,"MET from calo towers");} if (metopt1==0) {sprintf(met1,"MET from jets");} if (metopt1==1) {sprintf(met1,"MET type 1");} if (metopt1==3) {sprintf(met1,"MET from true Nu");} Char_t met2[80]; if (metopt2==-3) {sprintf(met2,"MET from MC, smeared");} if (metopt2==-2) {sprintf(met2,"MET from MC");} if (metopt2==-1) {sprintf(met2,"MET from calo towers");} if (metopt2==0) {sprintf(met2,"MET from jets");} if (metopt2==1) {sprintf(met2,"MET type 1");} if (metopt2==3) {sprintf(met1,"MET from true Nu");} Char_t strNeu1[80]; if (neutr1<3) {sprintf(strNeu1,"neutrino sol.: random");} if (neutr1==3) {sprintf(strNeu1,"neutrino sol.: best Mtop");} if (neutr1==4) {sprintf(strNeu1,"neutrino sol.: min Pz(W)");} if (neutr1==5) {sprintf(strNeu1,"neutrino sol.: min Pz(top)");} Char_t strNeu2[80]; if (neutr2<3) {sprintf(strNeu2,"neutrino sol.: random");} if (neutr2==3) {sprintf(strNeu2,"neutrino sol.: best Mtop");} if (neutr2==4) {sprintf(strNeu2,"neutrino sol.: min Pz(W)");} if (neutr2==5) {sprintf(strNeu2,"neutrino sol.: min Pz(top)");} Double_t ymax=max(histo1->GetMaximum(),histo2->GetMaximum()); Double_t xmax; Double_t xmin; if (var=="mtw") {xmin=5; xmax=45;} if (var=="met") {xmin=80; xmax=120;} if (var=="deltamet") {xmin=1.7; xmax=4.9;} if (var=="mtop") {xmin=30; xmax=115;} if (var=="njet") {xmin=7; xmax=9.8;} if (var=="ngjet") {xmin=7; xmax=9.8;} if (var=="nuncalib") {xmin=7; xmax=9.8;} if (var=="ncalib") {xmin=7; xmax=9.8;} if (var=="njetbarrel") {xmin=7; xmax=9.8;} if (var=="njetendcap") {xmin=7; xmax=9.8;} if (var=="nbjet") {xmin=7; xmax=9.8;} if (var=="njetveto") {xmin=7; xmax=9.8;} if (var=="nlep") {xmin=3; xmax=5;} if (var=="nele") {xmin=3; xmax=4;} if (var=="ht") {xmin=430; xmax=600;} if (var=="loght") {xmin=1; xmax=3;} if (var=="discriminant") {xmin=5000; xmax=20000;} if (var=="eiso") {xmin=0.6; xmax=1.85;} if (var=="coslj") {xmin=-1; xmax=-0.4;} if (var=="delta") {xmin=4; xmax=6;} if (var=="deltab") {xmin=4; xmax=6;} if (var=="nucheck") {xmin=75; xmax=100;} if (var=="highestbtag") {xmin=1; xmax=2; ymax=ymax/2.;} if (var=="ptj1") {xmin=170; xmax=250;} if (var=="ptj2") {xmin=170; xmax=250;} if (var=="btag1") {xmin=2; xmax=4.9;} if (var=="btag2") {xmin=2; xmax=4.9;} if (var=="pttot") {xmin=160; xmax=240;} if (var=="mlp1") {xmin=0.65; xmax=1.;} if (var=="mlp2") {xmin=0.65; xmax=1.;} if (var=="mlp3") {xmin=0.65; xmax=1.;} /* TLegend *legend = new TLegend(xmin,0.65*ymax,xmax,0.95*ymax,"",""); sprintf(temp,"%s, lept.: %s, neutr.sol.:%d, %s",sample_1,idlep1,neutr1,met1); legend->AddEntry(histo1,temp,"l"); sprintf(temp,"%s, lept.: %s, neutr.sol.:%d, %s",sample_2,idlep2,neutr2,met2); legend->AddEntry(histo2,temp,"l"); legend->Draw(); */ float mean1=histo1->GetMean(); float rms1=histo1->GetRMS(); float mean2=histo2->GetMean(); float rms2=histo2->GetRMS(); char strMean[20]; char strRMS[20]; TPavesText *tpt1 = new TPavesText(xmin,0.75*ymax,xmax,0.99*ymax,1,"br"); tpt1->SetTextSize(0.02); tpt1->SetTextColor(kBlack); tpt1->AddText(sample_1); sprintf(temp,"lept.: %s, %s",idlep1,met1); tpt1->AddText(temp); tpt1->AddText(strNeu1); sprintf(strMean,"mean= %f",mean1); // sprintf(strRMS,"rms= %f",rms1); tpt1->AddText(strMean); // tpt1->AddText(strRMS); tpt1->Draw(); TPavesText *tpt2 = new TPavesText(xmin,0.50*ymax,xmax,0.75*ymax,1,"br"); tpt2->SetTextSize(0.02); tpt2->SetTextColor(kRed); tpt2->AddText(sample_2); sprintf(temp,"lept.: %s, %s",idlep2,met2); tpt2->AddText(temp); tpt2->AddText(strNeu2); sprintf(strMean,"mean= %f",mean2); // sprintf(strRMS,"rms= %f",rms2); tpt2->AddText(strMean); // tpt2->AddText(strRMS); tpt2->Draw(); // c1->Update(); ///// salva su file PS: sprintf(temp,"%s_%d_%d_%d_%s_%s_%d_%d_%d_%s_%s.ps",var,metopt1,neutr1,bchoice1,idlep1,sample_1,metopt2,neutr2,bchoice2,idlep2,sample_2); c1->Print(temp); } // to-do-list: /* - void plot_allsamples(input: var, opzioni met/jet/neutrino) confronta tutte le categorie di eventi in un solo grafico - void plot_alljetmetopt(input: var, nomesample) confronta diverse opzioni jet/met - void plot_allnuopt(input: var, nomesample) confronta diverse soluzioni di neutrino */ void compareMET (Char_t var[10]="mtw",unsigned int neutr=4, unsigned int bchoice=7, unsigned int whichlepton=3, Char_t sample[40]="WWjets_inclusive") { gStyle->SetOptStat(0000); TFile *f1; TFile *f2; TFile *f3; TFile *f4; Char_t temp[80]; Char_t idlep[2]; if (whichlepton==1) {sprintf(idlep,"e");} else if (whichlepton==2) {sprintf(idlep,"m");} else {sprintf(idlep,"em");} // sprintf(temp,"histo_2_%d_%s_%s.root",neutr,idlep,sample);// calibrated met sprintf(temp,"histo_3_%d_%d_%s_%s.root",neutr,bchoice,idlep,sample);// W->nu f1 = new TFile(temp); sprintf(temp,"histo_-1_%d_%d_%s_%s.root",neutr,bchoice,idlep,sample);// met from towers f2 = new TFile(temp); sprintf(temp,"histo_0_%d_%d_%s_%s.root",neutr,bchoice,idlep,sample);// met from IC jets f3 = new TFile(temp); sprintf(temp,"histo_-2_%d_%d_%s_%s.root",neutr,bchoice,idlep,sample);//met from mc particles f4 = new TFile(temp); sprintf(temp,"histo_-3_%d_%d_%s_%s.root",neutr,bchoice,idlep,sample);//met from mc particles, smeared in energy and phi f5 = new TFile(temp); f1->ls(); f2->ls(); f3->ls(); f4->ls(); f5->ls(); TH1F *histo1 = (TH1F*)f1->Get(var); TH1F *histo2 = (TH1F*)f2->Get(var); TH1F *histo3 = (TH1F*)f3->Get(var); TH1F *histo4 = (TH1F*)f4->Get(var); TH1F *histo5 = (TH1F*)f5->Get(var); ///// normalizzazione: Double_t scale=1.; if (histo1->Integral()!=0) scale = 1/histo1->Integral(); histo1->Scale(scale); if (histo2->Integral()!=0) scale = 1/histo2->Integral(); histo2->Scale(scale); if (histo3->Integral()!=0) scale = 1/histo3->Integral(); histo3->Scale(scale); if (histo4->Integral()!=0) scale = 1/histo4->Integral(); histo4->Scale(scale); if (histo5->Integral()!=0) scale = 1/histo5->Integral(); histo5->Scale(scale); //// plot degli istogrammi: // c1=new TCanvas("c1",var,0,0,600,400); c1=new TCanvas("c1",var,2); c1->cd(); histo1->SetLineWidth(2); histo2->SetLineWidth(2); histo3->SetLineWidth(2); histo4->SetLineWidth(2); histo5->SetLineWidth(2); histo4->SetLineColor(kBlue); histo4->DrawCopy(); histo2->SetLineColor(kRed); histo2->DrawCopy("same"); histo3->SetLineColor(41); histo3->DrawCopy("same"); histo5->SetLineColor(kGreen); histo5->DrawCopy("same"); histo1->SetLineColor(kBlack); histo1->DrawCopy("same"); ///// legenda: Double_t ymax=max(histo1->GetMaximum(),histo2->GetMaximum(),histo3->GetMaximum(),histo4->GetMaximum(),histo5->GetMaximum()); Double_t xmax; Double_t xmin; if (var=="mtw") {xmin=5; xmax=50;} if (var=="mtop") {xmin=30; xmax=115;} if (var=="ht") {xmin=530; xmax=690;} if (var=="met") {xmin=80; xmax=120;} sprintf(temp,"%s, lept.: %s",sample,idlep); TLegend *legend = new TLegend(xmin,0.65*ymax,xmax,0.95*ymax,temp,""); legend->AddEntry(histo1,"calibrated MET","l"); legend->AddEntry(histo2,"raw MET","l"); legend->AddEntry(histo3,"MET from jets","l"); legend->AddEntry(histo4,"MET from MC truth","l"); legend->AddEntry(histo5,"MET from MC truth, smeared","l"); legend->Draw(); c1->Update(); ///// salva su file PS: sprintf(temp,"%s_metcomparison_%d_%s_%s.ps",var,neutr,idlep,sample); c1->Print(temp); } void run(unsigned int verbosity=1, int metopt=1, unsigned int neutrino=3, unsigned int bchoice=7, unsigned int whichlepton=2, Char_t nomesample[40]="TTbar_inclusive", unsigned int numfiles=10, Int_t maxevent=-1, bool plots=true, bool savetree=true) { gRandom->SetSeed(); W* myW; const Int_t maxtop = 2; Top* myT[maxtop]; //Reset ROOT gROOT->Reset(); // Load shared libraries gSystem->Load("$ROOTSYS/lib/libPhysics.so"); gSystem->Load("libExRootAnalysisReader"); // Create chain of root trees bool giovanni=false; // temporary variable TChain chain("Analysis"); for (Int_t i=1;i<=numfiles;i++){ // TString name("~/scratch0/rootfiles/"); Int_t castor=1; // temporary TString name; if (castor==0) name = new TString("~/scratch0/rootfiles/"); if (castor==1) name = new TString("rfio:/castor/cern.ch/user/g/giamman/"); name+=nomesample; if (!giovanni) name+="/test_"; if (giovanni) name+="/famos_mccard_ewtop_wt_ll_"; name+=i; if (giovanni) name+=".0"; name+=".root"; // TFile tempfile(name); TFile* tempfile; if (castor==0) { tempfile = new TFile(name); } else { tempfile = new TRFIOFile(name); } if (tempfile.IsOpen()) { cout << name << " aperto " << endl; chain.Add(name); } else { cout << "################# " << name << " non aperto " << endl; } } // Create object of class ExRootTreeReader ExRootTreeReader *myTree = new ExRootTreeReader(&chain); // New variables & cuts: Int_t Njets; Double_t pj_x[50]; Double_t pj_y[50]; Double_t pj_z[50]; Double_t pj_e[50]; Int_t jetcode[50]; Double_t jetcharge[50]; Double_t jetchargewl[50]; Double_t btag[50]; // bool istagged[50]; Double_t pe_x[15]; Double_t pe_y[15]; Double_t pe_z[15]; Double_t pe_e[15]; Double_t echarge[15]; Double_t eiso[15]; Int_t electrontype[15]; Double_t pmu_x[6]; Double_t pmu_y[6]; Double_t pmu_z[6]; Double_t pmu_e[6]; Double_t mucharge[6]; Int_t muontype[15]; Double_t pl_x[16]; Double_t pl_y[16]; Double_t pl_z[16]; Double_t pl_e[16]; Double_t lcharge[16]; Double_t isA[16]; Double_t PtMu; Double_t PtMuCut=10; Double_t PtE; Double_t PtECut=10; Double_t PtL; // Double_t PtLCut=20; Double_t PtJCut=30; Double_t missingEx, missingEy, missingET; Double_t missingET_MC, missingEx_MC, deltaMET, deltaMET_MET; Double_t METCut=20; Double_t BTagCut=0.; Double_t antiBTagCut=0.; Double_t EtaECut=2.5; // Anne-Sylvie usa 2.0 Double_t eopCut=0.8; // Anne-Sylvie Double_t invEPCut=0.8; // Anne-Sylvie Double_t HcalEcalCut=0.05; // Anne-Sylvie Double_t DRCut=0.15; // Anne-Sylvie Double_t DEtaTrkClusterCut=0.005; // Anne-Sylvie Double_t EleIsoCut=0.05; // Anne-Sylvie Double_t isoR=0.35; // Anne-Sylvie Double_t isoRveto=0.015; // Domenico Double_t EtaMuCut=2.1; Double_t cos_lj; Int_t whichtop; Int_t nontop; Double_t discriminant; Double_t PtJ1; Double_t PtJ2; Double_t PtTot; Double_t EtaJ1; Double_t EtaJ2; Double_t PhiJ1; Double_t PhiJ2; Double_t DeltaJ1J2; Double_t PtJ3=0.; // will be used only if njveto=1 is allowed Double_t btag1; Double_t btag2; Double_t mtw; Double_t mtop; Double_t ht; Double_t mtb; Double_t pttop; Double_t topcharge; Int_t nbins=20; Float_t low; Float_t high; low=-100; high=100; hpmux = new TH1F("mux","Px",nbins,low,high); hpmuy = new TH1F("muy","Py",nbins,low,high); hpmuz = new TH1F("muz","Pz",nbins,low,high); hpWx = new TH1F("Wx","Px",nbins,low,high); hpWy = new TH1F("Wy","Py",nbins,low,high); hpWz1 = new TH1F("Wz1","Pz",nbins,low,high); hpWz2 = new TH1F("Wz2","Pz",nbins,low,high); low=0; high=120; hMT = new TH1F("mtw","mtw",nbins,low,high); hMTcut = new TH1F("mtwcut","mtwcut",nbins,low,high); low=0; high=120; hMET = new TH1F("met","met",nbins,low,high); low=-5; high=5; hDeltaMET = new TH1F("deltamet","deltamet",nbins,low,high); low=-1; high=1; hCLJ = new TH1F("coslj","coslj",nbins,low,high); low=0; high=250; hPT = new TH1F("pttop","pttop",nbins,low,high); low=70; high=300; hMass = new TH1F("mtop","mtop",nbins,low,high); hMasscut = new TH1F("mtopcut","mtopcut",nbins,low,high); low=-1; high=1; hTopcharge = new TH1F("topch","topch",nbins,low,high); low=-1; high=1; hJetcharge1 = new TH1F("jch1","jch1",nbins,low,high); hJetcharge2 = new TH1F("jch2","jch2",nbins,low,high); hJetcharge3 = new TH1F("jch3","jch3",nbins,low,high); low=-1; high=1; hJetchargewl1 = new TH1F("jchw1","jchw1",nbins,low,high); hJetchargewl2 = new TH1F("jchw2","jchw2",nbins,low,high); hJetchargewl3 = new TH1F("jchw3","jchw3",nbins,low,high); low=-10; high=600; hHt = new TH1F("ht","ht",nbins,low,high); low=3; high=8; hLogHt = new TH1F("loght","loght",nbins,low,high); low=0; high=10; nbins=high; hNj = new TH1F("njet","njet",nbins,low,high); hNjveto = new TH1F("njetveto","njetveto",nbins,low,high); hNgj = new TH1F("ngjet","ngjet",nbins,low,high); hNjBarrel = new TH1F("njetbarrel","njetbarrel",nbins,low,high); hNjEndcap = new TH1F("njetendcap","njetendcap",nbins,low,high); hNuncalib = new TH1F("nuncalib","nuncalib",nbins,low,high); hNcalib = new TH1F("ncalib","ncalib",nbins,low,high); low=0; high=7; nbins=high; hNbj = new TH1F("nbjet","nbjet",nbins,low,high); low=0; high=5; nbins=high; hNl = new TH1F("nlep","nlep",nbins,low,high); low=0; high=5; nbins=high; hNe = new TH1F("nele","nele",nbins,low,high); low=0; high=6.3; nbins=50; hDelta = new TH1F("delta","delta",nbins,low,high); low=0; high=6.3; nbins=50; hDeltaB = new TH1F("deltab","deltab",nbins,low,high); low=-10000; high=20000; nbins=50; hDiscriminant = new TH1F("discriminant","discriminant",nbins,low,high); low=-0.1; high=2; nbins=50; hEleIso = new TH1F("eiso","eiso",nbins,low,high); low=1; high=3; nbins=2; hHighestBtag = new TH1F("highestbtag","highestbtag",nbins,low,high); low=-3; high=5; nbins=50; hBtag1 = new TH1F("btag1","btag1",nbins,low,high); hBtag2 = new TH1F("btag2","btag2",nbins,low,high); low=0; high=100; nbins=50; hNeutrinoCheck = new TH1F("nucheck","nucheck",nbins,low,high); low=0; high=250; hPtJ1 = new TH1F("ptj1","ptj1",nbins,low,high); low=0; high=250; hPtJ2 = new TH1F("ptj2","ptj2",nbins,low,high); low=0; high=250; hPtTotal = new TH1F("pttot","pttot",nbins,low,high); low=0; high=6.3; hPhiTB = new TH1F("phitb","phitb",nbins,low,high); low=-0.1; high=1.1; hMLP1 = new TH1F("mlp1","mlp1",nbins,low,high); hMLP2 = new TH1F("mlp2","mlp2",nbins,low,high); hMLP3 = new TH1F("mlp3","mlp3",nbins,low,high); //Declaration of leaves types Long64_t nEntries = myTree->GetEntries(); TClonesArray *tracks = myTree->UseBranch("TrkCmb"); TClonesArray *btagJets = myTree->UseBranch("BTagJet"); TClonesArray *genJets = myTree->UseBranch("GenJetIC5A"); // TClonesArray *ktJets = myTree->UseBranch("JetKT4"); TClonesArray *itConeJets = myTree->UseBranch("JetIC5A"); TClonesArray *itCalibConeJets = myTree->UseBranch("JetIC5C"); TClonesArray *muons = myTree->UseBranch("MuonGLB"); TClonesArray *electrons = myTree->UseBranch("Elec"); TClonesArray *electronTracks = myTree->UseBranch("TrkEP"); if (metopt<-1) { // MET from MC truth (-2: as is; -3: smeared) TClonesArray *met = myTree->UseBranch("METMC"); } else if (metopt==-1) { // raw (from calo towers) TClonesArray *met = myTree->UseBranch("METCT"); } else if (metopt==0) { // uncalibrated; use MET from jets TClonesArray *met = myTree->UseBranch("METIC"); } else if (metopt==1) { // calibrated; use METCT for towers TClonesArray *met = myTree->UseBranch("METCT"); // } else if (metopt==2) { // calibrated; use RawMET for towers // TClonesArray *met = myTree->UseBranch("RawMET"); } else if (metopt==2) { // calibrated; use METIC TClonesArray *met = myTree->UseBranch("METIC"); } else if (metopt==3) { // true neutrino from W cout << "*** No RecCollection is used for missingET ***" << endl; } TClonesArray *metmc = myTree->UseBranch("METMC"); // for reference, I take it anyway TClonesArray *partons = myTree->UseBranch("Gen"); // TClonesArray *simparticles = myTree->UseBranch("SimPart"); // Analysis loop Int_t nmax; if (maxevent<0 || maxevent>nEntries) { nmax=nEntries; } else { nmax=maxevent; } // counters: Int_t nSelected[3]={0,0,0}; Int_t nPreselected=0; Int_t nPreselectedElectrons=0; Int_t nPreselectedMuons=0; Int_t nCut1=0; Int_t nCut2=0; Int_t nCut3=0; Int_t nFinal=0; Int_t nosol=0; Int_t goodsol=0; Int_t plusTop=0; Int_t minusTop=0; Int_t plusLepton=0; Int_t minusLepton=0; Int_t plusLeptonM1=0; Int_t minusLeptonM1=0; Int_t plusLeptonM2=0; Int_t minusLeptonM2=0; Int_t plusLeptonM3=0; Int_t minusLeptonM3=0; Int_t plusLeptonPt1=0; Int_t minusLeptonPt1=0; Int_t plusLeptonPt2=0; Int_t minusLeptonPt2=0; Int_t plusLeptonPt3=0; Int_t minusLeptonPt3=0; Int_t plusLeptonPt4=0; Int_t minusLeptonPt4=0; Int_t countnubest=0; Int_t countnuworst=0; Int_t countbbest=0; Int_t countbworst=0; Int_t countbout=0; // define the roottree with the NN variables: Int_t use4train; TTree *nntree; if (savetree) { nntree = new TTree("T","nn inputs"); nntree->Branch("btag1",&btag1,"btag1/D"); nntree->Branch("btag2",&btag2,"btag2/D"); nntree->Branch("ptj1",&PtJ1,"ptj1/D"); nntree->Branch("ptj2",&PtJ2,"ptj2/D"); nntree->Branch("delta12",&DeltaJ1J2,"delta12/D"); nntree->Branch("eta2",&EtaJ2,"eta2/D"); nntree->Branch("mtop",&mtop,"mtop/D"); nntree->Branch("pttop",&pttop,"pttop/D"); nntree->Branch("pttot",&PtTot,"pttot/D"); nntree->Branch("ht",&ht,"ht/D"); nntree->Branch("mtb",&mtb,"mtb/D"); nntree->Branch("nmax",&nmax,"nmax/I");//useful for efficiency rescaling in nn.C nntree->Branch("use4train",&use4train,"use4train/I");// 1: use this event for NN training or cut optimization or GA or so on; 0: use for test/validation/efficiencies /* nntree->Branch("",&,"/D"); nntree->Branch("",&,"/D"); nntree->Branch("",&,"/D"); nntree->Branch("",&,"/D"); nntree->Branch("",&,"/D"); nntree->Branch("",&,"/D"); nntree->Branch("",&,"/D"); nntree->Branch("",&,"/D"); */ } // other NN stuff: /* Double_t inputs[8]; if (!gROOT->GetClass("TMultiLayerPerceptron")) { gSystem->Load("libMLP"); cout << endl << "ROOT neural nets library has just been loaded." << endl << endl; } TFile *fnn1 = new TFile("sch_vs_ttincl.root"); TTree *temptree1 = (TTree*)gDirectory->Get("T"); // metto "T" perche' e' il nome del tree contenuto nel file TFile *fnn2 = new TFile("sch_vs_wbb.root"); TTree *temptree2 = (TTree*)gDirectory->Get("T"); // metto "T" perche' e' il nome del tree contenuto nel file TFile *fnn3 = new TFile("sch_vs_ttincl.root"); TTree *temptree3 = (TTree*)gDirectory->Get("T"); // metto "T" perche' e' il nome del tree contenuto nel file Int_t neuro1=4; Int_t neuro2=2; TString struttura("ptj1,ptj2,btag1,btag2,mtw,mtop,ht,pttop"); // inputs struttura+=":"; struttura+=neuro1; if (neuro2>0) { // due strati nascosti struttura+=":"; struttura+=neuro2; } struttura+=":"; struttura+="@truth"; //target output; la @ significa che la rete lo considera normalizzato Double_t truth; TMultiLayerPerceptron *mlp1 = new TMultiLayerPerceptron(struttura,temptree1); mlp1->LoadWeights("sch_vs_ttincl_Batch_4_2_4500.txt"); TMultiLayerPerceptron *mlp2 = new TMultiLayerPerceptron(struttura,temptree2); mlp2->LoadWeights("sch_vs_wbb_Batch_4_0_750.txt"); TMultiLayerPerceptron *mlp3 = new TMultiLayerPerceptron(struttura,temptree3); mlp3->LoadWeights("sch_vs_ttincl_Batch_4_2_4500.txt"); */ // EVENT LOOP: for (Int_t i=0; iReadEntry(i); int nBtagJet=btagJets->GetEntries(); // int nKtJet=ktJets->GetEntries(); // int nItConeJet=itConeJets->GetEntries(); int nTracks=tracks->GetEntries(); int nMuon=muons->GetEntries(); int nElec=electrons->GetEntries(); int nMet=met->GetEntries(); int nPartons=partons->GetEntries(); // int nSimPart=simparticles->GetEntries(); Njets=nBtagJet; //// GENERATOR LEVEL: True tops and W's and b's: Int_t nTrueTops=0; Int_t nTrueWs=0; Int_t nTrueBs=0; Int_t whichparton=-1; // Int_t array_whichW[2]={-2,-2}; Int_t whichW=-1; Int_t whichB=-1; // Int_t array_motherofW[2]={-2,-2}; Int_t motherofW=0; Int_t motherofB=0; for (Int_t q=0; qAt(q); if (abs(myParton->PID)==6) { // top whichparton=q; nTrueTops++; } if (abs(myParton->PID)==24 && myParton->Status==3) { // W // array_whichW[nTrueWs]=q; // array_motherofW[nTrueWs]=myParton->M1; whichW=q; motherofW=myParton->M1; nTrueWs++; // in s-channel you have 2 Ws, but only the second from t->Wb } } Int_t nTrueNus=0; Int_t nTrueLs=0; Int_t whichNu=-1; Int_t whichL=-1; Int_t motherofNu=0; Int_t motherofL=0; Int_t pidNu=0; Int_t pidL=0; bool useNu=false; for (Int_t q=0; qAt(q); if ((abs(myParton->PID)==12 || abs(myParton->PID)==14 || abs(myParton->PID)==16) && myParton->Status==3) { // neutrino pidNu=myParton->PID; whichNu=q; motherofNu=myParton->M1; nTrueNus++; } if ((abs(myParton->PID)==11 || abs(myParton->PID)==13 || abs(myParton->PID)==15) && myParton->Status==3) { // lepton pidL=myParton->PID; whichL=q; motherofL=myParton->M1; nTrueLs++; } } if (nTrueNus==1 && nTrueLs==1 && motherofL==motherofNu) { // nu-lepton pairing if (pidNu==-1*(pidL+1)) { // further check useNu=true; } } for (Int_t q=0; qAt(q); if (abs(myParton->PID)==5 && myParton->Status==3) { // b from t if (verbosity==1) cout << " trovato b" << endl; motherofB=myParton->M1; if (motherofB>0) { if (motherofB==motherofW) { whichB=q; nTrueBs++; if (verbosity==1) cout << " t->Wb" << endl; } else { if (verbosity==1) cout << " nessun W e' fratello di questo b" << endl; } } } } if (nTrueTops==1) { TRootGenParticle *myTrueTop = (TRootGenParticle*) partons->At(whichparton); TrueTop = new TLorentzVector(myTrueTop->Px,myTrueTop->Py,myTrueTop->Pz,myTrueTop->E); if (nTrueBs==1) { TRootGenParticle *myTrueB = (TRootGenParticle*) partons->At(whichB); TrueB = new TLorentzVector(myTrueB->Px,myTrueB->Py,myTrueB->Pz,myTrueB->E); trueBeta=myTrueB->Eta; trueBphi=myTrueB->Phi; // also the W has to be put here, since t->Wq decays are not included TRootGenParticle *myTrueW = (TRootGenParticle*) partons->At(whichW); TrueW = new TLorentzVector(myTrueW->Px,myTrueW->Py,myTrueW->Pz,myTrueW->E); trueWz=myTrueW->Pz; } } if (nTrueNus==1 && useNu) { if (verbosity == 1) cout << "uso il nu" << endl; TRootGenParticle *myTrueNu = (TRootGenParticle*) partons->At(whichNu); TrueNeutrino = new TLorentzVector(myTrueNu->Px,myTrueNu->Py,myTrueNu->Pz,myTrueNu->E); } else { TrueNeutrino = new TLorentzVector(0,0,1,1); } //// Electron preselection: Int_t nGoodElec=0; Double_t tkCharge; for (Int_t j=0; jAt(j); Double_t eop=myElectron->EoverP; if (eopEHadOverE; if (he>HcalEcalCut) good=false; Double_t deltaeta=myElectron->DEtaTrkCluster; if (deltaeta>DEtaTrkClusterCut) good=false; if (deltaeta>0.15) good=false; Double_t pixel=myElectron->TrkPixLines; // cout << "pixel lines: " << pixel << endl; // arrivano fino a 7-8! Int_t tkiso=myElectron->TrkIsolation; // cout << "---tkiso: " << tkiso << endl; if (!tkiso) good=false; // ma e' sempre 1!!!! int tkstatus=myElectron->TrkStatus; Double_t tkP=myElectron->TrkP; Double_t tkPt=myElectron->TrkPT; // if (tkPtCalE; if (fabs((1./calE)-(1./tkP))>invEPCut) good=false; Double_t calEt=myElectron->CalET; if (calEtTrkPhi; Double_t tkEta=myElectron->TrkEta; // if (fabs(tkEta)>EtaECut) good=false; Double_t calEta=myElectron->CalEta; Double_t calPhi=myElectron->CalPhi; if (fabs(calEta)>EtaECut) good=false; if (deltaR(calEta,calPhi,tkEta,tkPhi)>DRCut) good=false; //// remove gamma conversions: for (Int_t j2=0; j2At(j2); Double_t tkPhi2=myElectron2->TrkPhi; Double_t tkEta2=myElectron2->TrkEta; Double_t deltaee=deltaR(tkEta,tkPhi,tkEta2,tkPhi2); if (deltaee<0.01) { good=false; //// matching with MC: /* for (Int_t k=0;kAt(k); Int_t pid=mySim->PID; if (abs(pid)==11) { Double_t etasim=mySim->Eta; Double_t phisim=mySim->Phi; Double_t deltaemc=deltaR(tkEta2,tkPhi2,etasim,phisim); cout << " deltaee = " << deltaee << ", deltaemc = " << deltaemc << endl; } } */ } } } //// associated EP track: track_elec = (TRootTrack*) myElectron->Trk.GetObject(); if (!track_elec) good=false; tkCharge=track_elec->Charge; //// isolation: Double_t dca_closest, z_closest; dca_closest=sqrt((track_elec->dcaX)**2+(track_elec->dcaY)**2); z_closest=track_elec->dcaZ; Double_t dca_tk, z_tk; Double_t minDelta=999.; Double_t iso=0.; for (Int_t l=0; lAt(l); // if (&myTrack == &track_elec) cout "************" << endl << " trovato! " << "************" << endl; Double_t etatk=myTrack->Eta; Double_t phitk=myTrack->Phi; // if (deltaR(calEta,calPhi,etatk,phitk)isoRveto) { // Domenico dca_tk=sqrt((myTrack->dcaX)**2+(myTrack->dcaY)**2); z_tk=myTrack->dcaZ; if (myTrack->PT > 0.9 && fabs(z_closest-z_tk)<0.4 && fabs(dca_closest-dca_tk)<0.1 && myTrack->Hits > 4) iso=iso+myTrack->PT; // Anne-Sylvie } } // iso=(iso-myElectron->TrkPT)/myElectron->CalET; // Anne-Sylvie iso=iso/myElectron->CalET; // Domenico (the leading track has been removed with the isoRveto cone) if (good) hEleIso->Fill(iso); if (iso>EleIsoCut) good=false; //// Properties of the preselected electron: if (good) { pe_x[nGoodElec]=eop*(myElectron->TrkPx); pe_y[nGoodElec]=eop*(myElectron->TrkPy); pe_z[nGoodElec]=eop*(myElectron->TrkPz); pe_e[nGoodElec]=myElectron->CalE; echarge[nGoodElec]=tkCharge; eiso[nGoodElec]=iso; nGoodElec++; } } hNe->Fill(nGoodElec); //// Muon preselection: Int_t nGoodMuon=0; for (Int_t j=0; jAt(j); Double_t pt=myMuon->PT; Double_t eta=myMuon->Eta; if (ptEtaMuCut) good=false; if (good) { pmu_x[nGoodMuon]=myMuon->Px; pmu_y[nGoodMuon]=myMuon->Py; pmu_z[nGoodMuon]=myMuon->Pz; pmu_e[nGoodMuon]=myMuon->P; mucharge[nGoodMuon]=myMuon->Charge; nGoodMuon++; } } //// Generic leptons: Int_t nLept; if (whichlepton==1) { // only electrons nLept=nGoodElec; for (Int_t j=0; jGetEntries(); int nConeJets=itConeJets->GetEntries(); int nCalibConeJets=itCalibConeJets->GetEntries(); int ngjet=0; TRootJet *myJ; for (Int_t j=0; jAt(j); Double_t PtJ=myJ->PT; if (PtJ>PtJCut) { ngjet++; } } hNgj->Fill(ngjet); int nuncalib=0; TRootJet *myJ; for (Int_t j=0; jAt(j); Double_t PtJ=myJ->PT; if (PtJ>PtJCut) { nuncalib++; } } hNuncalib->Fill(nuncalib); int ncalib=0; for (Int_t j=0; jAt(j); Double_t PtJ=myJ->PT; if (PtJ>PtJCut) { ncalib++; } } hNcalib->Fill(ncalib); // MET calibration for (Int_t j=0; jAt(j); if (metopt<1 && metopt>-3) { // no calibration missingEx=myMet->MEx; missingEy=myMet->MEy; missingET=myMet->MET; } else if (metopt==-3) { // smeared MC Float_t random; random= gRandom->Rndm(1); Double_t corr_energy=0.4*random-0.2; // random= gRandom->Rndm(1); // Double_t corr_eta=0.2*random-0.1; random= gRandom->Rndm(1); Double_t corr_phi=0.2*random-0.1; missingEx=myMet->MEx; missingEy=myMet->MEy; missingET=myMet->MET; Double_t phi=PhiFun(missingEx,missingEy); phi=phi+corr_phi; missingEx=missingET*cos(phi)*(1.+corr_energy); missingEy=missingET*sin(phi)*(1.+corr_energy); missingET=sqrt(missingEx**2+missingEy**2); // missingET=missingET*(1.+corr_energy); } else if (metopt==1) { // type 1 MET Double_t corr_x=0.; Double_t corr_y=0.; Double_t rumenta_x=0.; Double_t rumenta_y=0.; for (Int_t k=0; kAt(k); if ((notCalJet->PT)<=20) { // measurement of the low Et "garbage" rumenta_x = rumenta_x + notCalJet->Px; rumenta_y = rumenta_y + notCalJet->Py; } if ((notCalJet->PT)>20) { Double_t eta_nocal=notCalJet->Eta; Double_t phi_nocal=notCalJet->Phi; Double_t mindelta=999.; TRootJet *closestCalJet; Int_t safetycheckcount=0; for (Int_t l=0; lAt(l); Double_t eta_cal=calJet->Eta; Double_t phi_cal=calJet->Phi; /* if (deltaR(eta_cal,phi_cal,eta_nocal,phi_nocal)Px) - (notCalJet->Px); corr_y = corr_y + (calJet->Py) - (notCalJet->Py); // corr_x = corr_x + (calJet->PT)*cos(phi_cal) - (notCalJet->PT)*cos(phi_nocal); // corr_y = corr_y + (calJet->PT)*sin(phi_cal)-(notCalJet->PT)*cos(phi_nocal); } } if (safetycheckcount!=1) cout << " TROUBLE! safetycheckcount = " << safetycheckcount << endl; // corr_x = corr_x + (closestCalJet->Px)-(notCalJet->Px); // corr_y = corr_y + (closestCalJet->Py)-(notCalJet->Py); } } Double_t corrmu_x=0; Double_t corrmu_y=0; for (Int_t k=0; kMEx)+corr_x;//-corrmu_x;//+(resc20-1)*rumenta_x; missingEy=(myMet->MEy)+corr_y;//-corrmu_y;//+(resc20-1)*rumenta_y; missingET=sqrt(missingEx**2+missingEy**2); } else if (metopt==2) { // my personal tests Double_t corrmu_x=0; Double_t corrmu_y=0; for (Int_t k=0; kAt(k); if ((notCalJet->PT)<=20) { // measurement of the low Et "garbage" rumenta_x = rumenta_x + notCalJet->Px; rumenta_y = rumenta_y + notCalJet->Py; } if ((notCalJet->PT)>20) { Double_t eta_nocal=notCalJet->Eta; Double_t phi_nocal=notCalJet->Phi; Double_t mindelta=999.; for (Int_t l=0; lAt(l); Double_t eta_cal=calJet->Eta; Double_t phi_cal=calJet->Phi; if ( (fabs(eta_nocal-eta_cal)<0.001) && (fabs(phi_nocal-phi_cal)<0.001) ){ corr_x = corr_x + (calJet->Px) - (notCalJet->Px); corr_y = corr_y + (calJet->Py) - (notCalJet->Py); } } } } missingEx=(myMet->MEx)+corr_x; missingEy=(myMet->MEy)+corr_y; missingET=sqrt(missingEx**2+missingEy**2); } } // end of loop on nMet if (metopt==3) { // true W->nu missingEx=TrueNeutrino->Px(); missingEy=TrueNeutrino->Py(); missingET=TrueNeutrino->Pt(); } TRootMissingET *myMetMC = (TRootMissingET*) metmc->At(0); // for reference missingET_MC=myMetMC->MET; // for reference missingEx_MC=myMetMC->MEx; // for reference // deltaMET=missingET-missingET_MC; // deltaMET_MET=deltaMET/missingET_MC; deltaMET=missingEx-missingEx_MC; deltaMET_MET=deltaMET/missingEx_MC; PtTotX+=missingEx; PtTotY+=missingEy; bool selected = true; if (nLept!=1 || Njets==0) selected=false; if (verbosity==1) cout << "MissingEx: " << missingEx<< " missingEy: " << missingEy << endl; if (verbosity==1) cout << "nGoodElec=" << nGoodElec << ", nGoodMuon=" << nGoodMuon << endl; Int_t nl=0; mtw=-99.; bool isMuon=false; bool isElectron=false; // LEPTON LOOP: Int_t nuindex=0; for (Int_t j=0; jPtLCut) { nl++; } if (PtL>PtLCut && missingET>METCut) { myW = new W(pl_x[j],pl_y[j],pl_z[j],missingEx,missingEy); discriminant=myW->GetDiscriminant(); hDiscriminant->Fill(discriminant); mtw=myW->GetMT(); if (mtw > 100) selected=false; // this cut replaces the old Discriminant >= 0 if (!(myW->GetQuality())) { nosol++; // selected=false; // for safety } else { goodsol++; if (verbosity==1) cout << " soluzioni di Pz per il W: " << myW->GetPZ1() << ", "<< myW->GetPZ2() << endl; } recoWz[1]=myW->GetPZ1(); recoWz[2]=myW->GetPZ2(); } else if (missingET10) { // lepton veto (but I don't want to veto b->l, c->l in jets) selected=false; if (verbosity==1) cout << "too soft lepton Pt, " << PtL << endl; } else { selected=false; cout << "missingET = " << missingET << ", PtL = " << PtL << endl; } if (selected && (fabs(myW->GetPX()) > 9999. || fabs(myW->GetPY()) > 9999. || fabs(myW->GetPZ1()) > 9999. || fabs(myW->GetPZ2()) > 9999. || fabs(myW->GetE1()) > 9999. || fabs(myW->GetE2()) > 9999. || fabs(myW->GetPX()) < 0.001 || fabs(myW->GetPY()) < 0.001 || fabs(myW->GetPZ1()) < 0.001 || fabs(myW->GetPZ2()) < 0.001 || fabs(myW->GetE1()) < 0.001 || fabs(myW->GetE2()) < 0.001)) selected=false; // for safety } hNl->Fill(nl); hMET->Fill(missingET); hDeltaMET->Fill(deltaMET_MET); hMT->Fill(mtw); if (selected) hMTcut->Fill(mtw); // JET COUNTING AND JET VETO: if (verbosity==1) cout << "Njets=" << Njets << endl; Int_t nj=0; Int_t njveto=0; Int_t njbarrel=0; Int_t njendcap=0; Int_t nbj=0; for (Int_t k=0; kAt(k); } else { myJet = (TRootJet*) itConeJets->At(k); }*/ myJet = (TRootBTagJet*) btagJets->At(k); Double_t PtJ=myJet->PT; Double_t etajet=myJet->Eta; if (PtJ>PtJCut) { nj++; if (fabs(etajet)<3) { njbarrel++; } else { njendcap++; } } else if (PtJ20) { //veto region (20 GeV - PtJCut) njveto++; } Double_t btagJet=myJet->discValue; if (PtJ>PtJCut && btagJet>BTagCut) { nbj++; } } hNj->Fill(nj); hNjBarrel->Fill(njbarrel); hNjEndcap->Fill(njendcap); hNbj->Fill(nbj); hNjveto->Fill(njveto); // JET SELECTION: (jets are ordered in Pt, so I use only 1 and 2) Int_t indexHighestBtag=0; if (Njets<2) selected=false; // not needed anymore if (nj!=2) selected=false; // only two HARD (Pt>PtJCut GeV) jets if (njveto>0) selected=false; //jet veto in the SOFT region // if (selected) { if (Njets>1) { TRootBTagJet *myJet1; TRootBTagJet *myJet2; TRootBTagJet *myJet3; myJet1 = (TRootBTagJet*) btagJets->At(0); myJet2 = (TRootBTagJet*) btagJets->At(1); if (njveto>0) myJet3 = (TRootBTagJet*) btagJets->At(2); PtJ1=myJet1->PT; PtJ2=myJet2->PT; EtaJ1=myJet1->Eta; EtaJ2=myJet2->Eta; PhiJ1=myJet1->Phi; PhiJ2=myJet2->Phi; DeltaJ1J2=deltaR(EtaJ1,PhiJ1,EtaJ2,PhiJ2); PtTotX+=(myJet1->Px)+(myJet2->Px); PtTotY+=(myJet1->Py)+(myJet2->Py); PtTotZ+=(myJet1->Pz)+(myJet2->Pz); if (myJet3 && myJet3->PT > 20) PtJ3=myJet3->PT; hPtJ1->Fill(PtJ1); hPtJ2->Fill(PtJ2); if (PtJ2discValue; btag2=myJet2->discValue; hBtag1->Fill(btag1); hBtag2->Fill(btag2); if (btag1Px,myJet1->Py,myJet1->Pz,myJet1->E); LowestTagJet = new TLorentzVector(myJet2->Px,myJet1->Py,myJet2->Pz,myJet2->E); } else { HighestTagJet = new TLorentzVector(myJet2->Px,myJet2->Py,myJet2->Pz,myJet2->E); LowestTagJet = new TLorentzVector(myJet1->Px,myJet1->Py,myJet1->Pz,myJet1->E); } } hHighestBtag->Fill(1.*indexHighestBtag); ht=PtL+missingET+PtJ1+PtJ2+PtJ3; if (Njets>1) { hHt->Fill(ht); hLogHt->Fill(log(ht)); } PtTot=sqrt(PtTotX*PtTotX+PtTotY*PtTotY); hPtTotal->Fill(PtTot); // JET LOOP (highest Pt jets only): if (selected) { for (Int_t k=0; k<2; k++) { // I take the first two jets, independently from the number of allowed jets TRootBTagJet *myJet; /* if (JetAlgo==1) { myJet = (TRootJet*) ktJets->At(k); } else { myJet = (TRootJet*) itConeJets->At(k); }*/ myJet = (TRootBTagJet*) btagJets->At(k); pj_x[k]=myJet->Px; pj_y[k]=myJet->Py; pj_z[k]=myJet->Pz; pj_e[k]=myJet->E; btag[k]=myJet->discValue; // istagged[k]=myJet->isTagged; jetcode[k]=0; Double_t etajet=myJet->Eta; Double_t phijet=myJet->Phi; Double_t PtJ=myJet->PT; for (Int_t q=0; qAt(q); if (myParton->PID==21 || abs(myParton->PID)<6) { // partoni, escluso top if (myParton->PT>20) { // considera solo partoni da un certo pT in su Double_t etaparton=myParton->Eta; Double_t phiparton=myParton->Phi; if (deltaR(etajet,phijet,etaparton,phiparton) < 0.4) { jetcode[k]=myParton->PID; } } } } if (verbosity==1) cout << " Px: " << pj_x[k]<< " Py: " << pj_y[k]<< " Pz: " << pj_z[k]<< " E: " << pj_e[k] << " | PDG code: " << jetcode[k] << " BTag: " << btag[k] << endl; /* if (istagged[k]) { if (verbosity==1) cout << " isTagged: yes" << endl; } else { if (verbosity==1) cout << " isTagged: no" << endl; } */ ///////////////// studi sulla jet charge: jetcharge[k]=0; Int_t nTracksInJet=0; Double_t sumpp=0; if (btag[k]>BTagCut) { // qui andra' messo abs(jetcode[k])==5, per vedere se riesco a distinguere tra +5 e -5 for (Int_t l=0; lAt(l); Double_t etatk=myTrack->Eta; Double_t phitk=myTrack->Phi; if (deltaR(etajet,phijet,etatk,phitk) < 0.5) { nTracksInJet++; // Double_t ptk=myTrack->P; Double_t ptkx=myTrack->Px; Double_t ptky=myTrack->Py; Double_t ptkz=myTrack->Pz; Double_t ptkpar=(ptkx*pj_x[k]+ptky*pj_y[k]+ptkz*pj_z[k])/(pj_x[k]*pj_x[k]+pj_y[k]*pj_y[k]+pj_z[k]*pj_z[k]); Double_t ch=myTrack->Charge; sumpp = sumpp+ptkpar; jetcharge[k]=jetcharge[k]+ptkpar*ch; } } if (sumpp==0) { cout << "sumpp=0, num. of tracks: " << nTracksInJet << endl; } else { jetcharge[k]=jetcharge[k]/sumpp; if (verbosity==1) cout << "jetcharge: " << jetcharge[k] << ", num. of tracks: " << nTracksInJet << endl; if (jetcode[k]==5) { hJetcharge1->Fill(jetcharge[k]); } else if (jetcode[k]==-5) { hJetcharge2->Fill(jetcharge[k]); } else { hJetcharge3->Fill(jetcharge[k]); } } jetchargewl[k]=jetcharge[k]; } ///////////////// assemble a top candidate: if (btag[k]>BTagCut) { cout << "ntop = " << ntop << ", pj_x[k] = " << pj_x[k] << ", pj_y[k] = " << pj_y[k] << ", pj_z[k] = " << pj_z[k] << ", pj_e[k] = " << pj_e[k] << ", btag[k] = " << btag[k] << ", jetcharge[k] = " << jetcharge[k] << ", myW->GetPX() = " << myW->GetPX() << ", myW->GetPY() = " << myW->GetPY() << ", myW->GetPZ1() = " << myW->GetPZ1() << ", myW->GetPZ2() = " << myW->GetPZ2() << ", myW->GetE1() = " << myW->GetE1() << ", myW->GetE2() = " << myW->GetE2() << endl; myT[ntop]->create(neutrino,pj_x[k],pj_y[k],pj_z[k],pj_e[k],btag[k],jetcharge[k],myW->GetPX(),myW->GetPY(),myW->GetPZ1(),myW->GetPZ2(),myW->GetE1(),myW->GetE2()); deltaTrueRecoBfromT[ntop]=deltaR(etajet,phijet,trueBeta,trueBphi); ntop++; } } } // choose "best" Top solution: mtop=-10.; pttop=-10.; topcharge=0.; if (ntop == 1) { mtop=myT[0]->GetMass(); pttop=myT[0]->GetPT(); topcharge=lcharge[0]*myT[0]->GetJetcharge(); whichtop=0; } else if (ntop == 2) { if (bchoice==0) { // random Float_t random = gRandom->Rndm(1); if (random<0.5) { mtop=myT[0]->GetMass(); pttop=myT[0]->GetPT(); topcharge=lcharge[0]*myT[0]->GetJetcharge(); whichtop=0; } else { mtop=myT[1]->GetMass(); pttop=myT[1]->GetPT(); topcharge=lcharge[0]*myT[1]->GetJetcharge(); whichtop=1; } } if (bchoice==1) { mtop=myT[0]->GetMass(); pttop=myT[0]->GetPT(); topcharge=lcharge[0]*myT[0]->GetJetcharge(); whichtop=0; } if (bchoice==2) { mtop=myT[1]->GetMass(); pttop=myT[1]->GetPT(); topcharge=lcharge[0]*myT[1]->GetJetcharge(); whichtop=1; } if (bchoice==3) { // highest b-tag if (myT[0]->GetBtag()>myT[1]->GetBtag()) { mtop=myT[0]->GetMass(); pttop=myT[0]->GetPT(); topcharge=lcharge[0]*myT[0]->GetJetcharge(); whichtop=0; } else { mtop=myT[1]->GetMass(); pttop=myT[1]->GetPT(); topcharge=lcharge[0]*myT[1]->GetJetcharge(); whichtop=1; } } if (bchoice==4) { // jet charge "most opposite" to lepton charge Double_t charge0=-1.*(myT[0]->GetJetcharge())*lcharge[0]; Double_t charge1=-1.*(myT[1]->GetJetcharge())*lcharge[0]; if (charge0>charge1) { mtop=myT[0]->GetMass(); pttop=myT[0]->GetPT(); topcharge=lcharge[0]*myT[0]->GetJetcharge(); whichtop=0; } else { mtop=myT[1]->GetMass(); pttop=myT[1]->GetPT(); topcharge=lcharge[0]*myT[1]->GetJetcharge(); whichtop=1; } } if (bchoice==5) { // highest pt(top) if (myT[0]->GetPT()>myT[1]->GetPT()) { mtop=myT[0]->GetMass(); pttop=myT[0]->GetPT(); topcharge=lcharge[0]*myT[0]->GetJetcharge(); whichtop=0; } else { mtop=myT[1]->GetMass(); pttop=myT[1]->GetPT(); topcharge=lcharge[0]*myT[1]->GetJetcharge(); whichtop=1; } } if (bchoice==6) { // closest m(top) to 175 GeV Double_t nominal=175.; Double_t dm0=fabs(myT[0]->GetMass()-nominal); Double_t dm1=fabs(myT[1]->GetMass()-nominal); if (dm0>dm1) { mtop=myT[0]->GetMass(); pttop=myT[0]->GetPT(); topcharge=lcharge[0]*myT[0]->GetJetcharge(); whichtop=0; } else { mtop=myT[1]->GetMass(); pttop=myT[1]->GetPT(); topcharge=lcharge[0]*myT[1]->GetJetcharge(); whichtop=1; } } if (bchoice==7) { // jet charge opposite to lepton charge, otherwise highest PT(top) Double_t charge0=-1.*(myT[0]->GetJetcharge())*lcharge[0]; Double_t charge1=-1.*(myT[1]->GetJetcharge())*lcharge[0]; if (charge0*charge1<0.) { if (charge0>charge1) { mtop=myT[0]->GetMass(); pttop=myT[0]->GetPT(); topcharge=lcharge[0]*myT[0]->GetJetcharge(); whichtop=0; } else { mtop=myT[1]->GetMass(); pttop=myT[1]->GetPT(); topcharge=lcharge[0]*myT[1]->GetJetcharge(); whichtop=1; } } else { if (myT[0]->GetPT()>myT[1]->GetPT()) { mtop=myT[0]->GetMass(); pttop=myT[0]->GetPT(); topcharge=lcharge[0]*myT[0]->GetJetcharge(); whichtop=0; } else { mtop=myT[1]->GetMass(); pttop=myT[1]->GetPT(); topcharge=lcharge[0]*myT[1]->GetJetcharge(); whichtop=1; } } } } else if (ntop > 2) { selected=false; cout << "#### Too many top candidates: " << ntop << endl; } else if (ntop == 0) { selected=false; cout << "#### No top candidate!" << endl; } nontop=(whichtop+1)%2; RecoTop = new TLorentzVector(myT[whichtop]->GetPX(),myT[whichtop]->GetPY(),myT[whichtop]->GetPZ(),myT[whichtop]->GetE()); if (nTrueTops==1 && ntop>0 && ntop<3) { Double_t etaReco=RecoTop->Eta(); Double_t etaTrue=TrueTop->Eta(); Double_t phiReco=RecoTop->Phi(); Double_t phiTrue=TrueTop->Phi(); Double_t deltaTrueReco=deltaR(etaTrue,phiTrue,etaReco,phiReco); hDelta->Fill(deltaTrueReco); } // Double_t phiTB=acos(PtXY/PtTot); if (selected) { Double_t etb=myT[whichtop]->GetE()+myT[nontop]->GetE(); Double_t pztb=myT[whichtop]->GetPZ()+myT[nontop]->GetPZ(); mtb=sqrt(etb**2-PtTot**2-pztb**2); cout << "mtb = " << mtb << endl; } else { mtb=-1.; } if (mtop>0.) { // if (selection == 2) { // if (ntop==1) hMass->Fill(mtop); // } else { hMass->Fill(mtop); // } hPT->Fill(pttop); } // check neutrino solution quality (distance from truth): nuindex=myT[whichtop]->GetNuSolutionIndex(); if (discriminant>0 && selected && nTrueWs>0 && nuindex>0) { Double_t neutrino_quality=fabs(trueWz-recoWz[nuindex]); hNeutrinoCheck->Fill(neutrino_quality); Int_t nuwrongindex; if (nuindex==1) { nuwrongindex=2; } else if (nuindex==2) { nuwrongindex=1; } else { cout << " -- this shouldn't happen --" << endl; } if (fabs(trueWz-recoWz[nuindex]) < fabs(trueWz-recoWz[nuwrongindex]) ) { countnubest++; } else { countnuworst++; } } // check b-choice quality (distance from truth): if (selected && nTrueBs==1 && ntop==2) { Int_t wrongtop; if (whichtop==0) { wrongtop=1; } else if (whichtop==1) { wrongtop=0; } else { cout << " -- this shouldn't happen --" << endl; } if (verbosity==1) cout << "deltaTrueRecoBfromT[whichtop] = " << deltaTrueRecoBfromT[whichtop] << endl; hDeltaB->Fill(deltaTrueRecoBfromT[whichtop]); Double_t deltaRjet=2; if (deltaTrueRecoBfromT[whichtop]>deltaRjet && deltaTrueRecoBfromT[wrongtop]>deltaRjet) { // the true b from t isn't in one the two considered jets countbout++; } else { if (deltaTrueRecoBfromT[whichtop]mtwHighCut) selected=false; // inutile: c'e' gia' Delta>0 if (mtwhtHighCut) selected=false; Double_t ptj1HighCut=999; Double_t ptj1LowCut=50; Double_t ptj2HighCut=999; Double_t ptj2LowCut=30; if (PtJ1>ptj1HighCut || PtJ1ptj1HighCut || PtJ2EtaJCut || fabs(EtaJ2)>EtaJCut) schannel=false; // entrambi centrali Double_t PtTotCut=150; if (PtTot > PtTotCut) schannel=false; //t-channel: if (ntop!=1) tchannel=false; if (!(max(btag1,btag2)>BTagCut && min(btag1,btag2)Fill(mtop); if (mtopmtopHighCut) { selected=false; schannel=false; tchannel=false; } if (preselected) { nSelected[0]++; hTopcharge->Fill(topcharge); if (topcharge>0) { plusTop++; } else if (topcharge<0) { minusTop++; } if (lcharge[0]>0) { plusLepton++; if (pttop<50) plusLeptonPt1++; if (pttop>50 && pttop<100) plusLeptonPt2++; if (pttop>100 && pttop<150) plusLeptonPt3++; if (pttop>150) plusLeptonPt4++; if (mtop<140) plusLeptonM1++; if (mtop>140 && mtop<200) plusLeptonM2++; if (mtop>200) plusLeptonM3++; } else if (lcharge[0]<0) { minusLepton++; if (pttop<50) minusLeptonPt1++; if (pttop>50 && pttop<100) minusLeptonPt2++; if (pttop>100 && pttop<150) minusLeptonPt3++; if (pttop>150) minusLeptonPt4++; if (mtop<140) minusLeptonM1++; if (mtop>140 && mtop<200) minusLeptonM2++; if (mtop>200) minusLeptonM3++; } } if (schannel) { nSelected[1]++; } if (tchannel) { nSelected[2]++; } //// For polarization studies: if (tchannel) { TVector3 boostvector = -RecoTop.BoostVector(); Lepton->Boost(boostvector); LowestTagJet->Boost(boostvector); TVector3 vlight = LowestTagJet->Vect(); cos_lj=cos(Lepton->Angle(vlight)); cout << "--- cos_lj = " << cos_lj << ", angle = " << Lepton->Angle(vlight) << endl; hCLJ->Fill(cos_lj); } // fill the roottree with the NN variables: if (preselected && savetree) { nntree->Fill(); } // use the NN for further selection: if (preselected) { nPreselected++; if (isElectron) nPreselectedElectrons++; if (isMuon) nPreselectedMuons++; } /* inputs[0]=PtJ1; inputs[1]=PtJ2; inputs[2]=btag1; inputs[3]=btag2; inputs[4]=mtop; inputs[5]=pttop; inputs[6]=mtw; inputs[7]=ht; Double_t nn1=mlp1->Evaluate(0,inputs); Double_t nn2=((mlp2->Evaluate(0,inputs))-3.)*10.; Double_t nn3=mlp3->Evaluate(0,inputs); cout << "mlp1 = " << nn1 << endl; cout << "mlp2 = " << nn2 << endl; cout << "mlp3 = " << nn3 << endl; hMLP1->Fill(nn1); hMLP2->Fill(nn2); hMLP3->Fill(nn3); if (nn1 < nncut1) { schannel=false; } else { nCut1++; } if (nn2 < nncut2) { schannel=false; } else { nCut2++; } if (nn3 < nncut3) { schannel=false; } else { nCut3++; } if (schannel) nFinal++; } */ } //end of event loop cout << "-----------------------" << endl; Double_t f=ratio(nPreselected,nmax); Double_t df=sigmaratio(f,nmax); cout << "Preselected events: " << nPreselected << " (electrons:" << nPreselectedElectrons << ", muons:" << nPreselectedMuons << "), total events: " << nmax << ", ratio: " << f << " +- " << df << endl; /* Double_t f=ratio(nCut1,nmax); Double_t df=sigmaratio(f,nmax); cout << "NN1 cut: " << nCut1 << ", total events: " << nmax << ", ratio: " << f << " +- " << df << endl; Double_t f=ratio(nCut2,nmax); Double_t df=sigmaratio(f,nmax); cout << "NN2 cut: " << nCut2 << ", total events: " << nmax << ", ratio: " << f << " +- " << df << endl; Double_t f=ratio(nCut3,nmax); Double_t df=sigmaratio(f,nmax); cout << "NN3 cut: " << nCut3 << ", total events: " << nmax << ", ratio: " << f << " +- " << df << endl; Double_t f=ratio(nFinal,nmax); Double_t df=sigmaratio(f,nmax); cout << "Selected events in NN analysis: " << nFinal << ", total events: " << nmax << ", ratio: " << f << " +- " << df << endl; */ Double_t f=ratio(nSelected[1],nmax); Double_t df=sigmaratio(f,nmax); cout << "Selected events in cut analysis: " << nSelected[1] << ", total events: " << nmax << ", ratio: " << f << " +- " << df << endl; cout << "-----------------------" << endl; cout << "# Charge asymmetry, integrated:" << endl; Double_t f=ratio(plusLepton-minusLepton,abs(minusLepton)); Double_t df=sigmaratio(f,abs(minusLepton)); cout << "Positive lepton: " << plusLepton << ", negative lepton: " << minusLepton << " (excess: " << f << "+-" << df << ")" << endl; /* Double_t f=ratio(nSelected[0],nmax); Double_t df=sigmaratio(f,nmax); cout << "Selected events (s+t): " << nSelected[0] << ", total events: " << nmax << ", ratio: " << f << " +- " << df << endl; // s-ch.: Double_t f=ratio(nSelected[1],nmax); Double_t df=sigmaratio(f,nmax); cout << "Selected events (s-channel): " << nSelected[1] << ", total events: " << nmax << ", ratio: " << f << " +- " << df << endl; // t-ch.: Double_t f=ratio(nSelected[2],nmax); Double_t df=sigmaratio(f,nmax); cout << "Selected events (t-channel): " << nSelected[2] << ", total events: " << nmax << ", ratio: " << f << " +- " << df << endl; ///// the following is for generic (s+t) selection: cout << "###" << endl; Double_t f=ratio(minusTop,nSelected[0]); Double_t df=sigmaratio(f,nSelected[0]); cout << "Concord lj charge: " << plusTop << ", discord lj charge: " << minusTop << " (" << f << "+-" << df << ")" << endl; cout << "###" << endl; */ /* cout << "# Charge asymmetry, Pt(top)<50 GeV:" << endl; Double_t f=ratio(plusLeptonPt1-minusLeptonPt1,abs(minusLeptonPt1)); Double_t df=sigmaratio(f,abs(minusLeptonPt1)); cout << "Positive lepton: " << plusLeptonPt1 << ", negative lepton: " << minusLeptonPt1 << " (excess: " << f << "+-" << df << ")" << endl; cout << "# Charge asymmetry, 50150 GeV:" << endl; Double_t f=ratio(plusLeptonPt4-minusLeptonPt4,abs(minusLeptonPt4)); Double_t df=sigmaratio(f,abs(minusLeptonPt4)); cout << "Positive lepton: " << plusLeptonPt4 << ", negative lepton: " << minusLeptonPt4 << " (excess: " << f << "+-" << df << ")" << endl; cout << "# Charge asymmetry, M(top)<140 GeV:" << endl; Double_t f=ratio(plusLeptonM1-minusLeptonM1,abs(minusLeptonM1)); Double_t df=sigmaratio(f,abs(minusLeptonM1)); cout << "Positive lepton: " << plusLeptonM1 << ", negative lepton: " << minusLeptonM1 << " (excess: " << f << "+-" << df << ")" << endl; cout << "# Charge asymmetry, 140200 GeV:" << endl; Double_t f=ratio(plusLeptonM3-minusLeptonM3,abs(minusLeptonM3)); Double_t df=sigmaratio(f,abs(minusLeptonM3)); cout << "Positive lepton: " << plusLeptonM3 << ", negative lepton: " << minusLeptonM3 << " (excess: " << f << "+-" << df << ")" << endl; */ cout << "-----------------------" << endl; Double_t f=ratio(nosol,nosol+goodsol); Double_t df=sigmaratio(f,nosol+goodsol); cout << "Discriminant < 0: " << nosol << " (" << f*100 << "+-" << df*100 << " %), Discriminant > 0: " << goodsol << endl; cout << "-----------------------" << endl; Double_t f=ratio(countnubest,countnubest+countnuworst); Double_t df=sigmaratio(f,countnubest+countnuworst); cout << "Our neutrino choice criterium gave the best answer " << countnubest << " times against " << countnuworst << " (" << f*100 << "+-" << df*100 << " %)" << endl; cout << "-----------------------" << endl; Double_t f=ratio(countbbest,countbbest+countbworst+countbout); Double_t df=sigmaratio(f,countbbest+countbworst+countbout); cout << "Our b-jet choice criterium gave the best answer " << countbbest << " times against " << countbworst << " (" << f*100 << "+-" << df*100 << " %); no b-jet included the true b in " << countbout << " cases" << endl; cout << "-----------------------" << endl; ////////// Write on file: // (http://www.disp.uniroma2.it/users/buttarazzi/lezionisvolte/CPP3/CPP3.HTM) TString outname("out_"); outname+=nomesample; outname+="_"; outname+=metopt; outname+="_"; outname+=neutrino; outname+="_"; outname+=bchoice; outname+="_"; outname+=whichlepton; outname+=".txt"; fstream f1; f1.open(outname,ios::out); f1 << nosol << " " << goodsol; f1 << endl; f1 << countnubest << " " << countnuworst; f1 << endl; f1 << countbbest << " " << countbworst; f1 << endl; f1 << nmax; f1.close(); cout << outname << " saved." << endl; ////////// Visualization: Char_t temp[80]; Char_t idlep[2]; if (whichlepton==1) {sprintf(idlep,"e");} else if (whichlepton==2) {sprintf(idlep,"m");} else {sprintf(idlep,"em");} if (!plots) goto printfile; bool escape=false; while (!escape) { cout << endl << "Choose the histogram to plot:" << endl; cout << "0: none (exit)" << endl; cout << "1: W transverse mass" << endl; cout << "2: jet charge" << endl; cout << "3: jet charge times lepton charge" << endl; cout << "4: top Pt" << endl; cout << "5: top reco mass" << endl; cout << "6: number of jets" << endl; cout << "7: angular distrib. of lepton wrt light jet in top r.f. " << endl; cout << "8: discriminant of the W mass constraint equation" << endl; cout << "9: number of electrons" << endl; cout << "10: angular distance between true and reco top" << endl; cout << "11: angular distance between b(<-t) and reco b" << endl; cout << "12: electron isolation" << endl; cout << "13: highest btag jet" << endl; cout << "14: W true-reco z-distance" << endl; cout << "15: Ht" << endl; cout << "16: Pt jets #1 and #2" << endl; cout << "17: MET" << endl; cout << "18: btag jets #1 and #2" << endl; cout << "19: Pt of top+bottom system" << endl; // cout << "20: NN output 1 (s-ch. vs ttbar)" << endl; // cout << "21: NN output 2 (s-ch. vs Wbbbar)" << endl; // cout << "22: NN output 3 (s-ch. vs t-ch.)" << endl; Int_t scelta; cin >> scelta; Int_t delta=20*(scelta-1); if (scelta==0) escape=true; if (scelta==1) { cMT=new TCanvas("cMT","Mt W",0+delta,0+delta,600,400); cMT->cd(); pad1 = new TPad("pad1","before cut",0.05,0.50,0.95,0.95,21); pad2 = new TPad("pad2","after cut",0.05,0.05,0.95,0.45,21); pad1->Draw(); pad2->Draw(); pad1->cd(); hMT->DrawCopy(); pad2->cd(); hMTcut->DrawCopy(); sprintf(temp,"mtw_%d_%s_%s.ps",metopt,idlep,nomesample); cMT->Print(temp); } if (scelta==2) { cJetcharge=new TCanvas("cJetcharge","Jet charge",0+delta,0+delta,600,400); cJetcharge->cd(); Double_t scale; scale = 1/hJetcharge1->Integral(); hJetcharge1->Scale(scale); scale = 1/hJetcharge2->Integral(); hJetcharge2->Scale(scale); scale = 1/hJetcharge3->Integral(); hJetcharge3->Scale(scale); hJetcharge1->SetLineWidth(2); hJetcharge2->SetLineWidth(2); hJetcharge3->SetLineWidth(2); hJetcharge1->SetLineColor(kRed); hJetcharge1->DrawCopy(); hJetcharge2->SetLineColor(kBlue); hJetcharge2->DrawCopy("same"); hJetcharge3->SetLineColor(kBlack); // hJetcharge3->DrawCopy("same"); Double_t xmin=-0.99; Double_t xmax=-0.6; Double_t ymax=hJetcharge1->GetMaximum(); TLegend *legend = new TLegend(xmin,0.75*ymax,xmax,0.95*ymax,"flavour","br"); legend->AddEntry(hJetcharge1,"b","l"); legend->AddEntry(hJetcharge2,"bbar","l"); legend->Draw(); sprintf(temp,"jetcharge_%s.ps",nomesample); cJetcharge->Print(temp); } if (scelta==3) { cTopcharge=new TCanvas("cTopcharge","Top charge",0+delta,0+delta,600,400); cTopcharge->cd(); hTopcharge->DrawCopy(); sprintf(temp,"topcharge_%s.ps",nomesample); cTopcharge->Print(temp); } if (scelta==4) { cPT=new TCanvas("cPT","PT top",0+delta,0+delta,600,400); cPT->cd(); hPT->DrawCopy(); sprintf(temp,"pttop_%d_%d_%s_%s.ps",metopt,neutrino,idlep,nomesample); cPT->Print(temp); } if (scelta==5) { cMass=new TCanvas("cMass","Mass top",0+delta,0+delta,600,400); cMass->cd(); pad1 = new TPad("pad1","before cut",0.05,0.50,0.95,0.95,21); pad2 = new TPad("pad2","after cut",0.05,0.05,0.95,0.45,21); pad1->Draw(); pad2->Draw(); pad1->cd(); hMass->DrawCopy(); pad2->cd(); hMasscut->DrawCopy(); sprintf(temp,"mtop_%d_%d_%s_%s.ps",metopt,neutrino,idlep,nomesample); cMass->Print(temp); } if (scelta==6) { cJETS=new TCanvas("cJETS","Jets",0+delta,0+delta,600,400); cJETS->cd(); hNgj->SetLineWidth(2); hNuncalib->SetLineWidth(2); hNcalib->SetLineWidth(2); hNgj->DrawCopy(); hNuncalib->SetLineColor(kRed); hNuncalib->DrawCopy("same"); hNcalib->SetLineColor(kBlue); hNcalib->DrawCopy("same"); Double_t xmin=6.5; Double_t xmax=10; Double_t ymax=hNgj->GetMaximum(); TLegend *legend = new TLegend(xmin,0.75*ymax,xmax,0.95*ymax,"","br"); legend->AddEntry(hNgj,"generated jets","l"); legend->AddEntry(hNuncalib,"reco. uncalib. jets","l"); legend->AddEntry(hNcalib,"reco. calib. jets","l"); legend->Draw(); sprintf(temp,"njets_%s.ps",nomesample); cJETS->Print(temp); } if (scelta==7) { cCLJ=new TCanvas("cCLJ","cos alpha_LJ",0+delta,0+delta,600,400); cCLJ->cd(); hCLJ->DrawCopy(); sprintf(temp,"coslj_%d_%s_%s.ps",neutrino,idlep,nomesample); cCLJ->Print(temp); } if (scelta==8) { cDiscriminant=new TCanvas("cDiscriminant","Discriminant",0+delta,0+delta,600,400); cDiscriminant->cd(); hDiscriminant->DrawCopy(); sprintf(temp,"discriminant_%d_%s_%s.ps",metopt,idlep,nomesample); cDiscriminant->Print(temp); } if (scelta==9) { cNe=new TCanvas("cNe","N. electrons",0+delta,0+delta,600,400); cNe->cd(); hNe->DrawCopy(); } if (scelta==10) { cDelta=new TCanvas("cDelta","Delta top true / top reco",0+delta,0+delta,600,400); cDelta->cd(); hDelta->DrawCopy(); sprintf(temp,"delta_%d_%d_%s_%s.ps",neutrino,bchoice,idlep,nomesample); cDelta->Print(temp); } if (scelta==11) { cDeltaB=new TCanvas("cDeltaB","Delta b true / b reco",0+delta,0+delta,600,400); cDeltaB->cd(); hDeltaB->DrawCopy(); sprintf(temp,"deltab_%d_%d_%s_%s.ps",neutrino,bchoice,idlep,nomesample); cDeltaB->Print(temp); } if (scelta==12) { cEleIso=new TCanvas("cEleIso","Electron isolation",0+delta,0+delta,600,400); cEleIso->cd(); hEleIso->DrawCopy(); sprintf(temp,"eleiso_%s.ps",nomesample); cEleIso->Print(temp); } if (scelta==13) { cHighestBtag=new TCanvas("cHighestBtag","Highest btag jet",0+delta,0+delta,600,400); cHighestBtag->cd(); hHighestBtag->DrawCopy(); sprintf(temp,"highestbtag_%s.ps",nomesample); cHighestBtag->Print(temp); } if (scelta==14) { cNeutrinoCheck=new TCanvas("cNeutrinoCheck","W true-reco z-distance",0+delta,0+delta,600,400); cNeutrinoCheck->cd(); hNeutrinoCheck->DrawCopy(); sprintf(temp,"nucheck_%s.ps",nomesample); cNeutrinoCheck->Print(temp); } if (scelta==15) { cHt=new TCanvas("cHt","Ht",0+delta,0+delta,600,400); cHt->cd(); pad1 = new TPad("pad1","Ht",0.05,0.50,0.95,0.95,21); pad2 = new TPad("pad2","log(Ht)",0.05,0.05,0.95,0.45,21); pad1->Draw(); pad2->Draw(); pad1->cd(); hHt->DrawCopy(); pad2->cd(); hLogHt->DrawCopy(); sprintf(temp,"ht_%d_%s.ps",metopt,nomesample); cHt->Print(temp); } if (scelta==16) { cPTjets=new TCanvas("cPTjets","Pt jets",0+delta,0+delta,600,400); cPTjets->cd(); pad1 = new TPad("pad1","Pt jet 1",0.05,0.50,0.95,0.95,21); pad2 = new TPad("pad2","Pt jet 2",0.05,0.05,0.95,0.45,21); pad1->Draw(); pad2->Draw(); pad1->cd(); hPtJ1->DrawCopy(); pad2->cd(); hPtJ2->DrawCopy(); sprintf(temp,"ptjets_%d_%s.ps",metopt,nomesample); cPTjets->Print(temp); } if (scelta==17) { cMET=new TCanvas("cMET","Missing Et",0+delta,0+delta,600,400); cMET->cd(); pad1 = new TPad("pad1","MET",0.05,0.50,0.95,0.95,21); pad2 = new TPad("pad2","deltaMET/MET",0.05,0.05,0.95,0.45,21); pad1->Draw(); pad2->Draw(); pad1->cd(); hMET->DrawCopy(); pad2->cd(); hDeltaMET->DrawCopy(); sprintf(temp,"met_%d_%s.ps",metopt,nomesample); cMET->Print(temp); } if (scelta==18) { cBtag=new TCanvas("cBtag","b-tag",0+delta,0+delta,600,400); cBtag->cd(); pad1 = new TPad("pad1","b-tag jet 1",0.05,0.50,0.95,0.95,21); pad2 = new TPad("pad2","b-tag jet 2",0.05,0.05,0.95,0.45,21); pad1->Draw(); pad2->Draw(); pad1->cd(); hBtag1->DrawCopy(); pad2->cd(); hBtag2->DrawCopy(); sprintf(temp,"btag_%d_%s.ps",metopt,nomesample); cBtag->Print(temp); } if (scelta==19) { cPtTotal=new TCanvas("cPtTotal","Pt of top+bottom system",0+delta,0+delta,600,400); cPtTotal->cd(); hPtTotal->DrawCopy(); sprintf(temp,"pttotal_%d_%s_%s.ps",metopt,idlep,nomesample); cPtTotal->Print(temp); } /* if (scelta==20) { cMLP1=new TCanvas("cMLP1","NN output 1",0+delta,0+delta,600,400); cMLP1->cd(); hMLP1->DrawCopy(); sprintf(temp,"mlp1_%d_%s_%s.ps",metopt,idlep,nomesample); cMLP1->Print(temp); } if (scelta==21) { cMLP2=new TCanvas("cMLP2","NN output 2",0+delta,0+delta,600,400); cMLP2->cd(); hMLP2->DrawCopy(); sprintf(temp,"mlp2_%d_%s_%s.ps",metopt,idlep,nomesample); cMLP2->Print(temp); } if (scelta==22) { cMLP3=new TCanvas("cMLP3","NN output 3",0+delta,0+delta,600,400); cMLP3->cd(); hMLP3->DrawCopy(); sprintf(temp,"mlp3_%d_%s_%s.ps",metopt,idlep,nomesample); cMLP3->Print(temp); } */ } // Write histos in a new rootfile: printfile: sprintf(temp,"histo_%d_%d_%d_%s_%s.root",metopt,neutrino,bchoice,idlep,nomesample); TFile fout(temp,"RECREATE"); hMT->Write(); hMTcut->Write(); hMET->Write(); hDeltaMET->Write(); hMass->Write(); hMasscut->Write(); hHt->Write(); hLogHt->Write(); hNj->Write(); hNjveto->Write(); hNjBarrel->Write(); hNjEndcap->Write(); hNbj->Write(); hNgj->Write(); hNuncalib->Write(); hNcalib->Write(); hNl->Write(); hNe->Write(); hDelta->Write(); hDeltaB->Write(); hCLJ->Write(); hEleIso->Write(); hHighestBtag->Write(); hBtag1->Write(); hBtag2->Write(); hDiscriminant->Write(); hNeutrinoCheck->Write(); hPtJ1->Write(); hPtJ2->Write(); hPtTotal->Write(); hMLP1->Write(); hMLP2->Write(); hMLP3->Write(); fout.Write(); cout << "### histogram file " << temp << " created." << endl; // write the roottree with the NN variables: if (savetree) { sprintf(temp,"nn_%d_%d_%d_%s_%s.root",metopt,neutrino,bchoice,idlep,nomesample); TFile foutnn(temp,"RECREATE"); nntree->Write(); foutnn.Write(); cout << "$$$ NN input file " << temp << " created." << endl; } // Clean histos: delete hpmux; delete hpmuy; delete hpmuz; delete hpWx; delete hpWy; delete hpWz1; delete hpWz2; delete hMT; delete hMTcut; delete hMET; delete hDeltaMET; delete hCLJ; delete hPT; delete hPtJ1; delete hPtJ2; delete hPtTotal; delete hPhiTB; delete hMass; delete hMasscut; delete hTopcharge; delete hJetcharge1; delete hJetcharge2; delete hJetcharge3; delete hJetchargewl1; delete hJetchargewl2; delete hJetchargewl3; delete hHt; delete hLogHt; delete hNj; delete hNjveto; delete hNjBarrel; delete hNjEndcap; delete hNbj; delete hNgj; delete hNuncalib; delete hNcalib; delete hNl; delete hNe; delete hDelta; delete hDeltaB; delete hDiscriminant; delete hEleIso; delete hHighestBtag; delete hBtag1; delete hBtag2; delete hNeutrinoCheck; delete hMLP1; delete hMLP2; delete hMLP3; } Double_t max(Double_t a, Double_t b) { Double_t res=a; if (b>a) res=b; return res; } Double_t max(Double_t a, Double_t b, Double_t c, Double_t d, Double_t e) { Double_t res; if (a>b && a>c && a>d && a>e) res=a; if (b>a && b>c && b>d && b>e) res=b; if (c>a && c>b && c>d && c>e) res=c; if (d>a && d>b && d>c && d>e) res=d; if (e>a && e>b && e>c && e>d) res=e; return res; } Double_t min(Double_t a, Double_t b) { Double_t res=a; if (b (2*TMath::Pi()-delta)) delta=(2*TMath::Pi()-delta); return delta; } Double_t PhiFun(Double_t x,Double_t y) // in CMSJET: real function phifun { Double_t phi; phi=atan2(y,x); if (phi<0) phi=phi+2*TMath::Pi(); // cout << " phi = " << phi << endl; return phi; } Double_t ratio(Int_t n1,Int_t n2) { Double_t x1=1.*n1; Double_t x2=1.*n2; Double_t r; if (x2!=0){ r=x1/x2; } else { cout << "**** dividing by zero!!!!!!!" << endl; r=-99999.; } return r; } Double_t sigmaratio(Double_t r,Int_t n2) { Double_t x2=1.*n2; Double_t sigma; if (x2!=0){ sigma=sqrt(r*(1.-r)/x2); } else { cout << "**** dividing by zero!!!!!!!" << endl; sigma=-99999.; } return sigma; } //////////////////// class Top { public: Top(); // Top(unsigned int str, Double_t jx, Double_t jy, Double_t jz, Double_t je, Double_t bt, Double_t jetch, Double_t wx, Double_t wy, Double_t wz1, Double_t wz2, Double_t we1, Double_t we2); ~Top(); void create(unsigned int str, Double_t jx, Double_t jy, Double_t jz, Double_t je, Double_t bt, Double_t jetch, Double_t wx, Double_t wy, Double_t wz1, Double_t wz2, Double_t we1, Double_t we2); Int_t GetNuSolutionIndex() const {return nusol;} Double_t GetPX() const {return x;} Double_t GetPY() const {return y;} Double_t GetPZ() const {return z;} Double_t GetPZ1() const {return z1;} Double_t GetPZ2() const {return z2;} Double_t GetE() const {return e;} Double_t GetE1() const {return e1;} Double_t GetE2() const {return e2;} Double_t GetM1() const {return m1;} Double_t GetM2() const {return m2;} Double_t GetPT() const {return pt;} Double_t GetBtag() const {return btag;} Double_t GetJetcharge() const {return jetcharge;} bool GetExist() const {return exist;} Double_t GetMass(); void swap(Double_t, Double_t); Double_t GetJetPX() const {return jx;} Double_t GetJetPY() const {return jy;} Double_t GetJetPZ() const {return jz;} Double_t GetJetE() const {return je;} private: Int_t nusol; Double_t wwz1; Double_t wwz2; Double_t x; Double_t y; Double_t z; Double_t z1; Double_t z2; Double_t e; Double_t e1; Double_t e2; Double_t pt; Double_t m1; Double_t m2; Double_t btag; Double_t jetcharge; Double_t mass; unsigned int strategy; bool exist; }; Top::Top() { exist=false; // the constructor is called for every event, but the top starts to exist only when the "create" method is used. } //Top::Top(unsigned int str, Double_t jx, Double_t jy, Double_t jz, Double_t je, Double_t bt, Double_t jetch, Double_t wx, Double_t wy, Double_t wz1, Double_t wz2, Double_t we1, Double_t we2) void Top::create(unsigned int str, Double_t jx, Double_t jy, Double_t jz, Double_t je, Double_t bt, Double_t jetch, Double_t wx, Double_t wy, Double_t wz1, Double_t wz2, Double_t we1, Double_t we2) { // cout << "*** Creato t ***" << endl; exist=true; strategy=str; btag=bt; jetcharge=jetch; wwz1=wz1; wwz2=wz2; x=jx+wx; y=jy+wy; z=0.; z1=jz+wz1; z2=jz+wz2; e=0.; e1=je+we1; e2=je+we2; pt=sqrt(x*x+y*y); m1=sqrt(e1*e1-x*x-y*y-z1*z1); m2=sqrt(e2*e2-x*x-y*y-z2*z2); nusol=0; } Double_t Top::GetMass() { Double_t Mtrue=175.; mass=-1.; if (strategy==0) { // random Float_t random = gRandom->Rndm(1); // cout << " random =" << random << endl; if (random<0.5) { nusol=1; mass=m1; z=z1; e=e1; } else { nusol=2; mass=m2; z=z2; e=e2; } } if (strategy==1) { // 1st neutrino solution nusol=1; mass=m1; z=z1; e=e1; } if (strategy==2) { // 2nd neutrino solution nusol=2; mass=m2; z=z2; e=e2; } if (strategy==3) { // closest to Mtrue if (fabs(m1-Mtrue)