PWGPP/analysisQA/processJETriggerQA_V2.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_V2(TString strFileIn    = "AnalysisResults.root",
  12                            TString strFileIn2   = "",
  13                            TString suftype      = "eps",
  14                            Float_t jetR         = 0.2,
  15                            Float_t minTrkPT     = 0.15,
  16                            Float_t minClusterET = 0.3,
  17                            Int_t run            = 0,
  18                            TString trigsuffix   = "",
  19                            const char* outfile  = "JETriggerQA_outfile.root"
  20                            ){
  21
  22   gStyle->SetOptStat(0);
  23   gStyle->SetOptTitle(0);
  24
  25   TString prefix = "fig_je_TriggerQA_";
  26
  27   TFile * f1 = TFile::Open(strFileIn.Data());
  28   Bool_t drawComp = kFALSE;
  29   if(strFileIn2!="") {
  30     TFile * f2 = TFile::Open(strFileIn2.Data());
  31     drawComp = kTRUE;
  32   }
  33
  34   if(!trigsuffix=="" || !trigsuffix=="EJE" || !trigsuffix=="EGA")
  35   {cout << "Unknown trigger suffix. It should be either empty, EJE or EGA." << endl; return;}
  36
  37   //Load histogram list
  38   TString folder = 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%s/TriggerQA_Jet_AKTFullR%02d0_PicoTracks_pT%04d_CaloClustersCorr_ET%04d_pt_scheme_Jet_AKTChargedR%02d0_PicoTracks_pT%04d_CaloClustersCorr_ET%04d_pt_scheme_TC%s",
  39        TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),
  40        TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),trigsuffix.Data(),
  41        TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),
  42        TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),trigsuffix.Data());
  43   TList *histList[2];
  44   histList[0] =  (TList*)f1->Get(folder);
  45   if(histList[0]==0){
  46     cout << "Could not find " << folder << " in " << strFileIn << endl;
  47     return;
  48   }
  49   //Load histogram list for comparison
  50   if(drawComp) {
  51     TList *histList[1] = 0x0;
  52     histList[1] =  (TList*)f2->Get(folder);
  53     if(histList==0){
  54       cout << "Could not find " << folder << " in " << strFileIn2 << endl;
  55       return;
  56     }
  57   }
  58
  59   //---------------------------------------------------------------------------------------------------
  60   //       jet histograms
  61   //---------------------------------------------------------------------------------------------------
  62
  63   const Int_t kJetType = 2;
  64   TString suffix[kJetType] = {"Charged","Full"};
  65   TString suffix2[2] = {"","_old"};
  66   TH3F *h3PtEtaPhiJet[kJetType][2];
  67   TH1F *hPtJet[kJetType][2];
  68   TH2F *hEtaPhiJet[kJetType][2];
  69   TH2F *hRhoCent[kJetType][2];
  70   //Number [0] is for the original data, [1] is for the comparison graphs.
  71
  72   TH1F *hNEventSel[2]; Float_t nEvents[2];
  73   hNEventSel[0] = (TH1F*) histList[0]->FindObject("fhNEvents");
  74   Float_t nEvents[0]  = hNEventSel[0]->GetBinContent(2);
  75   if(drawComp) {
  76     TH1F *hNEventSel[1] = (TH1F*) histList[1]->FindObject("fhNEvents");
  77     Float_t nEvents[1]  = hNEventSel[1]->GetBinContent(2);
  78   }
  79
  80   TProfile *hTriggerbit = histList[0]->FindObject("fhTriggerbit");
  81   UInt_t triggerbit = hTriggerbit->GetBinContent(1);
  82
  83   for(Int_t itype = 0; itype < kJetType; itype++){
  84     for(Int_t comp = 0; comp < 2; comp++){
  85       if((comp==1)&&(drawComp==kFALSE)) continue;
  86       h3PtEtaPhiJet[itype][comp] = (TH3F*)  histList[comp]->FindObject(Form("fh3PtEtaPhiJet%s",suffix[itype].Data()));
  87       if(! h3PtEtaPhiJet[itype][comp]) continue;
  88
  89       //jet pt spectra
  90       Int_t binMin = 1;
  91
  92       if(ptmin>0.) binMin = h3PtEtaPhiJet[itype][comp]->GetXaxis()->FindBin(ptmin+0.00001);
  93       h3PtEtaPhiJet[itype][comp]->GetXaxis()->SetRange(binMin, h3PtEtaPhiJet[itype][comp]->GetNbinsX());
  94
  95       hPtJet[itype][comp] = (TH1F*) h3PtEtaPhiJet[itype][comp]->Project3D("x");
  96       hPtJet[itype][comp]->Scale(1./nEvents[comp],"width");
  97       SetHist((TH1F*) hPtJet[itype][comp],"p_{T,corr}^{jet} [GeV/c]","1/N_{evt} dN/dp_{T,corr}^{jet} [(GeV/c)^{-1}]");
  98       hPtJet[itype][comp]->SetName(Form("hPtJet%s%s",suffix[itype].Data(),suffix2[comp].Data()));
  99
 100       //eta versus phi
 101       hEtaPhiJet[itype][comp] = (TH2F*) h3PtEtaPhiJet[itype][comp]->Project3D("yz");
 102       SetHist((TH1F*) hEtaPhiJet[itype][comp],"#varphi^{jet} (rad)","#eta^{jet}");
 103       hEtaPhiJet[itype][comp]->SetName(Form("hEtaPhiJet%s%s",suffix[itype].Data(),suffix2[comp].Data()));
 104     }
 105   }
 106
 107   for(Int_t itype = 0; itype < kJetType; itype++){
 108     for(Int_t comp = 0; comp < 2; comp++){
 109       if((comp==1)&&(drawComp==kFALSE)) continue;
 110       //rho versus centrality
 111       hRhoCent[itype][comp] = (TH2F*) histList[comp]->FindObject(Form("fHistRhovsCent%s",suffix[itype].Data()));
 112       if(!hRhoCent[itype][comp]) continue;
 113       SetHist((TH1F*) hRhoCent[itype][comp], "Centrality (%)", Form("%s#rho [GeV/c*rad^{-1}]", (suffix[itype]=="Full")?"s":"" ));
 114       hRhoCent[itype][comp]->SetName(Form("hRhoCent%s%s",suffix[itype].Data(),suffix2[comp].Data()));
 115     }
 116   }
 117
 118   //______________
 119   //Draw histograms
 120
 121   TCanvas *c[100];
 122   TH1F *frame[100];
 123   Int_t nCan = 0;
 124   TLegend *leg;
 125
 126   TString longtrigname = "";
 127   UInt_t emcalPhysSelbit = 1<<31;
 128   if(triggerbit==AliVEvent::kAny) longtrigname = "kAny";
 129   else if(triggerbit==AliVEvent::kAnyINT) longtrigname = "kAnyINT";
 130   else if(triggerbit==emcalPhysSelbit) longtrigname = "EmcalPhysicsSelectionTask";
 131   else {
 132     if(triggerbit & AliVEvent::kCentral) longtrigname += "kCentral";
 133     if(triggerbit & AliVEvent::kSemiCentral) longtrigname += "kSemiCentral";
 134     if(triggerbit & AliVEvent::kMB)  longtrigname += "kMB";
 135   }
 136
 137   if(trigsuffix=="EJE")  longtrigname += " and kEMCEJE";
 138   if(trigsuffix=="EGA")  longtrigname += " and kEMCEGA";
 139
 140   for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
 141
 142      //draw pt spectrum
 143      if(!hPtJet[itype][0]) continue;
 144      c[nCan] = new TCanvas(Form("c%d",nCan),Form("c%d: Pt %s jets",nCan,suffix[itype].Data()),600,450);
 145      SetCanvas((TCanvas*) c[nCan]);
 146      c[nCan]->SetLogy();
 147
 148      frame[nCan] = gPad->DrawFrame(hPtJet[itype][0]->GetBinLowEdge(1),
 149                                    1e-7,
 150                                    hPtJet[itype][0]->GetBinLowEdge(hPtJet[itype][0]->GetNbinsX()+1),
 151                                    hPtJet[itype][0]->GetBinContent(hPtJet[itype][0]->GetMaximumBin())*2.);
 152
 153      SetHist((TH1F*) frame[nCan],hPtJet[itype][0]->GetXaxis()->GetTitle(),hPtJet[itype][0]->GetYaxis()->GetTitle());
 154
 155      if(drawComp) {
 156        hPtJet[itype][1]->SetLineColor(2);
 157        hPtJet[itype][1]->SetLineWidth(5);
 158        hPtJet[itype][1]->DrawCopy("same");
 159      }
 160      hPtJet[itype][0]->DrawCopy("same");
 161
 162      leg = new TLegend(0.15,0.5,0.88,0.88);
 163      SetLeg(leg);
 164
 165      TString txt =  Form("%s jets AKT R=%.1f",suffix[itype].Data(),jetR);
 166      if(run>0) txt += Form(" run:%d",run);
 167      leg->AddEntry((TObject*) 0, txt.Data(),"");
 168      leg->AddEntry((TObject*) 0, Form("Trigger: %s",longtrigname.Data()),"");
 169      leg->AddEntry((TObject*) 0, Form("p_{T,trk}> %d MeV/c",TMath::Nint(minTrkPT*1000)),"");
 170      if(itype==1) leg->AddEntry((TObject*) 0, Form("E_{T}>%d MeV/c",TMath::Nint(minClusterET*1000)),"");
 171      if(drawComp){
 172        leg->AddEntry(hPtJet[itype][0], "Recent data", "l");
 173        leg->AddEntry(hPtJet[itype][1], "Previous data", "l");
 174      }
 175      leg->AddEntry((TObject*) 0, Form("#it{N}_{events} = %.0f",nEvents),"");
 176      leg->Draw();
 177
 178      c[nCan]->SaveAs(Form("%s_Pt_AKT%02d_pT%04d_ET%04d_Run%d_Trigger%s_%s.%s",prefix.Data(),TMath::Nint(jetR*10),
 179                           TMath::Nint(minTrkPT*1000),TMath::Nint(minClusterET*1000), run, longtrigname.Data(), suffix[itype].Data(), suftype.Data()));
 180
 181      nCan++;
 182   }//end of loop over charged and full jets
 183   //_________________
 184   //draw eta versus phi
 185   for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
 186     if(!hEtaPhiJet[itype][0]) continue;
 187     c[nCan] = new TCanvas(Form("c%d",nCan),Form("c%d: eta-phi %s jets",nCan,suffix[itype].Data()),600,450);
 188     if(drawComp) c[nCan]->Divide(2,1);
 189
 190     for(Int_t comp = 0; comp < 2; comp++) {
 191       if((comp==1)&&(!drawComp)) continue;
 192       c[nCan]->cd(comp+1);
 193       SetCanvas(gPad);
 194       gPad->SetRightMargin(0.14);
 195 //      c[nCan]->SetRightMargin(0.15);
 196
 197       frame[nCan] = gPad->DrawFrame(hEtaPhiJet[itype][comp]->GetXaxis()->GetBinLowEdge(1),
 198                                      hEtaPhiJet[itype][comp]->GetYaxis()->GetBinLowEdge(1),
 199                                      hEtaPhiJet[itype][comp]->GetXaxis()->GetBinLowEdge(hEtaPhiJet[itype][comp]->GetNbinsX()),
 200                                      hEtaPhiJet[itype][comp]->GetYaxis()->GetBinLowEdge(hEtaPhiJet[itype][comp]->GetNbinsY()));
 201
 202       SetHist((TH1F*) frame[nCan],hEtaPhiJet[itype][comp]->GetXaxis()->GetTitle(),hEtaPhiJet[itype][comp]->GetYaxis()->GetTitle());
 203
 204       hEtaPhiJet[itype][comp]->DrawCopy("colz");
 205
 206       leg = new TLegend(0.15,0.5,0.88,0.88);
 207       SetLeg(leg);
 208
 209       TString txt =  Form("%s jets AKT R=%.1f",suffix[itype].Data(),jetR);
 210       if(run>0) txt += Form(" run:%d",run);
 211       leg->AddEntry((TObject*) 0, txt.Data(),"");
 212       if(drawComp) leg->AddEntry((TObject*) 0, comp ? "Recent production" : "Previous production", "");
 213       leg->AddEntry((TObject*) 0, Form("Trigger: %s",longtrigname.Data()),"");
 214       leg->AddEntry((TObject*) 0, Form("p_{T,trk}> %d MeV/c",TMath::Nint(minTrkPT*1000)),"");
 215       if(itype==1) leg->AddEntry((TObject*) 0, Form("E_{T}>%d MeV/c",TMath::Nint(minClusterET*1000)),"");
 216       leg->AddEntry((TObject*) 0, Form("#it{N}_{events} = %.0f",nEvents),"");
 217       leg->AddEntry((TObject*) 0, Form("p_{T,corr}^{jet} > %.1f GeV/c", ptmin),"");
 218       leg->Draw();
 219     }
 220
 221     c[nCan]->SaveAs(Form("%s_EtaPhi_AKT%02d_pT%04d_ET%04d_Run%d_Trigger%s_%s.%s",prefix.Data(),TMath::Nint(jetR*10),
 222                           TMath::Nint(minTrkPT*1000),TMath::Nint(minClusterET*1000), run, longtrigname.Data(), suffix[itype].Data(), suftype.Data()));
 223
 224     nCan++;
 225   }//end of the loop over charged and full jets
 226   //_________________
 227   //draw rho versus centrality
 228   for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
 229      if(!hRhoCent[itype][0]) continue;
 230      c[nCan] = new TCanvas(Form("c%d",nCan),Form("c%d: Rho-Cent",nCan),600,450);
 231      if(drawComp) c[nCan]->Divide(2,1);
 232
 233      for(Int_t comp = 0; comp < 2; comp++) {
 234        if((comp==1)&&(!drawComp)) continue;
 235        c[nCan]->cd(comp+1);
 236        SetCanvas(gPad);
 237        gPad->SetRightMargin(0.15);
 238        c[nCan]->SetLogz();
 239
 240        frame[nCan] = gPad->DrawFrame(hRhoCent[itype][comp]->GetXaxis()->GetBinLowEdge(1),
 241                                      hRhoCent[itype][comp]->GetYaxis()->GetBinLowEdge(1),
 242                                      hRhoCent[itype][comp]->GetXaxis()->GetBinLowEdge(hRhoCent[itype][comp]->GetNbinsX()),
 243                                      hRhoCent[itype][comp]->GetYaxis()->GetBinLowEdge(hRhoCent[itype][comp]->GetNbinsY()));
 244
 245        SetHist((TH1F*) frame[nCan],hRhoCent[itype][comp]->GetXaxis()->GetTitle(),hRhoCent[itype][comp]->GetYaxis()->GetTitle());
 246
 247        hRhoCent[itype][comp]->DrawCopy("colz");
 248
 249        leg = new TLegend(-0.05,0.5,0.88,0.88);
 250        SetLeg(leg);
 251
 252        TString txt =  Form("%s jets AKT R=%.1f",suffix[itype].Data(),jetR);
 253        if(run>0) txt += Form(" run:%d",run);
 254
 255        leg->AddEntry((TObject*) 0, txt.Data(),"");
 256        if(drawComp) leg->AddEntry((TObject*) 0, comp ? "Recent production" : "Previous production", "");
 257        leg->AddEntry((TObject*) 0, Form("Trigger: %s",longtrigname.Data()),"");
 258        leg->AddEntry((TObject*) 0, Form("#it{N}_{events} = %.0f",nEvents),"");
 259        leg->Draw();
 260     }
 261
 262     c[nCan]->SaveAs(Form("%s_RhoCent_AKT%02d_pT%04d_ET%04d_Run%d_Trigger%s_%s.%s",prefix.Data(),TMath::Nint(jetR*10),
 263                          TMath::Nint(minTrkPT*1000),TMath::Nint(minClusterET*1000), run, longtrigname.Data(), suffix[itype].Data(), suftype.Data()));
 264     nCan++;
 265   }//end of the loop over charged and full jets
 266
 267   //---------------------------------------------------------------------------------------------------
 268   //                       WRITE OUTPUT TO ROOT FILE
 269   //---------------------------------------------------------------------------------------------------
 270
 271
 272   /* Standalone output
 273     TFile *histOut = new TFile(Form("%s_AKT%02d_pT%04d_ET%04d_jetPtMin%.1f_Run%d.root",prefix.Data(),
 274     TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),ptmin,run),"RECREATE");
 275
 276     for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
 277
 278     if(hPtJet[itype])     hPtJet[itype]->Write();
 279     if(hEtaPhiJet[itype]) hEtaPhiJet[itype]->Write();
 280     if(hRhoCent[itype])   hRhoCent[itype]->Write();
 281     }
 282
 283     histOut->Close();
 284
 285   */
 286
 287   // Common output -
 288   // Added by sjena
 289
 290   TFile *fout = TFile::Open(outfile,"UPDATE");
 291   fout->ls();
 292
 293   TDirectoryFile *cdd = NULL;
 294   cdd = (TDirectoryFile*)fout->Get("JE");
 295   if(!cdd) {
 296     Printf("Warning: JE <dir> doesn't exist, creating a new one");
 297     cdd = (TDirectoryFile*)fout->mkdir("JE");
 298   }
 299   cdd->cd();
 300   cdd->ls();
 301
 302   for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
 303
 304     if(hPtJet[itype][0])     hPtJet[itype][0]->Write(Form("%s%d_%s",prefix.Data(), itype, hPtJet[itype][0]->GetName()));
 305     if(drawComp) if(hPtJet[itype][1])     hPtJet[itype][1]->Write(Form("%s%d_%s",prefix.Data(), itype, hPtJet[itype][1]->GetName()));
 306     if(hEtaPhiJet[itype][0]) hEtaPhiJet[itype][0]->Write(Form("%s%d_%s",prefix.Data(), itype, hEtaPhiJet[itype][0]->GetName()));
 307     if(drawComp) if(hEtaPhiJet[itype][1]) hEtaPhiJet[itype][1]->Write(Form("%s%d_%s",prefix.Data(), itype, hEtaPhiJet[itype][1]->GetName()));
 308     if(hRhoCent[itype][0])   hRhoCent[itype][0]->Write(Form("%s%d_%s",prefix.Data(), itype, hRhoCent[itype][0]->GetName()));
 309     if(drawComp) if(hRhoCent[itype][1])   hRhoCent[itype][1]->Write(Form("%s%d_%s",prefix.Data(), itype, hRhoCent[itype][1]->GetName()));
 310   }
 311
 312   fout->Close();
 313
 314
 315
 316 }
 317 //__________________________________________________________
 318
 319 void SetHist(TH1* h,TString titx, TString tity){
 320
 321    h->GetXaxis()->SetTitle(titx.Data());
 322    h->GetYaxis()->SetTitle(tity.Data());
 323    h->GetXaxis()->SetTitleSize(0.06);
 324    h->GetYaxis()->SetTitleSize(0.06);
 325    h->GetYaxis()->SetTitleOffset(1.);
 326    h->GetXaxis()->SetTitleOffset(1.);
 327    h->SetLineWidth(3);
 328
 329 }
 330 //_____________________________________________________________________
 331
 332 void SetCanvas(TVirtualPad* c){
 333    c->SetLeftMargin(0.15);
 334    c->SetBottomMargin(0.15);
 335    c->SetRightMargin(0.05);
 336    c->SetTopMargin(0.05);
 337    c->SetTickx();
 338    c->SetTicky();
 339 }
 340 //_____________________________________________________________________
 341
 342 void SetLeg(TLegend* le){
 343    le->SetFillColor(10);
 344    le->SetBorderSize(0);
 345    le->SetFillStyle(0);
 346    le->SetTextSize(0.05);
 347 }