From 13b0d3c1e6adff3493b2629e379c8915b9f301bb Mon Sep 17 00:00:00 2001 From: dsilverm Date: Fri, 10 Sep 2010 13:46:50 +0000 Subject: [PATCH] histo binning standardization, some updates for MC --- PWG4/totEt/AliAnalysisEt.cxx | 47 ++++++++++++++--------- PWG4/totEt/AliAnalysisEtCuts.h | 2 +- PWG4/totEt/AliAnalysisEtMonteCarlo.cxx | 36 +++++++++++++++-- PWG4/totEt/AliAnalysisEtMonteCarlo.h | 8 +++- PWG4/totEt/AliAnalysisEtReconstructed.cxx | 11 +++--- PWG4/totEt/macros/runCaloEt.C | 13 ++++--- 6 files changed, 80 insertions(+), 37 deletions(-) diff --git a/PWG4/totEt/AliAnalysisEt.cxx b/PWG4/totEt/AliAnalysisEt.cxx index c7db4fa732f..61503f3da95 100644 --- a/PWG4/totEt/AliAnalysisEt.cxx +++ b/PWG4/totEt/AliAnalysisEt.cxx @@ -128,79 +128,88 @@ void AliAnalysisEt::Init() void AliAnalysisEt::CreateHistograms() { + // histogram binning for E_T, p_T and Multiplicity: defaults for p+p + Int_t nbinsEt = 1000; + Double_t minEt = 0.0001; + Double_t maxEt = 100; + Int_t nbinsPt = 200; + Double_t minPt = 0; + Double_t maxPt = 20; + Int_t nbinsMult = 200; + Double_t minMult = -0.5; // offset -0.5 to have integer bins centered around 0 + Double_t maxMult = nbinsMult + minMult; // 1 bin per integer value TString histname = "fHistEt" + fHistogramNameSuffix; - - fHistEt = new TH1F(histname.Data(), "Total E_{T} Distribution", 1000, 0.00, 99); + fHistEt = new TH1F(histname.Data(), "Total E_{T} Distribution", nbinsEt, minEt, maxEt); fHistEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})"); fHistEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)"); histname = "fHistChargedEt" + fHistogramNameSuffix; - fHistChargedEt = new TH1F(histname.Data(), "Total Charged E_{T} Distribution", 1000, 0.00, 99); + fHistChargedEt = new TH1F(histname.Data(), "Total Charged E_{T} Distribution", nbinsEt, minEt, maxEt); fHistChargedEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})"); fHistChargedEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)"); histname = "fHistNeutralEt" + fHistogramNameSuffix; - fHistNeutralEt = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution", 1000, 0.00, 99); + fHistNeutralEt = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution", nbinsEt, minEt, maxEt); fHistNeutralEt->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})"); fHistNeutralEt->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)"); histname = "fHistEtAcc" + fHistogramNameSuffix; - fHistEtAcc = new TH1F(histname.Data(), "Total E_{T} Distribution in Acceptance", 1000, 0.00, 99); + fHistEtAcc = new TH1F(histname.Data(), "Total E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt); fHistEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})"); fHistEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)"); histname = "fHistChargedEtAcc" + fHistogramNameSuffix; - fHistChargedEtAcc = new TH1F(histname.Data(), "Total Charged E_{T} Distribution in Acceptance", 1000, 0.00, 99); + fHistChargedEtAcc = new TH1F(histname.Data(), "Total Charged E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt); fHistChargedEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})"); fHistChargedEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)"); histname = "fHistNeutralEtAcc" + fHistogramNameSuffix; - fHistNeutralEtAcc = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution in Acceptance", 1000, 0.00, 99); + fHistNeutralEtAcc = new TH1F(histname.Data(), "Total Neutral E_{T} Distribution in Acceptance", nbinsEt, minEt, maxEt); fHistNeutralEtAcc->GetXaxis()->SetTitle("E_{T} (GeV/c^{2})"); fHistNeutralEtAcc->GetYaxis()->SetTitle("dN/dE_{T} (c^{2}/GeV)"); std::cout << histname << std::endl; histname = "fHistMult" + fHistogramNameSuffix; - fHistMult = new TH1F(histname.Data(), "Total Multiplicity", 200, 0, 199); + fHistMult = new TH1F(histname.Data(), "Total Multiplicity", nbinsMult, minMult, maxMult); fHistMult->GetXaxis()->SetTitle("N"); fHistMult->GetYaxis()->SetTitle("Multiplicity"); histname = "fHistChargedMult" + fHistogramNameSuffix; - fHistChargedMult = new TH1F(histname.Data(), "Charged Multiplicity", 200, 0, 199); + fHistChargedMult = new TH1F(histname.Data(), "Charged Multiplicity", nbinsMult, minMult, maxMult); fHistChargedMult->GetXaxis()->SetTitle("N"); fHistChargedMult->GetYaxis()->SetTitle("Multiplicity"); histname = "fHistNeutralMult" + fHistogramNameSuffix; - fHistNeutralMult = new TH1F(histname.Data(), "Charged Multiplicity", 200, 0, 199); + fHistNeutralMult = new TH1F(histname.Data(), "Neutral Multiplicity", nbinsMult, minMult, maxMult); fHistNeutralMult->GetXaxis()->SetTitle("N"); fHistNeutralMult->GetYaxis()->SetTitle("Multiplicity"); histname = "fHistPhivsPtPos" + fHistogramNameSuffix; - fHistPhivsPtPos = new TH2F(histname.Data(), "Phi vs pT of positively charged tracks hitting the calorimeter", 200, 0, 2*TMath::Pi(), 2000, 0, 100); + fHistPhivsPtPos = new TH2F(histname.Data(), "Phi vs pT of positively charged tracks hitting the calorimeter", 200, 0, 2*TMath::Pi(), nbinsPt, minPt, maxPt); histname = "fHistPhivsPtNeg" + fHistogramNameSuffix; - fHistPhivsPtNeg = new TH2F(histname.Data(), "Phi vs pT of negatively charged tracks hitting the calorimeter", 200, 0, 2*TMath::Pi(), 2000, 0, 100); + fHistPhivsPtNeg = new TH2F(histname.Data(), "Phi vs pT of negatively charged tracks hitting the calorimeter", 200, 0, 2*TMath::Pi(), nbinsPt, minPt, maxPt); histname = "fHistBaryonEt" + fHistogramNameSuffix; - fHistBaryonEt = new TH1F(histname.Data(), "E_{T} for baryons", 1000, 0.0001, 100); + fHistBaryonEt = new TH1F(histname.Data(), "E_{T} for baryons", nbinsEt, minEt, maxEt); histname = "fHistAntiBaryonEt" + fHistogramNameSuffix; - fHistAntiBaryonEt = new TH1F(histname.Data(), "E_{T} for anti baryons", 1000, 0.0001, 100); + fHistAntiBaryonEt = new TH1F(histname.Data(), "E_{T} for anti baryons", nbinsEt, minEt, maxEt); histname = "fHistMesonEt" + fHistogramNameSuffix; - fHistMesonEt = new TH1F(histname.Data(), "E_{T} for mesons", 1000, 0.0001, 100); + fHistMesonEt = new TH1F(histname.Data(), "E_{T} for mesons", nbinsEt, minEt, maxEt); histname = "fHistBaryonEtAcc" + fHistogramNameSuffix; - fHistBaryonEtAcc = new TH1F(histname.Data(), "E_{T} for baryons in calorimeter acceptance", 1000, 0.0001, 100); + fHistBaryonEtAcc = new TH1F(histname.Data(), "E_{T} for baryons in calorimeter acceptance", nbinsEt, minEt, maxEt); histname = "fHistAntiBaryonEtAcc" + fHistogramNameSuffix; - fHistAntiBaryonEtAcc = new TH1F(histname.Data(), "E_{T} for anti baryons in calorimeter acceptance", 1000, 0.0001, 100); + fHistAntiBaryonEtAcc = new TH1F(histname.Data(), "E_{T} for anti baryons in calorimeter acceptance", nbinsEt, minEt, maxEt); histname = "fHistMesonEtAcc" + fHistogramNameSuffix; - fHistMesonEtAcc = new TH1F(histname.Data(), "E_{T} for mesons in calorimeter acceptance", 1000, 0.0001, 100); + fHistMesonEtAcc = new TH1F(histname.Data(), "E_{T} for mesons in calorimeter acceptance", nbinsEt, minEt, maxEt); histname = "fHistEtRecvsEtMC" + fHistogramNameSuffix; - fHistEtRecvsEtMC = new TH2F(histname.Data(), "Reconstructed E_{t} vs MC E_{t}", 1000, 0.000, 100, 1000, 0.0001, 100); + fHistEtRecvsEtMC = new TH2F(histname.Data(), "Reconstructed E_{t} vs MC E_{t}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); histname = "fHistTMDeltaR" + fHistogramNameSuffix; fHistTMDeltaR = new TH1F(histname.Data(), "#Delta R for calorimeter clusters", 200, 0, 50); diff --git a/PWG4/totEt/AliAnalysisEtCuts.h b/PWG4/totEt/AliAnalysisEtCuts.h index f2c87207a80..e8385dbd82d 100644 --- a/PWG4/totEt/AliAnalysisEtCuts.h +++ b/PWG4/totEt/AliAnalysisEtCuts.h @@ -58,7 +58,7 @@ namespace EtReconstructedCutsEmcal { const Char_t kClusterType = -1; - const Double_t kClusterEnergyCut = 0.0; + const Double_t kClusterEnergyCut = 0.1; // GeV const Double_t kSingleCellEnergyCut = 0.5; const Double_t kTrackDistanceCut = 15.0; diff --git a/PWG4/totEt/AliAnalysisEtMonteCarlo.cxx b/PWG4/totEt/AliAnalysisEtMonteCarlo.cxx index 50c479e0251..6fa95816385 100644 --- a/PWG4/totEt/AliAnalysisEtMonteCarlo.cxx +++ b/PWG4/totEt/AliAnalysisEtMonteCarlo.cxx @@ -1,8 +1,9 @@ #include "AliAnalysisEtMonteCarlo.h" #include "AliAnalysisEtCuts.h" - +#include "AliESDtrack.h" #include "AliStack.h" #include "AliMCEvent.h" +#include "TH2F.h" Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev) { @@ -43,7 +44,7 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev) if (TMath::Abs(part->Eta()) < fEtaCut) { - TParticlePDG *pdgCode = part->GetPDG(0); + TParticlePDG *pdgCode = part->GetPDG(0); if ( TMath::Abs(pdgCode->PdgCode()) == ProtonCode || TMath::Abs(pdgCode->PdgCode()) == NeutronCode || @@ -76,8 +77,15 @@ Int_t AliAnalysisEtMonteCarlo::AnalyseEvent(AliVEvent* ev) fTotChargedEtAcc += part->Energy()*TMath::Sin(part->Theta()); fTotEtAcc += part->Energy()*TMath::Sin(part->Theta()); } - } - } + + // if (TrackHitsCalorimeter(part, event->GetMagneticField())) + if (TrackHitsCalorimeter(part)) // magnetic field info not filled? + { + if (pdgCode->Charge() > 0) fHistPhivsPtPos->Fill(part->Phi(),part->Pt()); + else fHistPhivsPtNeg->Fill(part->Phi(), part->Pt()); + } + } + } } fTotNeutralEtAcc = fTotNeutralEt; @@ -101,3 +109,23 @@ void AliAnalysisEtMonteCarlo::Init() fIPzCut = EtReconstructedCuts::kIPzCut; } + +bool AliAnalysisEtMonteCarlo::TrackHitsCalorimeter(TParticle* part, Double_t magField) +{ + // printf(" TrackHitsCalorimeter - magField %f\n", magField); + AliESDtrack *esdTrack = new AliESDtrack(part); + // Printf("MC Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt()); + + Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField); + + // if(prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt()); + + bool status = prop && + TMath::Abs(esdTrack->Eta()) < fEtaCutAcc && + esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. && + esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.; + delete esdTrack; + + return status; +} + diff --git a/PWG4/totEt/AliAnalysisEtMonteCarlo.h b/PWG4/totEt/AliAnalysisEtMonteCarlo.h index 7ea6db6dfec..8e4d3dae8e5 100644 --- a/PWG4/totEt/AliAnalysisEtMonteCarlo.h +++ b/PWG4/totEt/AliAnalysisEtMonteCarlo.h @@ -2,7 +2,7 @@ #define ALIANALYSISETMONTECARLO_H #include "AliAnalysisEt.h" - +#include "TParticle.h" class AliAnalysisEtMonteCarlo : public AliAnalysisEt { @@ -14,7 +14,11 @@ public: virtual Int_t AnalyseEvent(AliVEvent* event); virtual void Init(); - + +protected: + + virtual bool TrackHitsCalorimeter(TParticle *part, Double_t magField=0.5); + }; #endif // ALIANALYSISETMONTECARLO_H diff --git a/PWG4/totEt/AliAnalysisEtReconstructed.cxx b/PWG4/totEt/AliAnalysisEtReconstructed.cxx index 81d22aba2e4..08cb28e7599 100644 --- a/PWG4/totEt/AliAnalysisEtReconstructed.cxx +++ b/PWG4/totEt/AliAnalysisEtReconstructed.cxx @@ -74,12 +74,12 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev) } } - Double_t phi = track->Phi(); - Double_t pt = track->Pt(); if (TrackHitsCalorimeter(track, event->GetMagneticField())) { - if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi,pt); - else fHistPhivsPtNeg->Fill(phi, pt); + Double_t phi = track->Phi(); + Double_t pt = track->Pt(); + if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt); + else fHistPhivsPtNeg->Fill(phi, pt); } } @@ -123,12 +123,13 @@ Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev) float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2; // float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) ); fTotNeutralEt += cluster->E() * TMath::Sin(theta); + fNeutralMultiplicity++; fMultiplicity++; - } fTotNeutralEtAcc = fTotNeutralEt; + fTotEt = fTotChargedEt + fTotNeutralEt; fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc; diff --git a/PWG4/totEt/macros/runCaloEt.C b/PWG4/totEt/macros/runCaloEt.C index a0d278e3b91..e801cf98903 100644 --- a/PWG4/totEt/macros/runCaloEt.C +++ b/PWG4/totEt/macros/runCaloEt.C @@ -5,7 +5,7 @@ //With the argument true this submits jobs to the grid //As written this requires an xml script tag.xml in the ~/et directory on the grid to submit jobs void runCaloEt(bool submit = true, // true or false - const char *dataType="sim", // "sim" or "real" or "simPbPb" + const char *dataType="simPbPb", // "sim" or "real" or "simPbPb" const char *det = "EMCAL") // "PHOS" or "EMCAL" { TStopwatch timer; @@ -39,7 +39,6 @@ void runCaloEt(bool submit = true, // true or false gROOT->ProcessLine(".L AliAnalysisTaskTotEt.cxx+g"); } - char *kTreeName = "esdTree" ; TChain * chain = new TChain(kTreeName,"myESDTree") ; @@ -62,7 +61,7 @@ void runCaloEt(bool submit = true, // true or false cout << " taskName " << taskName << " outputName " << outputName << endl; - if(submit){ + if (submit) { gROOT->LoadMacro("CreateAlienHandlerCaloEtSim.C"); AliAnalysisGrid *alienHandler = CreateAlienHandlerCaloEtSim(outputDir, outputName); if (!alienHandler) return; @@ -74,16 +73,18 @@ void runCaloEt(bool submit = true, // true or false AliMCEventHandler* handler = new AliMCEventHandler; if ( dataStr.Contains("sim") ) { cout << " MC " << endl; - chain->Add("/home/dsilverm/data/E_T/sim/LHC10d1/117222/100/AliESDs.root"); // link to local test file if ( dataStr.Contains("PbPb") ) { cout << " PbPb " << endl; chain->Add("/home/dsilverm/data/E_T/sim/LHC10e11/191001/001/AliESDs.root"); // link to local test file } + else { // pp + chain->Add("/home/dsilverm/data/E_T/sim/LHC10d1/117222/100/AliESDs.root"); // link to local test file + } handler->SetReadTR(kFALSE); mgr->SetMCtruthEventHandler(handler); } - else { - chain->Add("/home/dsilverm/data/E_T/data/2010/LHC10b/000117222/ESDs/pass2/10000117222021.30/AliESDs.root"); // link to local test file + else { // real data + chain->Add("/home/dsilverm/data/E_T/data/2010/LHC10b/000117222/ESDs/pass2/10000117222021.30/AliESDs.root"); // link to local test file cout << " not MC " << endl; } -- 2.43.0