ClassImp(AliAnalysisTaskJetSpectrum)
+const Float_t AliAnalysisTaskJetSpectrum::fgkJetNpartCut[AliAnalysisTaskJetSpectrum::kMaxCorrelation] = {5,10,1E+09};
+
AliAnalysisTaskJetSpectrum::AliAnalysisTaskJetSpectrum(): AliAnalysisTaskSE(),
fJetHeaderRec(0x0),
fJetHeaderGen(0x0),
fh1PtHard_Trials(0x0),
fh1NGenJets(0x0),
fh1NRecJets(0x0),
- fHistList(0x0)
+ fHistList(0x0) ,
+ ////////////////
+ fh1JetMultiplicity(0x0) ,
+ fh2ERecZRec(0x0) ,
+ fh2EGenZGen(0x0) ,
+ fh2Efficiency(0x0) ,
+ fh3EGenERecN(0x0)
+ ////////////////
{
// Default constructor
for(int ij = 0;ij<kMaxJets;++ij){
- fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
- fh2PtFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = 0;
- fh3PtRecGenHard[ij] = fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
- }
+ fh1E[ij] = fh1PtRecIn[ij] = fh1PtRecOut[ij] = fh1PtGenIn[ij] = fh1PtGenOut[ij] = 0;
+ fh2PtFGen[ij] = fh2Frag[ij] = fh2FragLn[ij] = 0;
+ fh3PtRecGenHard[ij] = fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
+ }
+ for(int ic = 0;ic < kMaxCorrelation;ic++){
+ fhnCorrelation[ic] = 0;
+ }
}
fh1PtHard_Trials(0x0),
fh1NGenJets(0x0),
fh1NRecJets(0x0),
- fHistList(0x0)
+ fHistList(0x0) ,
+ ////////////////
+ fh1JetMultiplicity(0x0) ,
+ fh2ERecZRec(0x0) ,
+ fh2EGenZGen(0x0) ,
+ fh2Efficiency(0x0) ,
+ fh3EGenERecN(0x0)
+ ////////////////
{
// Default constructor
for(int ij = 0;ij<kMaxJets;++ij){
fh3PtRecGenHard[ij] = fh3PtRecGenHard_NoW[ij] = fh3RecEtaPhiPt[ij] = fh3RecEtaPhiPt_NoGen[ij] =fh3GenEtaPhiPt_NoFound[ij] = fh3GenEtaPhiPt[ij] = 0;
}
+ for(int ic = 0;ic < kMaxCorrelation;ic++){
+ fhnCorrelation[ic] = 0;
+ }
+
DefineOutput(1, TList::Class());
}
nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt);
}
+
+ /////////////////////////////////////////////////////////////////
+ fh1JetMultiplicity = new TH1F("fh1JetMultiplicity", "Jet Multiplicity", 51, 0., 50.);
- const Int_t saveLevel = 2; // large save level more histos
+ fh2ERecZRec = new TH2F("fh2ERecZRec", " ; E^{jet}_{rec} [GeV]; z^{lp}_{rec}", 100, 0., 250., 100, 0., 2.);
+ fh2EGenZGen = new TH2F("fh2EGenZGen", " ; E^{jet}_{gen} [GeV]; z^{lp}_{gen}", 100, 0., 250., 100, 0., 2.);
+ fh2Efficiency = new TH2F("fh2Efficiency", "ERec/EGen;E^{jet}_{gen} [GeV];E^{jet}_{rec}/E^{jet}_{gen}", 100, 0., 250., 100, 0., 1.5);
+ fh3EGenERecN = new TH3F("fh3EGenERecN", "Efficiency vs. Jet Multiplicity", 100., 0., 250., 100, 0., 250., 51, 0., 50.);
+
+ // Response map
+ //arrays for bin limits
+ const Int_t nbin[4] = {100, 100, 100, 100};
+ Double_t LowEdge[4] = {0.,0.,0.,0.};
+ Double_t UpEdge[4] = {250., 250., 1., 1.};
+
+ for(int ic = 0;ic < kMaxCorrelation;ic++){
+ fhnCorrelation[ic] = new THnSparseF(Form("fhnCorrelation_%d",ic), "Response Map", 4, nbin, LowEdge, UpEdge);
+ if(ic==0) fhnCorrelation[ic]->SetTitle(Form("ResponseMap 0 <= npart <= %.0E",fgkJetNpartCut[ic]));
+ else fhnCorrelation[ic]->SetTitle(Form("ResponseMap %.0E < npart <= %.0E",fgkJetNpartCut[ic-1],fgkJetNpartCut[ic]));
+ }
+ const Int_t saveLevel = 1; // large save level more histos
if(saveLevel>0){
fHistList->Add(fh1PtHard);
fHistList->Add(fh1PtHard_NoW);
fHistList->Add(fh1PtHard_Trials);
fHistList->Add(fh1NGenJets);
fHistList->Add(fh1NRecJets);
+ ////////////////////////
+ fHistList->Add(fh1JetMultiplicity);
+ fHistList->Add(fh2ERecZRec);
+ fHistList->Add(fh2EGenZGen);
+ fHistList->Add(fh2Efficiency);
+ fHistList->Add(fh3EGenERecN);
+
+ for(int ic = 0;ic < kMaxCorrelation;++ic){
+ fHistList->Add(fhnCorrelation[ic]);
+ }
+ ////////////////////////
for(int ij = 0;ij<kMaxJets;++ij){
fHistList->Add(fh1E[ij]);
fHistList->Add(fh1PtRecIn[ij]);
Int_t nGenJets = 0;
AliAODJet recJets[kMaxJets];
Int_t nRecJets = 0;
+ ///////////////////////////
+ Int_t nTracks = 0;
+ //////////////////////////
Double_t eventW = 1;
Double_t ptHard = 0;
nRecJets = aodRecJets->GetEntries();
fh1NRecJets->Fill(nRecJets);
nRecJets = TMath::Min(nRecJets,kMaxJets);
+ //////////////////////////////////////////
+ nTracks = fAOD->GetNumberOfTracks();
+ ///////////////////////////////////////////
for(int ir = 0;ir < nRecJets;++ir){
AliAODJet *tmp = dynamic_cast<AliAODJet*>(aodRecJets->At(ir));
fh3RecEtaPhiPt[ir]->Fill(etaRec,phiRec,ptRec,eventW);
// Fill Correlation
Int_t ig = iGenIndex[ir];
- if(ig>=0&&ig<nGenJets){
+ if(ig>=0 && ig<nGenJets){
fh1PtRecOut[ir]->Fill(ptRec,eventW);
- Double_t ptGen = genJets[ig].Pt();
+ Double_t ptGen = genJets[ig].Pt();
+ Double_t phiGen = genJets[ig].Phi();
+ if(phiGen<0)phiGen+=TMath::Pi()*2.;
+ Double_t etaGen = genJets[ig].Eta();
fh2PtFGen[ir]->Fill(ptRec,ptGen,eventW);
fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW);
fh3PtRecGenHard_NoW[ir]->Fill(ptRec,ptGen,ptHard,1);
- }
+ /////////////////////////////////////////////////////
+ Double_t ERec = recJets[ir].E();
+ Double_t EGen = genJets[ig].E();
+ fh2Efficiency->Fill(EGen, ERec/EGen);
+ if (EGen>=0. && EGen<=250.){
+ Double_t Eleading = -1;
+ Double_t ptleading = -1;
+ Int_t nPart=0;
+ // loop over tracks
+ for (Int_t it = 0; it< nTracks; it++){
+ if (fAOD->GetTrack(it)->E() > EGen) continue;
+ Double_t phiTrack = fAOD->GetTrack(it)->Phi();
+ if (phiTrack<0) phiTrack+=TMath::Pi()*2.;
+ Double_t etaTrack = fAOD->GetTrack(it)->Eta();
+ Float_t deta = etaRec - etaTrack;
+ Float_t dphi = TMath::Abs(phiRec - phiTrack);
+ Float_t R = TMath::Sqrt(deta*deta + dphi*dphi);
+ // find leading particle
+ if (R<0.4 && fAOD->GetTrack(it)->E()>Eleading){
+ Eleading = fAOD->GetTrack(it)->E();
+ ptleading = fAOD->GetTrack(it)->Pt();
+ }
+ if (fAOD->GetTrack(it)->Pt()>0.03*EGen && fAOD->GetTrack(it)->E()<=EGen && R<0.7)
+ nPart++;
+ }
+ // fill Response Map (4D histogram) and Energy vs z distributions
+ Double_t var[4] = {EGen, ERec, ptleading/EGen, ptleading/ERec};
+ fh2ERecZRec->Fill(var[1],var[3]); // this has to be filled always in the real case...
+ fh2EGenZGen->Fill(var[0],var[2]);
+ fh1JetMultiplicity->Fill(nPart);
+ fh3EGenERecN->Fill(EGen, ERec, nPart);
+ for(int ic = 0;ic <kMaxCorrelation;ic++){
+ if (nPart<=fgkJetNpartCut[ic]){
+ fhnCorrelation[ic]->Fill(var);
+ break;
+ }
+ }
+ }
+ }
+ ////////////////////////////////////////////////////
else{
fh3RecEtaPhiPt_NoGen[ir]->Fill(etaRec,phiRec,ptRec,eventW);
}
* See cxx source for full Copyright notice */
#include "AliAnalysisTaskSE.h"
+////////////////
+#include <THnSparse.h>
+////////////////
class AliJetHeader;
class AliESDEvent;
class AliAODEvent;
//
enum {kAnaMC = 0x1};
+ enum {kMaxJets = 5};
+ enum {kMaxCorrelation = 3};
private:
AliAnalysisTaskJetSpectrum(const AliAnalysisTaskJetSpectrum&);
AliAnalysisTaskJetSpectrum& operator=(const AliAnalysisTaskJetSpectrum&);
+ static const Float_t fgkJetNpartCut[kMaxCorrelation];
- enum {kMaxJets = 5};
AliJetHeader *fJetHeaderRec;
AliJetHeader *fJetHeaderGen;
// ============================================
TList *fHistList; // Output list
+
+ ///////// For 2 dimensional unfolding //////////////////
+ TH1F* fh1JetMultiplicity;
+ TH2F* fh2ERecZRec;
+ TH2F* fh2EGenZGen;
+ TH2F* fh2Efficiency;
+ TH3F* fh3EGenERecN;
+ THnSparseF* fhnCorrelation[kMaxCorrelation];
+ ///////////////////////////////////////////////////////
ClassDef(AliAnalysisTaskJetSpectrum, 1) // Analysis task for standard jet analysis
fCurrentCorrelation = (THnSparseF*)fCorrelation->Clone("fCurrentCorrelation");
fCurrentCorrelation->Sumw2();
+ Printf("Correlation Matrix has %.0E filled bins", fCurrentCorrelation->GetNbins());
+
if (createBigBin)
{
Int_t maxBinE = 0, maxBinZ = 0;
//______________________________________________________________________________
-void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 10000)
+void AnalysisTrainCAF(Int_t nEvents = 10000, Int_t nOffset = 0, char *ds = "/PWG4/kleinb/LHC08v_jetjet15-50")
{
// Example of running analysis train in CAF. To run in debug mode:
// - export ROOTSYS=debug on your local client
if (iESDfilter) iAODhandler=1;
// Dataset from CAF
- // TString dataset = "/PWG4/kleinb/LHC08v_jetjet15-50";
- TString dataset = "/PWG4/kleinb/LHC08r_jetjet50";
+ TString dataset(ds);
// CKB quick hack for local analysis
- gROOT->LoadMacro("CreateESDChain.C");
- TChain *chain = CreateESDChain("jetjet15-50.txt",1000);
-
+ // gROOT->LoadMacro("CreateESDChain.C");
+ // TChain *chain = CreateESDChain("jetjet15-50.txt",1000);
+ TChain *chain = 0;
printf("==================================================================\n");
printf("=========== RUNNING ANALYSIS TRAIN IN CAF MODE =============\n");
TProof::Open(proofNode);
// Clear packages if changing ROOT version on CAF or local
- // gProof->ClearPackages();
+ // gProof->ClearPackages();
// Enable proof debugging if needed
// gProof->SetLogLevel(5);
// To debug the train in PROOF mode, type in a root session:
AliAODHandler* aodHandler = new AliAODHandler();
aodHandler->SetFillAOD(kFALSE);
mgr->SetOutputEventHandler(aodHandler);
- aodHandler->SetOutputFileName(Form("AliAODs_CKB_%07d-%07d.root",nOffset,nOffset+nEvents));
+ aodHandler->SetOutputFileName(Form("AliAODs_pwg4_%07d-%07d.root",nOffset,nOffset+nEvents));
cout_aod = mgr->CreateContainer("cAOD", TTree::Class(),
AliAnalysisManager::kOutputContainer, "default");
cout_aod->SetSpecialOutput();
mgr->AddTask(jetana);
// Output histograms list for jet analysis
AliAnalysisDataContainer *cout_jet = mgr->CreateContainer("jethist", TList::Class(),
- AliAnalysisManager::kOutputContainer,"jethist.root");
+ AliAnalysisManager::kOutputContainer,Form("jethist_%07d-%07d.root",nOffset,nOffset+nEvents));
// Dummy AOD output container for jet analysis (no client yet)
c_aodjet = mgr->CreateContainer("cAODjet", TTree::Class(),
AliAnalysisManager::kExchangeContainer);
mgr->AddTask(jetanaAOD);
// Output histograms list for jet analysis
AliAnalysisDataContainer *cout_jetAOD = mgr->CreateContainer("jethistAOD", TList::Class(),
- AliAnalysisManager::kOutputContainer, "jethistAOD.root");
+ AliAnalysisManager::kOutputContainer, Form("jethistAOD_%07d-%07d.root",nOffset,nOffset+nEvents));
// Dummy AOD output container for jet analysis (no client yet)
c_aodjet0 = mgr->CreateContainer("cAODjet0", TTree::Class(),
AliAnalysisManager::kExchangeContainer);
mgr->AddTask(jetanaMC);
// Output histograms list for jet analysis
AliAnalysisDataContainer *cout_jetMC = mgr->CreateContainer("jethistMC", TList::Class(),
- AliAnalysisManager::kOutputContainer, "jethistMC.root");
+ AliAnalysisManager::kOutputContainer,Form("jethistMC_%07d-%07d.root",nOffset,nOffset+nEvents));
// Dummy AOD output container for jet analysis (no client yet)
c_aodjetMC = mgr->CreateContainer("cAODjetMC", TTree::Class(),
AliAnalysisManager::kExchangeContainer);
// pwg4spec->SetBranchGen("jetsMC");
mgr->AddTask(pwg4spec);
- AliAnalysisDataContainer *coutput1_Spec = mgr->CreateContainer("pwg4spec", TList::Class(),AliAnalysisManager::kOutputContainer, "histos_pwg4spec.root");
+ AliAnalysisDataContainer *coutput1_Spec = mgr->CreateContainer("pwg4spec", TList::Class(),AliAnalysisManager::kOutputContainer,Form( "pwg4spec_%07d-%07d.root",nOffset,nOffset+nEvents));
// coutput1_Spec->SetSpecialOutput();
// Dummy AOD output container for jet analysis (no client yet)
c_aodSpec = mgr->CreateContainer("cAODjetSpec", TTree::Class(),
mgr->AddTask(ueana);
- AliAnalysisDataContainer *coutput1_UE = mgr->CreateContainer("histosUE", TList::Class(),AliAnalysisManager::kOutputContainer, "histosUE.root");
+ AliAnalysisDataContainer *coutput1_UE = mgr->CreateContainer("histosUE", TList::Class(),AliAnalysisManager::kOutputContainer, Form("pwg4UE_%07d-%07d.root",nOffset,nOffset+nEvents));
mgr->ConnectInput (ueana, 0, cinput);
// mgr->ConnectInput (ueana, 0, c_aodjet);
AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
- TFile::Open(fileNameRec);
+ TFile::Open(fileNameGen);
jetSpec->LoadHistograms();
- TFile::Open(fileNameGen);
- TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
- jetSpec->SetGenSpectrum(hist);
+ TFile::Open(fileNameRec);
+ TH2F* hist = (TH2F*) gFile->Get("unfolding/fRecSpectrum");
+ jetSpec->SetRecSpectrum(hist);
jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0);
// last parameter = calculateErrors <- this method to calculate the errors takes a lot of time
}
//_______________________________________________________________________________________________________________
-void FillSpecFromFile(const char* fileNameSpec = "histos_pwg4spec.root")
+void FillSpecFromFiles(const char* fileNameReal = "histos_pwg4spec.root",const char* fileNameSim = "histos_pwg4spec.root")
{
// This functions avoids problems when the number of bins or the bin limits
// in the histograms of the AliJetSpectrumUnfolding and AliAnalysisTaskJetSpectrum classes
gSystem->Load("libJETAN");
gSystem->Load("libPWG4JetTasks");
- file = new TFile(fileNameSpec,"read");
+ file = new TFile(fileNameSim,"read");
tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
- THnSparseF *fhCorrelation = (THnSparseF*)(tlist->FindObject("fhCorrelation_less5tracks"));
- THnSparseF *fhCorrelation2 = (THnSparseF*)(tlist->FindObject("fhCorrelation_5to10tracks"));
- THnSparseF *fhCorrelation3 = (THnSparseF*)(tlist->FindObject("fhCorrelation_more10tracks"));
- TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fhEGenZGen"));
- TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fhERecZRec"));
- fhCorrelation->Add(fhCorrelation2, 1);
- fhCorrelation->Add(fhCorrelation3, 1);
-
- delete fhCorrelation2;
- delete fhCorrelation3;
+ THnSparseF *fhCorrelation = 0;
+ for(int ic = 0;ic<3;++ic){
+ THnSparseF *hTmp = (THnSparseF*)(tlist->FindObject(Form("fhnCorrelation_%d",ic)));
+ if(fhCorrelation==0)fhCorrelation = (THnSparseF*)hTmp->Clone("fhnCorrelation");
+ else fhCorrelation->Add(hTmp);
+ }
+ TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen"));
AliJetSpectrumUnfolding *jetSpec = new AliJetSpectrumUnfolding("unfolding","unfolding");
jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
}
- file = TFile::Open("gen_pwg4spec.root", "RECREATE");
- jetSpec->SaveHistograms();
- file->Close();
- jetSpec->GetGenSpectrum()->Reset();
- printf("true distribution has been set\n");
- // reconstructed jets (measured distribution)
- for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++)
- for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++)
- {
- Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me);
- Float_t zm = fhERecZRec->GetYaxis()->GetBinCenter(mz);
- Int_t bine = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec);
- Int_t binz = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm);
- Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz);
- Float_t err = jetSpec->GetRecSpectrum()->GetBinError(bine, binz);
- jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz));
- jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
- }
- // Response function
- jetSpec->GetCorrelation()->Reset();
- jetSpec->GetCorrelation()->Sumw2();
-
+ Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++)
{
//printf("idx: %d\n",idx);
jetSpec->GetCorrelation()->SetBinContent(bin, cont + fhCorrelation->GetBinContent(idx));
jetSpec->GetCorrelation()->SetBinError(bin, err + fhCorrelation->GetBinError(idx));
}
+ // need number of entries for correct drawing
+ jetSpec->GetCorrelation()->SetEntries(fhCorrelation->GetEntries());
+
+
+ file = TFile::Open("gen_pwg4spec.root", "RECREATE");
+ jetSpec->SaveHistograms();
+ Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
+ file->Close();
+ jetSpec->GetGenSpectrum()->Reset();
+ printf("true distribution has been set\n");
+
+ // reconstructed jets (measured distribution)
+
+
+ // Response function
+ jetSpec->GetCorrelation()->Reset();
+ jetSpec->GetCorrelation()->Sumw2();
+ jetSpec->GetGenSpectrum()->Reset();
+ jetSpec->GetRecSpectrum()->Reset();
+
+
+ file = new TFile(fileNameReal,"read");
+ tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
+
+
+ TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fh2ERecZRec"));
+
+ for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++)
+ for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++)
+ {
+ Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me);
+ Float_t zm = fhERecZRec->GetYaxis()->GetBinCenter(mz);
+ Int_t bine = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec);
+ Int_t binz = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm);
+ Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz);
+ Float_t err = jetSpec->GetRecSpectrum()->GetBinError(bine, binz);
+ jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz));
+ jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
+ }
file = TFile::Open("rec_pwg4spec.root", "RECREATE");
jetSpec->SaveHistograms();
void correct(){
// simple steering to correct a given distribution;
loadlibs();
- FillSpecFromFile("/home/ydelgado/pcalice014.cern.ch/macros/jets/CAFdata/histos_pwg4spec.root");
+ FillSpecFromFiles("pwg4spec_0000000-0010000.root","pwg4spec_15-50.root");
char name[100];
sprintf(name, "unfolded_pwg4spec.root");
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff