[u/mrichter/AliRoot.git] / PWGPP / analysisQA / processJETriggerQA_V2.C

/***************************************************
processJETriggerQA:
To procees JE Triggere QA wagon's output

****************************************************/


const Float_t ptmin =  0. ; //lower cutoff of jet pt spectrum

void processJETriggerQA_V2(TString strFileIn    = "AnalysisResults.root", 
			   TString suftype      = "eps",
			   Float_t jetR         = 0.2, 
			   Float_t minTrkPT     = 0.15, 
			   Float_t minClusterET = 0.3, 
			   Int_t run            = 0,
			   TString trigsuffix   = "",
			   const char* outfile  = "JETriggerQA_outfile.root"
			   ){

  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);

  TString prefix = "fig_je_TriggerQA_";

  TFile * f1 = TFile::Open(strFileIn.Data());

  if(!trigsuffix=="" || !trigsuffix=="EJE" || !trigsuffix=="EGA")
  {cout << "Unknown trigger suffix. It should be either empty, EJE or EGA." << endl; return;}

  //Load histogram list
  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", 
       TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),
       TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),trigsuffix.Data(),
       TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),  
       TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),trigsuffix.Data());
  TList *histList = 0x0;
  histList =  (TList*)f1->Get(folder);
  if(histList==0){
    cout << "Could not find " << folder << " in " << strFileIn << endl;
    return;
  }

  //---------------------------------------------------------------------------------------------------
  //       jet histograms
  //---------------------------------------------------------------------------------------------------

  const Int_t kJetType = 2;
  TString suffix[kJetType] = {"Charged","Full"};
  TH3F *h3PtEtaPhiJet[kJetType];
  TH1F *hPtJet[kJetType];
  TH2F *hEtaPhiJet[kJetType];
  TH2F *hRhoCent[kJetType];

  TH1F *hNEventSel = histList->FindObject("fhNEvents");
  Float_t nEvents  = hNEventSel->GetBinContent(2);

  TProfile *hTriggerbit = histList->FindObject("fhTriggerbit");
  UInt_t triggerbit = hTriggerbit->GetBinContent(1);

  for(Int_t itype = 0; itype < kJetType; itype++){

     h3PtEtaPhiJet[itype] = (TH3F*)  histList->FindObject(Form("fh3PtEtaPhiJet%s",suffix[itype].Data()));
     if(! h3PtEtaPhiJet[itype]) continue;

     //jet pt spectra
     Int_t binMin = 1;
     if(ptmin>0.) binMin = h3PtEtaPhiJet[itype]->GetXaxis()->FindBin(ptmin+0.00001);
     h3PtEtaPhiJet[itype]->GetXaxis()->SetRange(binMin, h3PtEtaPhiJet[itype]->GetNbinsX());

     hPtJet[itype] = (TH1F*) h3PtEtaPhiJet[itype]->Project3D("x");
     if(nEvents>0)
     hPtJet[itype]->Scale(1./nEvents,"width");
     SetHist((TH1F*) hPtJet[itype],"p_{T,corr}^{jet} (GeV)","1/N_{evt} dN/dp_{T,corr}^{jet} (GeV^{-1})");
     hPtJet[itype]->SetName(Form("hPtJet%s",suffix[itype].Data())); 

     //eta versus phi
     hEtaPhiJet[itype] = (TH2F*) h3PtEtaPhiJet[itype]->Project3D("yz");
     SetHist((TH1F*) hEtaPhiJet[itype],"#varphi^{jet} (rad)","#eta^{jet}");
     hEtaPhiJet[itype]->SetName(Form("hEtaPhiJet%s",suffix[itype].Data()));
  }

  for(Int_t itype = 0; itype < kJetType; itype++){
  //rho versus centrality
    hRhoCent[itype] = (TH2F*) histList->FindObject(Form("fHistRhovsCent%s",suffix[itype].Data()));
     if(!hRhoCent[itype]) continue;
     SetHist((TH1F*) hRhoCent[itype], "Centrality (%)", Form("%s#rho (GeV/c*rad^{-1})", (suffix[itype]=="Full")?"s":"" ));
     hRhoCent[itype]->SetName(Form("hRhoCent%s",suffix[itype].Data()));
  }

  //______________
  //Draw histograms

  TCanvas *c[100];
  TH1F *frame[100];
  Int_t nCan = 0;
  TLegend *leg; 

  TString longtrigname = "";
  UInt_t emcalPhysSelbit = 1<<31;
  if(triggerbit==AliVEvent::kAny) longtrigname = "kAny";
  else if(triggerbit==AliVEvent::kAnyINT) longtrigname = "kAnyINT";
  else if(triggerbit==emcalPhysSelbit) longtrigname = "EmcalPhysicsSelectionTask";
  else {
    if(triggerbit & AliVEvent::kCentral) longtrigname += "kCentral";
    if(triggerbit & AliVEvent::kSemiCentral) longtrigname += "kSemiCentral";
    if(triggerbit & AliVEvent::kMB)  longtrigname += "kMB";
  }

  if(trigsuffix=="EJE")  longtrigname += "AndkEMCEJE";
  if(trigsuffix=="EGA")  longtrigname += "AndkEMCEGA";
 
  for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets

     //draw pt spectrum
     if(!hPtJet[itype]) continue; 
     c[nCan] = new TCanvas(Form("c%d",nCan),Form("c%d: Pt %s jets",nCan,suffix[itype].Data()),600,450);
     SetCanvas((TCanvas*) c[nCan]);
     c[nCan]->SetLogy();

     frame[nCan] = gPad->DrawFrame(hPtJet[itype]->GetBinLowEdge(1),  
                                   1e-7,
                                   hPtJet[itype]->GetBinLowEdge(hPtJet[itype]->GetNbinsX()+1), 
                                   hPtJet[itype]->GetBinContent(hPtJet[itype]->GetMaximumBin())*2.);

     SetHist((TH1F*) frame[nCan],hPtJet[itype]->GetXaxis()->GetTitle(),hPtJet[itype]->GetYaxis()->GetTitle());
 

     hPtJet[itype]->DrawCopy("same");
 
     leg = new TLegend(0.15,0.5,0.88,0.88);
     SetLeg(leg);
     
     TString txt =  Form("%s jets AKT R=%.1f",suffix[itype].Data(),jetR);
     if(run>0) txt += Form(" run:%d",run); 
     leg->AddEntry((TObject*) 0, txt.Data(),"");
     leg->AddEntry((TObject*) 0, Form("Trigger: %s",longtrigname.Data()),"");
     leg->AddEntry((TObject*) 0, Form("p_{T,trk}> %d MeV",TMath::Nint(minTrkPT*1000)),"");
     if(itype==1) leg->AddEntry((TObject*) 0, Form("E_{T}>%d MeV",TMath::Nint(minClusterET*1000)),""); 
     leg->AddEntry((TObject*) 0, Form("#it{N}_{events} = %.0f",nEvents),"");
     leg->Draw();

     c[nCan]->SaveAs(Form("%s_Pt_AKT%02d_pT%04d_ET%04d_Run%d_Trigger%s_%s.%s",prefix.Data(),TMath::Nint(jetR*10),
			  TMath::Nint(minTrkPT*1000),TMath::Nint(minClusterET*1000), run, longtrigname.Data(), suffix[itype].Data(), suftype.Data()));

     nCan++;
  }//end of loop over charged and full jets
  //_________________
  //draw eta versus phi
  for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
     if(!hEtaPhiJet[itype]) continue; 
     c[nCan] = new TCanvas(Form("c%d",nCan),Form("c%d: eta-phi %s jets",nCan,suffix[itype].Data()),600,450);
     SetCanvas((TCanvas*) c[nCan]);
     c[nCan]->SetRightMargin(0.15);

     frame[nCan] = gPad->DrawFrame(hEtaPhiJet[itype]->GetXaxis()->GetBinLowEdge(1),  
                                   hEtaPhiJet[itype]->GetYaxis()->GetBinLowEdge(1),
                                   hEtaPhiJet[itype]->GetXaxis()->GetBinLowEdge(hEtaPhiJet[itype]->GetNbinsX()),
                                   hEtaPhiJet[itype]->GetYaxis()->GetBinLowEdge(hEtaPhiJet[itype]->GetNbinsY()));

     SetHist((TH1F*) frame[nCan],hEtaPhiJet[itype]->GetXaxis()->GetTitle(),hEtaPhiJet[itype]->GetYaxis()->GetTitle());
 

     hEtaPhiJet[itype]->DrawCopy("same,colz");
 
     leg = new TLegend(0.15,0.5,0.88,0.88);
     SetLeg(leg);

     TString txt =  Form("%s jets AKT R=%.1f",suffix[itype].Data(),jetR);
     if(run>0) txt += Form(" run:%d",run); 
     leg->AddEntry((TObject*) 0, txt.Data(),"");
     leg->AddEntry((TObject*) 0, Form("Trigger: %s",longtrigname.Data()),"");
     leg->AddEntry((TObject*) 0, Form("p_{T,trk}> %d MeV",TMath::Nint(minTrkPT*1000)),"");
     if(itype==1) leg->AddEntry((TObject*) 0, Form("E_{T}>%d MeV",TMath::Nint(minClusterET*1000)),"");
     leg->AddEntry((TObject*) 0, Form("#it{N}_{events} = %.0f",nEvents),"");
     leg->AddEntry((TObject*) 0, Form("p_{T,corr}^{jet} > %.1f GeV", ptmin),"");
     leg->Draw();

     c[nCan]->SaveAs(Form("%s_EtaPhi_AKT%02d_pT%04d_ET%04d_Run%d_Trigger%s_%s.%s",prefix.Data(),TMath::Nint(jetR*10),
			  TMath::Nint(minTrkPT*1000),TMath::Nint(minClusterET*1000), run, longtrigname.Data(), suffix[itype].Data(), suftype.Data()));

     nCan++;
  }//end of the loop over charged and full jets
  //_________________
  //draw rho versus centrality
  for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
     if(hRhoCent[itype]) {
       c[nCan] = new TCanvas(Form("c%d",nCan),Form("c%d: Rho-Cent",nCan),600,450);
       SetCanvas((TCanvas*) c[nCan]);
       c[nCan]->SetRightMargin(0.15);
       c[nCan]->SetLogz();
  
       frame[nCan] = gPad->DrawFrame(hRhoCent[itype]->GetXaxis()->GetBinLowEdge(1),
                                     hRhoCent[itype]->GetYaxis()->GetBinLowEdge(1),
                                     hRhoCent[itype]->GetXaxis()->GetBinLowEdge(hRhoCent[itype]->GetNbinsX()),
                                     hRhoCent[itype]->GetYaxis()->GetBinLowEdge(hRhoCent[itype]->GetNbinsY()));
  
       SetHist((TH1F*) frame[nCan],hRhoCent[itype]->GetXaxis()->GetTitle(),hRhoCent[itype]->GetYaxis()->GetTitle());
  
       hRhoCent[itype]->DrawCopy("colz");
  
       leg = new TLegend(0.15,0.5,0.88,0.88);
       SetLeg(leg);
  
       TString txt =  Form("%s jets AKT R=%.1f",suffix[itype].Data(),jetR);
       if(run>0) txt += Form(" run:%d",run); 
  
       leg->AddEntry((TObject*) 0, txt.Data(),"");
       leg->AddEntry((TObject*) 0, Form("Trigger: %s",longtrigname.Data()),"");
       leg->AddEntry((TObject*) 0, Form("#it{N}_{events} = %.0f",nEvents),"");
       leg->Draw();
  
       c[nCan]->SaveAs(Form("%s_RhoCent_AKT%02d_pT%04d_ET%04d_Run%d_Trigger%s_%s.%s",prefix.Data(),TMath::Nint(jetR*10),
			    TMath::Nint(minTrkPT*1000),TMath::Nint(minClusterET*1000), run, longtrigname.Data(), suffix[itype].Data(), suftype.Data()));
       nCan++; 
    }
  }//end of the loop over charged and full jets

  //---------------------------------------------------------------------------------------------------
  //                       WRITE OUTPUT TO ROOT FILE
  //---------------------------------------------------------------------------------------------------


  /* Standalone output
    TFile *histOut = new TFile(Form("%s_AKT%02d_pT%04d_ET%04d_jetPtMin%.1f_Run%d.root",prefix.Data(),
    TMath::Nint(jetR*10), TMath::Nint(minTrkPT*1000), TMath::Nint(minClusterET*1000),ptmin,run),"RECREATE");
    
    for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
     
    if(hPtJet[itype])     hPtJet[itype]->Write();
    if(hEtaPhiJet[itype]) hEtaPhiJet[itype]->Write();
    if(hRhoCent[itype])   hRhoCent[itype]->Write();
    }
    
    histOut->Close();
  
  */

  // Common output - 
  // Added by sjena

  TFile *fout = TFile::Open(outfile,"UPDATE");
  fout->ls();
  
  TDirectoryFile *cdd = NULL;
  cdd = (TDirectoryFile*)fout->Get("JE");
  if(!cdd) {
    Printf("Warning: JE <dir> doesn't exist, creating a new one");
    cdd = (TDirectoryFile*)fout->mkdir("JE");
  }
  cdd->cd();
  cdd->ls();
  
  for(Int_t itype = 0; itype < kJetType; itype++){ //loop over charged and full jets
    
    if(hPtJet[itype])     hPtJet[itype]->Write(Form("%s%d_%s",prefix.Data(), itype, hPtJet[itype]->GetName()));
    if(hEtaPhiJet[itype]) hEtaPhiJet[itype]->Write(Form("%s%d_%s",prefix.Data(), itype, hEtaPhiJet[itype]->GetName()));
    if(hRhoCent[itype])   hRhoCent[itype]->Write(Form("%s%d_%s",prefix.Data(), itype, hRhoCent[itype]->GetName()));
  }
  
  fout->Close();


}
//__________________________________________________________

void SetHist(TH1* h,TString titx, TString tity){

   h->GetXaxis()->SetTitle(titx.Data());
   h->GetYaxis()->SetTitle(tity.Data());
   h->GetXaxis()->SetTitleSize(0.06);
   h->GetYaxis()->SetTitleSize(0.06);
   h->GetYaxis()->SetTitleOffset(1.);
   h->GetXaxis()->SetTitleOffset(1.);
   h->SetLineWidth(3);

}
//_____________________________________________________________________

void SetCanvas(TCanvas* c){
   c->SetLeftMargin(0.15);
   c->SetBottomMargin(0.15);
   c->SetRightMargin(0.05);
   c->SetTopMargin(0.05);
   c->SetTickx();
   c->SetTicky();
}
//_____________________________________________________________________

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