PWGPP/analysisQA/processJETriggerQA.C

   1 /***************************************************
   2 processJETriggerQA:
   3 To procees JE Triggere QA wagon's output
   4
   5 ****************************************************/
   6
   7
   8
   9 const Float_t ptmin =  0. ; //lower cutoff of jet pt spectrum
  10
  11 void processJETriggerQA(TString strFileIn = "AnalysisResults.root",
  12                         TString suftype="eps",
  13                         Float_t jetR         = 0.2,
  14                         Float_t minTrkPT     = 0.15,
  15                         Float_t minClusterET = 0.3,
  16                         Int_t run            = 0,
  17                         const char* outfile   = "JETriggerQA_outfile.root"
  18                         ){
  19
  20   gStyle->SetOptStat(0);
  21   gStyle->SetOptTitle(0);
  22
  23   TString prefix = "fig_je_TriggerQA_";
  24
  25   TFile * f1 = TFile::Open(strFileIn.Data());
  26
  27   //Load histogram list
  28   TList *histList = 0x0;
  29   histList =  (TList*)f1->Get(Form("TriggerQA_Jet_AKTFullR%02d0_PicoTracks_pT%04d_CaloClustersCorr_ET%04d_pt_scheme_Jet_AKTChargedR%02d0_PicoTracks_pT%04d_CaloClustersCorr_ET%04d_pt_scheme_TC/TriggerQA_Jet_AKTFullR%02d0_PicoTracks_pT%04d_CaloClustersCorr_ET%04d_pt_scheme_Jet_AKTChargedR%02d0_PicoTracks_pT%04d_CaloClustersCorr_ET%04d_pt_scheme_TC",
  30        TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),
  31        TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),
  32        TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),
  33        TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000)));
  34
  35   //Load a second histogram list in order to incorporate results from the Rho Task as well
  36   TList *histListRho = 0x0;
  37   histListRho = (TList*)f1->Get(Form("Rho_Jet_KTChargedR%02d0_PicoTracks_pT%04d_CaloClustersCorr_ET%04d_pt_scheme_TPC_histos",
  38        TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000)));
  39
  40   //---------------------------------------------------------------------------------------------------
  41   //       jet histograms
  42   //---------------------------------------------------------------------------------------------------
  43
  44   const Int_t kJetType = 2;
  45   TString suffix [kJetType] =  {"Charged","Full"};
  46   TH3F *h3PtEtaPhiJet[kJetType];
  47   TH1F *hPtJet[kJetType];
  48   TH2F *hEtaPhiJet[kJetType];
  49   TH2F *hRhoCent;
  50
  51   TH1F *hNEventSel = histList->FindObject("fhNEvents");
  52   Float_t nEvents  = hNEventSel->GetBinContent(2);
  53
  54   for(Int_t itype = 0; itype < kJetType; itype++){
  55
  56      h3PtEtaPhiJet[itype] = (TH3F*)  histList->FindObject(Form("fh3PtEtaPhiJet%s",suffix[itype].Data()));
  57      if(! h3PtEtaPhiJet[itype]) continue;
  58
  59      //jet pt spectra
  60      Int_t binMin = 1;
  61      if(ptmin>0.) binMin = h3PtEtaPhiJet[itype]->GetXaxis()->FindBin(ptmin+0.00001);
  62      h3PtEtaPhiJet[itype]->GetXaxis()->SetRange(binMin, h3PtEtaPhiJet[itype]->GetNbinsX());
  63
  64      hPtJet[itype] = (TH1F*) h3PtEtaPhiJet[itype]->Project3D("x");
  65      hPtJet[itype]->Scale(1./nEvents,"width");
  66      SetHist((TH1F*) hPtJet[itype],"p_{T,corr}^{jet} (GeV)","1/N_{evt} dN/dp_{T,corr}^{jet} (GeV^{-1})");
  67      hPtJet[itype]->SetName(Form("hPtJet%s",suffix[itype].Data()));
  68
  69      //eta versus phi
  70      hEtaPhiJet[itype] = (TH2F*) h3PtEtaPhiJet[itype]->Project3D("yz");
  71      SetHist((TH1F*) hEtaPhiJet[itype],"#varphi^{jet} (rad)","#eta^{jet}");
  72      hEtaPhiJet[itype]->SetName(Form("hEtaPhiJet%s",suffix[itype].Data()));
  73   }
  74
  75   //rho versus centrality
  76   hRhoCent = (TH2F*) histListRho->FindObject("fHistRhovsCent");
  77   if(hRhoCent) {
  78      SetHist((TH1F*) hRhoCent, "Centrality (%)", "#rho (GeV/c*rad^{-1})");
  79      hRhoCent->SetName("hRhoCent");
  80   }
  81
  82   //______________
  83   //Draw histograms
  84
  85   TCanvas *c[100];
  86   TH1F *frame[100];
  87   Int_t nCan = 0;
  88   TLegend *leg;
  89
  90   for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
  91
  92      //draw pt spectrum
  93      if(!hPtJet[itype]) continue;
  94      c[nCan] = new TCanvas(Form("c%d",nCan),Form("c%d: Pt %s jets",nCan,suffix[itype].Data()),600,450);
  95      SetCanvas((TCanvas*) c[nCan]);
  96      c[nCan]->SetLogy();
  97
  98      frame[nCan] = gPad->DrawFrame(hPtJet[itype]->GetBinLowEdge(1),
  99                                    1e-7,
 100                                    hPtJet[itype]->GetBinLowEdge(hPtJet[itype]->GetNbinsX()+1),
 101                                    hPtJet[itype]->GetBinContent(hPtJet[itype]->GetMaximumBin())*2.);
 102
 103      SetHist((TH1F*) frame[nCan],hPtJet[itype]->GetXaxis()->GetTitle(),hPtJet[itype]->GetYaxis()->GetTitle());
 104
 105
 106      hPtJet[itype]->DrawCopy("same");
 107
 108      leg = new TLegend(0.35,0.5,0.88,0.88);
 109      SetLeg(leg);
 110
 111      TString txt =  Form("%s jets AKT R=%.1f",suffix[itype].Data(),jetR);
 112      if(run>0) txt += Form(" run:%d",run);
 113                               ;
 114      leg->AddEntry((TObject*) 0, txt.Data(),"");
 115      leg->AddEntry((TObject*) 0, Form("p_{T,trk}> %d MeV",TMath::Nint(minTrkPT*1000)),"");
 116      if(itype==1) leg->AddEntry((TObject*) 0, Form("E_{T}>%d MeV",TMath::Nint(minClusterET*1000)),"");
 117      leg->AddEntry((TObject*) 0, Form("#it{N}_{events} = %.0f",nEvents),"");
 118      leg->Draw();
 119
 120      c[nCan]->SaveAs(Form("%s_Pt_AKT%02d_pT%04d_ET%04d_Run%d_%s.%s",prefix.Data(),TMath::Nint(jetR*10),
 121                           TMath::Nint(minTrkPT*1000),TMath::Nint(minClusterET*1000), run, suffix[itype].Data(), suftype.Data()));
 122
 123      nCan++;
 124
 125      //_________________
 126      //draw eta versus phi
 127
 128      if(!hEtaPhiJet[itype]) continue;
 129      c[nCan] = new TCanvas(Form("c%d",nCan),Form("c%d: eta-phi %s jets",nCan,suffix[itype].Data()),600,450);
 130      SetCanvas((TCanvas*) c[nCan]);
 131      c[nCan]->SetRightMargin(0.15);
 132
 133      frame[nCan] = gPad->DrawFrame(hEtaPhiJet[itype]->GetXaxis()->GetBinLowEdge(1),
 134                                    hEtaPhiJet[itype]->GetYaxis()->GetBinLowEdge(1),
 135                                    hEtaPhiJet[itype]->GetXaxis()->GetBinLowEdge(hEtaPhiJet[itype]->GetNbinsX()),
 136                                    hEtaPhiJet[itype]->GetYaxis()->GetBinLowEdge(hEtaPhiJet[itype]->GetNbinsY()));
 137
 138      SetHist((TH1F*) frame[nCan],hEtaPhiJet[itype]->GetXaxis()->GetTitle(),hEtaPhiJet[itype]->GetYaxis()->GetTitle());
 139
 140
 141      hEtaPhiJet[itype]->DrawCopy("same,colz");
 142
 143      leg = new TLegend(0.35,0.5,0.88,0.88);
 144      SetLeg(leg);
 145
 146      leg->AddEntry((TObject*) 0, txt.Data(),"");
 147      leg->AddEntry((TObject*) 0, Form("p_{T,trk}> %d MeV",TMath::Nint(minTrkPT*1000)),"");
 148      if(itype==1) leg->AddEntry((TObject*) 0, Form("E_{T}>%d MeV",TMath::Nint(minClusterET*1000)),"");
 149      leg->AddEntry((TObject*) 0, Form("#it{N}_{events} = %.0f",nEvents),"");
 150      leg->AddEntry((TObject*) 0, Form("p_{T,corr}^{jet} > %.1f GeV", ptmin),"");
 151      leg->Draw();
 152
 153      c[nCan]->SaveAs(Form("%s_EtaPhi_AKT%02d_pT%04d_ET%04d_Run%d_%s.%s",prefix.Data(),TMath::Nint(jetR*10),
 154                           TMath::Nint(minTrkPT*1000),TMath::Nint(minClusterET*1000), run, suffix[itype].Data(), suftype.Data()));
 155
 156      nCan++;
 157
 158   }//end of the loop over charged and full jets
 159
 160   //draw rho versus centrality
 161   if(hRhoCent) {
 162      c[nCan] = new TCanvas(Form("c%d",nCan),Form("c%d: Rho-Cent",nCan),600,450);
 163      SetCanvas((TCanvas*) c[nCan]);
 164      c[nCan]->SetRightMargin(0.15);
 165      c[nCan]->SetLogz();
 166
 167      frame[nCan] = gPad->DrawFrame(hRhoCent->GetXaxis()->GetBinLowEdge(1),
 168                                    hRhoCent->GetYaxis()->GetBinLowEdge(1),
 169                                    hRhoCent->GetXaxis()->GetBinLowEdge(hRhoCent->GetNbinsX()),
 170                                    hRhoCent->GetYaxis()->GetBinLowEdge(hRhoCent->GetNbinsY()));
 171
 172      SetHist((TH1F*) frame[nCan],hRhoCent->GetXaxis()->GetTitle(),hRhoCent->GetYaxis()->GetTitle());
 173
 174      hRhoCent->DrawCopy("colz");
 175
 176      leg = new TLegend(0.35,0.5,0.88,0.88);
 177      SetLeg(leg);
 178
 179      TString txt2 = Form("Jets KT R=%.1f", jetR);
 180      if(run>0) txt2 += Form(" run:%d",run);
 181
 182      leg->AddEntry((TObject*) 0, txt2.Data(),"");
 183      leg->AddEntry((TObject*) 0, Form("#it{N}_{events} = %.0f",nEvents),"");
 184      leg->Draw();
 185
 186      c[nCan]->SaveAs(Form("%s_RhoCent_AKT%02d_pT%04d_ET%04d_Run%d.%s",prefix.Data(),TMath::Nint(jetR*10),
 187                           TMath::Nint(minTrkPT*1000),TMath::Nint(minClusterET*1000), run, suftype.Data()));
 188
 189      nCan++;
 190   }
 191
 192   //---------------------------------------------------------------------------------------------------
 193   //                       WRITE OUTPUT TO ROOT FILE
 194   //---------------------------------------------------------------------------------------------------
 195
 196
 197   /* Standalone output
 198     TFile *histOut = new TFile(Form("%s_AKT%02d_pT%04d_ET%04d_jetPtMin%.1f_Run%d.root",prefix.Data(),
 199     TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),ptmin,run),"RECREATE");
 200
 201     for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
 202
 203     if(hPtJet[itype])     hPtJet[itype]->Write();
 204     if(hEtaPhiJet[itype]) hEtaPhiJet[itype]->Write();
 205     if(hRhoCent)          hRhoCent->Write();
 206     }
 207
 208     histOut->Close();
 209
 210   */
 211
 212   // Common output -
 213   // Added by sjena
 214
 215   TFile *fout = TFile::Open(outfile,"UPDATE");
 216   fout->ls();
 217
 218   TDirectoryFile *cdd = NULL;
 219   cdd = (TDirectoryFile*)fout->Get("JE");
 220   if(!cdd) {
 221     Printf("Warning: JE <dir> doesn't exist, creating a new one");
 222     cdd = (TDirectoryFile*)fout->mkdir("JE");
 223   }
 224   cdd->cd();
 225   cdd->ls();
 226
 227   for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
 228
 229     if(hPtJet[itype])     hPtJet[itype]->Write(Form("%s%d_%s",prefix.Data(), itype, hPtJet[itype]->GetName()));
 230     if(hEtaPhiJet[itype]) hEtaPhiJet[itype]->Write(Form("%s%d_%s",prefix.Data(), itype, hEtaPhiJet[itype]->GetName()));
 231     if(hRhoCent)          hRhoCent->Write(Form("%s%d_%s",prefix.Data(), itype, hRhoCent->GetName()));
 232   }
 233
 234   fout->Close();
 235
 236
 237
 238 }
 239 //__________________________________________________________
 240
 241 void SetHist(TH1* h,TString titx, TString tity){
 242
 243    h->GetXaxis()->SetTitle(titx.Data());
 244    h->GetYaxis()->SetTitle(tity.Data());
 245    h->GetXaxis()->SetTitleSize(0.06);
 246    h->GetYaxis()->SetTitleSize(0.06);
 247    h->GetYaxis()->SetTitleOffset(1.);
 248    h->GetXaxis()->SetTitleOffset(1.);
 249    h->SetLineWidth(3);
 250
 251 }
 252 //_____________________________________________________________________
 253
 254 void SetCanvas(TCanvas* c){
 255    c->SetLeftMargin(0.15);
 256    c->SetBottomMargin(0.15);
 257    c->SetRightMargin(0.05);
 258    c->SetTopMargin(0.05);
 259    c->SetTickx();
 260    c->SetTicky();
 261 }
 262 //_____________________________________________________________________
 263
 264 void SetLeg(TLegend* le){
 265    le->SetFillColor(10);
 266    le->SetBorderSize(0);
 267    le->SetFillStyle(0);
 268    le->SetTextSize(0.05);
 269 }