histo binning standardization, some updates for MC
authordsilverm <dsilverm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 10 Sep 2010 13:46:50 +0000 (13:46 +0000)
committerdsilverm <dsilverm@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 10 Sep 2010 13:46:50 +0000 (13:46 +0000)
PWG4/totEt/AliAnalysisEt.cxx
PWG4/totEt/AliAnalysisEtCuts.h
PWG4/totEt/AliAnalysisEtMonteCarlo.cxx
PWG4/totEt/AliAnalysisEtMonteCarlo.h
PWG4/totEt/AliAnalysisEtReconstructed.cxx
PWG4/totEt/macros/runCaloEt.C

index c7db4fa..61503f3 100644 (file)
@@ -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);
index f2c8720..e8385db 100644 (file)
@@ -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;
index 50c479e..6fa9581 100644 (file)
@@ -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;
+}
+
index 7ea6db6..8e4d3da 100644 (file)
@@ -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
index 81d22ab..08cb28e 100644 (file)
@@ -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;
 
index a0d278e..e801cf9 100644 (file)
@@ -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;
   }