/************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * * Author: The ALICE Off-line Project. * * Contributors are mentioned in the code where appropriate. * * * * Permission to use, copy, modify and distribute this software and its * * documentation strictly for non-commercial purposes is hereby granted * * without fee, provided that the above copyright notice appears in all * * copies and that both the copyright notice and this permission notice * * appear in the supporting documentation. The authors make no claims * * about the suitability of this software for any purpose. It is * * provided "as is" without express or implied warranty. * **************************************************************************/ //--------------------------------------------------------------------- // JetAnalysis class // Analyse Jets // Author: andreas.morsch@cern.ch //--------------------------------------------------------------------- #include "AliJetAnalysis.h" ClassImp(AliJetAnalysis) //////////////////////////////////////////////////////////////////////// #include #include #include #include #include #include #include #include #include "AliJetProductionDataPDC2004.h" #include "AliJet.h" #include "AliJetESDReaderHeader.h" #include "AliUA1JetHeader.h" #include "AliLeading.h" AliJetAnalysis::AliJetAnalysis() { // Constructor fDirectory = 0x0; fEventMin = 0; fEventMax = -1; fRunMin = 0; fRunMax = 11; } void AliJetAnalysis::Analyse() { // // Some histos // TH1F::AddDirectory(0); TProfile::AddDirectory(0); TH1F* e0H = new TH1F("e0H" ,"Jet Energy (reconstructed)", 40, 0., 200.); TH1F* e1H = new TH1F("e1H" , "Jet Energy (generated)", 40, 0., 200.); TH1F* e2H = new TH1F("e2H" , "Jet Energy (generated, nrec = 0", 40, 0., 200.); TH1F* e3H = new TH1F("e3H" , "Jet Energy (leading)", 40, 0., 200.); TH1F* e4H = new TH1F("e4H" , "Jet Energy (reconstructed: 105 < Egen < 125", 40, 0., 200.); TH1F* e5H = new TH1F("e5H" , "Jet Energy (generated)", 40, 0., 200.); TH1F* e6H = new TH1F("e6H" , "Jet Energy (generated)", 40, 0., 200.); TH1F* e7H = new TH1F("e7H" , "Jet Energy (generated)", 40, 0., 200.); TH1F* e8H = new TH1F("e8H" , "Jet Energy (generated)", 40, 0., 200.); TProfile* r5H = new TProfile("r5H" , "rec/generated", 20, 0., 200, 0., 1., "S"); TProfile* r6H = new TProfile("r6H" , "rec/generated", 20, 0., 200, 0., 1., "S"); TProfile* r7H = new TProfile("r7H" , "rec/generated", 20, 0., 200, 0., 1., "S"); TProfile* r8H = new TProfile("r8H" , "rec/generated", 20, 0., 200, 0., 1., "S"); TH1F* dr1H = new TH1F("dr1H", "delta R", 160, 0., 2.); TH1F* dr2H = new TH1F("dr2H", "delta R", 160, 0., 2.); TH1F* dr3H = new TH1F("dr4H", "delta R", 160, 0., 2.); TH1F* etaH = new TH1F("etaH", "eta", 160, -2., 2.); TH1F* eta1H = new TH1F("eta1H", "eta", 160, -2., 2.); TH1F* eta2H = new TH1F("eta2H", "eta", 160, -2., 2.); TH1F* phiH = new TH1F("phiH", "phi", 160, -3., 3.); TH1F* dphiH = new TH1F("dphiH", "phi", 160, 0., 3.1415); TH1F* phi1H = new TH1F("phi1H", "phi", 160, 0., 6.28); TH1F* phi2H = new TH1F("phi2H", "phi", 160, 0., 6.28); TProfile* drP1 = new TProfile("drP1" , "Delta_R", 20, 0., 200, -1., 1., "S"); TProfile* drP2 = new TProfile("drP2" , "Delta_R", 20, 0., 200, -1., 1., "S"); // Run data AliJetProductionDataPDC2004* runData = new AliJetProductionDataPDC2004(); // Loop over runs TFile* jFile = 0x0; for (Int_t iRun = fRunMin; iRun <= fRunMax; iRun++) { // Open file char fn[20]; sprintf(fn, "%s/%s.root", fDirectory, (runData->GetRunTitle(iRun)).Data()); jFile = new TFile(fn); printf(" Analyzing run: %d %s\n", iRun,fn); // Get jet header and display parameters AliUA1JetHeader* jHeader = (AliUA1JetHeader*) (jFile->Get("AliUA1JetHeader")); // jHeader->PrintParameters(); // Get reader header and events to be looped over AliJetESDReaderHeader *jReaderH = (AliJetESDReaderHeader*)(jFile->Get("AliJetKineReaderHeader")); if (fEventMin == -1) fEventMin = jReaderH->GetFirstEvent(); if (fEventMax == -1) { fEventMax = jReaderH->GetLastEvent(); } else { fEventMax = TMath::Min(fEventMax, jReaderH->GetLastEvent()); } // Calculate weight Float_t wgt = runData->GetWeight(iRun) / Float_t(fEventMax - fEventMin + 1); Float_t ptmin, ptmax; runData->GetPtHardLimits(iRun, ptmin, ptmax); // Loop over events AliJet *jets = 0x0; AliJet *gjets = 0x0; AliLeading *leading = 0x0; Float_t egen, etag, econ, erec; for (Int_t i = fEventMin; i < fEventMax; i++) { printf(" Analyzing run: %d Event %d / %d \n", iRun, i, fEventMax); // Het next tree with AliJet char nameT[100]; sprintf(nameT, "TreeJ%d",i); TTree *jetT =(TTree *)(jFile->Get(nameT)); jetT->SetBranchAddress("FoundJet", &jets); jetT->SetBranchAddress("GenJet", &gjets); jetT->SetBranchAddress("LeadingPart", &leading); jetT->GetEntry(0); // // Find the jet with the highest E_T within fiducial region // Int_t njets = jets->GetNJets(); Int_t imax = -1; Int_t jmax = -1; Float_t emax = 0.; for (Int_t ij = 0; ij < njets; ij++) { if (jets->GetPt(ij) > emax && TMath::Abs(jets->GetEta(ij)) < 0.60) { emax = jets->GetPt(ij); jmax = imax; imax = ij; } } if (imax == -1) { Int_t ngen = gjets->GetNJets(); if(ngen>0) e2H->Fill(gjets->GetPt(0), wgt); } else { // printf("Reconstructed Jet %5d %13.3f %13.3f %13.3f\n", // imax, emax, jets->GetEta(imax), jets->GetPhi(imax)); econ = jets->GetPt(imax); erec = jets->GetPt(imax) / 0.65; dr1H->Fill(jets->GetEta(imax) - gjets->GetEta(0)); // // Find the generated jet closest to the reconstructed // Float_t rmin; Float_t etamin = 1.e6; Int_t igen; Float_t etaj = jets->GetEta(imax); Float_t phij = jets->GetPhi(imax); Int_t ngen = gjets->GetNJets(); if (ngen != 0) { rmin = 1.e6; igen = -1; for (Int_t j = 0; j < ngen; j++) { etag = gjets->GetEta(j); Float_t phig = gjets->GetPhi(j); Float_t deta = etag - etaj; Float_t dphi = TMath::Abs(phig - phij); if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi; Float_t r = TMath::Sqrt(deta * deta + dphi * dphi); if (r < rmin) { rmin = r; etamin = deta; igen = j; } } egen = gjets->GetPt(igen); e1H->Fill(gjets->GetPt(igen), wgt); etag = gjets->GetEta(igen); Float_t phig = gjets->GetPhi(igen); Float_t dphi = phig - phij; // if (econ < ptmax) { e0H->Fill(erec, wgt); // } else { // e0H->Fill(erec, 6.7e-6); // } if (egen > 20. && egen < 40.) { phiH->Fill((dphi)); etaH->Fill(etag - etaj, wgt); phi1H->Fill(phig); phi2H->Fill(phij); e4H->Fill(jets->GetPt(imax)); } if (erec > 90. && erec < 110. && rmin < 0.1) { e5H->Fill(egen, wgt); dr2H->Fill(rmin); if (egen < 30.) { printf("Strange jet %6d %13.3f %13.3f %13.3f \n", imax, etaj, phij, erec); for (Int_t j = 0; j < ngen; j++) { printf("Generated %6d %13.3f %13.3f %13.3f\n", j, gjets->GetEta(j), gjets->GetPhi(j), gjets->GetPt(j)); } } } if (rmin < 0.1) { r5H->Fill(egen, jets->GetPt(imax)/egen, wgt); r6H->Fill(jets->GetPt(imax) / 0.4, jets->GetPt(imax)/egen, wgt); e8H->Fill(erec, wgt); } if (rmin < 0.1) { drP1->Fill(egen, etamin, wgt); } } // ngen !=0 } // has reconstructed jet // // Leading particle // if (leading->LeadingFound()) { Float_t etal = leading->GetLeading()->Eta(); Float_t phil = leading->GetLeading()->Phi(); Float_t el = leading->GetLeading()->E(); // printf("Leading %f %f %f \n", etal, phil, el); e3H->Fill(el, wgt); Float_t rmin = 1.e6; Float_t etamin = 1.e6; Int_t igen = -1; Int_t ngen = gjets->GetNJets(); for (Int_t j = 0; j < ngen; j++) { etag = gjets->GetEta(j); Float_t phig = gjets->GetPhi(j); Float_t deta = etag-etal; Float_t dphi = TMath::Abs(phig - phil); if (dphi > TMath::Pi()) dphi = 2. * TMath::Pi() - dphi; Float_t r = TMath::Sqrt(deta * deta + dphi * dphi); if (r < rmin) { rmin = r; igen = j; etamin = deta; } } if (egen > 20. && egen < 40.) { dr3H->Fill(rmin); eta1H->Fill(etag-etal, wgt); } if (el > 54. && el < 66.) e6H->Fill(egen, wgt); if (rmin < 0.3) { r7H->Fill(egen, el/egen, wgt); r8H->Fill(el / 0.2, el/egen, wgt); } if (rmin < 0.1) { drP2->Fill(egen, etamin, wgt); } } // if leading particle // // Generated Jet // Int_t ngen = gjets->GetNJets(); emax = 0.; imax = -1; for (Int_t j = 0; j < ngen; j++) { if (gjets->GetPt(j) > emax && TMath::Abs(gjets->GetEta(j)) < 0.5) { emax = gjets->GetPt(j); imax = j; } } if (imax != -1) e7H->Fill(emax, wgt); delete jetT; } // events if (jFile) jFile->Close(); delete jFile; } // runs // delete jFile; // if (jFile) jFile->Close(); /* TFile* f = new TFile("j.root", "recreate"); e0H->Write(); e1H->Write(); e2H->Write(); e3H->Write(); e4H->Write(); e7H->Write(); e8H->Write(); f->Close(); */ // Get Control Plots // gStyle->SetOptStat(0); TCanvas* c1 = new TCanvas("c1"); e0H->Draw(); e1H->SetLineColor(2); e2H->SetLineColor(4); e3H->SetLineColor(5); e1H->Draw("same"); e3H->Draw("same"); TCanvas* c2 = new TCanvas("c2"); // dr1H->Draw(); dr2H->SetLineColor(2); dr2H->Draw(""); TCanvas* c3 = new TCanvas("c3"); dr2H->Draw(); dr3H->Draw("same"); TCanvas* c4 = new TCanvas("c4"); e0H->Draw(); TCanvas* c5 = new TCanvas("c5"); etaH->SetXTitle("#eta_{rec} - #eta_{gen}"); etaH->Draw(); eta1H->SetLineColor(2); eta1H->Draw("same"); TCanvas* c5a = new TCanvas("c5a"); eta1H->Draw(); TCanvas* c5b = new TCanvas("c5b"); eta2H->Draw(); TCanvas* c6 = new TCanvas("c6"); e4H->Draw(); TCanvas* c7 = new TCanvas("c7"); phiH->Draw(); TCanvas* c7a = new TCanvas("c7a"); phi1H->Draw(); TCanvas* c7b = new TCanvas("c7b"); phi2H->Draw(); TCanvas* c8 = new TCanvas("c8"); e5H->SetXTitle("E_{gen} (GeV)"); e5H->Draw(); e6H->SetLineColor(2); e6H->Draw("same"); TCanvas* c9 = new TCanvas("c9"); e6H->Draw(); TCanvas* c10 = new TCanvas("c10"); r5H->SetMaximum(1.); r5H->Draw(); r5H->SetXTitle("E_{gen} (GeV)"); r5H->SetYTitle("E_{leading} / E_{gen}"); r6H->SetLineColor(2); r6H->Draw("same"); TCanvas* c11 = new TCanvas("c11"); // // Leading particle rec/gen // r7H->SetMaximum(1.); r7H->Draw(); r7H->SetXTitle("E_{gen} (GeV)"); r7H->SetYTitle("E_{leading} / E_{gen}"); r8H->SetLineColor(2); r8H->Draw("same"); TCanvas* c12 = new TCanvas("c12"); drP1->SetXTitle("E_{gen} (GeV)"); drP1->SetYTitle("#eta_{rec} - #eta_{gen}"); drP1->Draw(); TCanvas* c13 = new TCanvas("c13"); drP2->SetXTitle("E_{gen} (GeV)"); drP2->SetYTitle("#eta_{leading} - #eta_{gen}"); drP2->Draw(); TCanvas* c14 = new TCanvas("c14"); dphiH->Draw(); /* e1H->Write(); e2H->Write(); e3H->Write(); e4H->Write(); */ }