AliCutTask uses global cut initialization
authorjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Jan 2008 14:10:30 +0000 (14:10 +0000)
committerjgrosseo <jgrosseo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 22 Jan 2008 14:10:30 +0000 (14:10 +0000)
PWG0/esdTrackCuts/AliCutTask.cxx
PWG0/esdTrackCuts/AliCutTask.h
PWG0/esdTrackCuts/CreateCuts.C [deleted file]
PWG0/esdTrackCuts/draw.C [new file with mode: 0644]
PWG0/esdTrackCuts/run.C

index f64b40d..b408fa4 100644 (file)
 ClassImp(AliCutTask)
 
 //________________________________________________________________________
-AliCutTask::AliCutTask(const char *name) 
-  : AliAnalysisTask(name, ""), fESD(0), fTrackCuts(0), fTrackCutsPrimaries(0),
-    fTrackCutsSecondaries(0), fVertex(0), fTriggerStats(0), fOutput(0)
+AliCutTask::AliCutTask(const char *name) : 
+  AliAnalysisTask(name, ""), fESD(0), fTrackCuts(0), 
+  fAnalysisMode(AliPWG0Helper::kTPCITS),
+  fTrackCutsPrimaries(0),
+  fTrackCutsSecondaries(0), fVertex(0), fTriggerStats(0), fOutput(0),
+  fPrimStats(0)
 {
   // Constructor
 
@@ -50,7 +53,7 @@ void AliCutTask::ConnectInputData(Option_t *)
     Printf("ERROR: Could not read chain from input slot 0");
   } else {
     // Disable all branches and enable only the needed ones
-    tree->SetBranchStatus("*", kFALSE);
+    //tree->SetBranchStatus("*", kFALSE);
 
     // old esd
     tree->SetBranchStatus("fTriggerMask", kTRUE);
@@ -60,8 +63,8 @@ void AliCutTask::ConnectInputData(Option_t *)
     // new esd
     tree->SetBranchStatus("TriggerMask", kTRUE);
     tree->SetBranchStatus("AliESDHeader", kTRUE);
-    AliPWG0Helper::SetBranchStatusRecursive(tree, "Tracks", kTRUE);
-    AliPWG0Helper::SetBranchStatusRecursive(tree, "SPDVertex", kTRUE);
+    //AliPWG0Helper::SetBranchStatusRecursive(tree, "Tracks", kTRUE);
+    //AliPWG0Helper::SetBranchStatusRecursive(tree, "SPDVertex", kTRUE);
 
     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 
@@ -98,6 +101,16 @@ void AliCutTask::CreateOutputObjects()
   fTriggerStats->GetXaxis()->SetBinLabel(5, "VZERO_OR_LEFT | VZERO_OR_RIGHT");
   fOutput->Add(fTriggerStats);
 
+  fPrimStats = new TH1F("fPrimStats", "fPrimStats", 8, 0.5, 8.5);
+  fPrimStats->GetXaxis()->SetBinLabel(1, "Primary accepted once");
+  fPrimStats->GetXaxis()->SetBinLabel(2, "Primary accepted more than once");
+  fPrimStats->GetXaxis()->SetBinLabel(3, "Primary accepted more than once (count)");
+  fPrimStats->GetXaxis()->SetBinLabel(4, "Primary track rejected");
+  fPrimStats->GetXaxis()->SetBinLabel(5, "Primary track rejected (count)");
+  fPrimStats->GetXaxis()->SetBinLabel(6, "Primary track rejected, but other track accepted");
+  fPrimStats->GetXaxis()->SetBinLabel(7, "Primary track rejected, but other track accepted (count)");
+  fOutput->Add(fPrimStats);
+
   // disable info messages of AliMCEvent (per event)
   AliLog::SetClassDebugLevel("AliMCEvent", AliLog::kWarning - AliLog::kDebug + 1);
 }
@@ -133,15 +146,22 @@ void AliCutTask::Exec(Option_t *)
   if (fESD->GetTriggerMask() & 1 || fESD->GetTriggerMask() & 2)
     fTriggerStats->Fill(4);
 
-  //if (!AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask()), AliPWG0Helper::kMB1)
-  if (fESD->GetTriggerMask() & (0x1 << 14) == 0)
+  //if (fESD->GetTriggerMask() & (0x1 << 14) == 0)
+  if (!AliPWG0Helper::IsEventTriggered(fESD->GetTriggerMask(), AliPWG0Helper::kMB2))
     return;
 
-  if (!AliPWG0Helper::IsVertexReconstructed(fESD->GetVertex()))
+  // get the ESD vertex
+  const AliESDVertex* vtxESD = AliPWG0Helper::GetVertex(fESD, fAnalysisMode);
+  
+  if (!vtxESD)
+  {
+    Printf("No vertex. Skipping event");
     return;
+  }
 
-  //Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
-  fTrackCuts->CountAcceptedTracks(fESD);
+  Printf("There are %d tracks in this event", fESD->GetNumberOfTracks());
+  Int_t acceptedTracks = fTrackCuts->CountAcceptedTracks(fESD);
+  Printf("Accepted %d tracks", acceptedTracks);
   
   if (fTrackCutsPrimaries && fTrackCutsSecondaries)
   {
@@ -163,7 +183,17 @@ void AliCutTask::Exec(Option_t *)
       Printf("ERROR: Stack not available");
       return;
     }
-    
+
+    // count if primaries are counted several times
+    Int_t max = stack->GetNprimary();
+    Int_t* primAcc = new Int_t[max];
+    Int_t* primRej = new Int_t[max];
+    for (Int_t i=0; i<max; i++)
+    {
+      primAcc[i] = 0;
+      primRej[i] = 0;
+    }
+
     for (Int_t trackId = 0; trackId < fESD->GetNumberOfTracks(); trackId++)
     {
       AliESDtrack* track = fESD->GetTrack(trackId);
@@ -182,7 +212,19 @@ void AliCutTask::Exec(Option_t *)
       Int_t label = TMath::Abs(track->GetLabel());
       if (stack->IsPhysicalPrimary(label) == kTRUE)
       {
-        fTrackCutsPrimaries->AcceptTrack(track);
+       if (label >= max)
+       {
+         Printf("Warning label %d is higher than number of primaries %d", label, max);
+         continue;
+       }
+           
+        if (fTrackCutsPrimaries->AcceptTrack(track)) 
+       {
+         primAcc[label]++;
+       } 
+       else
+         primRej[label]++;
+
       }
       else
       {
@@ -208,10 +250,32 @@ void AliCutTask::Exec(Option_t *)
         fTrackCutsSecondaries->AcceptTrack(track);
       }
     }
+
+    for (Int_t i=0; i<max; i++) {
+      if (primAcc[i] == 1) {
+       fPrimStats->Fill(1);
+      } else if (primAcc[i] > 1) {
+       fPrimStats->Fill(2);
+       fPrimStats->Fill(3, primAcc[i]);
+      }
+
+      if (primRej[i] > 0) {
+       if (primAcc[i] == 0) {
+         fPrimStats->Fill(4);
+         fPrimStats->Fill(5, primRej[i]);
+       } else if (primAcc[i] > 0) {
+         fPrimStats->Fill(6);
+         fPrimStats->Fill(7, primRej[i]);
+       }
+      }
+    }
+
+    delete[] primAcc;
+    delete[] primRej;
   }
 
   // get the ESD vertex
-  fVertex->Fill(fESD->GetVertex()->GetZv());
+  fVertex->Fill(vtxESD->GetZv());
 }
 
 //________________________________________________________________________
@@ -258,6 +322,9 @@ void AliCutTask::Terminate(Option_t *)
 
   if (fTrackCutsSecondaries)
     fTrackCutsSecondaries->SaveHistograms();
+
+  if (fPrimStats)
+    fPrimStats->Write();
   
   file->Close();
   
index 238225e..a7b6c0c 100644 (file)
@@ -9,6 +9,7 @@ class AliESDEvent;
 class TList;
 
 #include "AliAnalysisTask.h"
+#include "AliPWG0Helper.h"
 
 class AliCutTask : public AliAnalysisTask {
  public:
@@ -21,11 +22,13 @@ class AliCutTask : public AliAnalysisTask {
   virtual void   Terminate(Option_t *);
 
   void SetTrackCuts(AliESDtrackCuts* cuts) { fTrackCuts = cuts; }
+  void SetAnalysisMode(AliPWG0Helper::AnalysisMode mode) { fAnalysisMode = mode; }
   void EnableSecondaryStudy();
   
  private:
   AliESDEvent *fESD;           //! ESD object
   AliESDtrackCuts* fTrackCuts; // track cuts
+  AliPWG0Helper::AnalysisMode fAnalysisMode; // detector that is used for analysis
   
   AliESDtrackCuts* fTrackCutsPrimaries; // cuts for tracks from primary particles
   AliESDtrackCuts* fTrackCutsSecondaries; // cuts for tracks from secondary particles
@@ -35,6 +38,8 @@ class AliCutTask : public AliAnalysisTask {
 
   TList* fOutput;                  //! list send on output slot 0
 
+  TH1F* fPrimStats; //! statistics about primaries, see bin names in CreateOutputData
+
   AliCutTask(const AliCutTask&); // not implemented
   AliCutTask& operator=(const AliCutTask&); // not implemented
   
diff --git a/PWG0/esdTrackCuts/CreateCuts.C b/PWG0/esdTrackCuts/CreateCuts.C
deleted file mode 100644 (file)
index e0e4435..0000000
+++ /dev/null
@@ -1,22 +0,0 @@
-/* $Id$ */
-
-// this macro creates the track and event cuts used in this analysis
-
-AliESDtrackCuts* CreateTrackCuts(Bool_t hists = kTRUE)
-{
-  AliESDtrackCuts* esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
-
-  if (hists)
-    esdTrackCuts->DefineHistograms(1);
-
-  esdTrackCuts->SetMinNClustersTPC(50);
-  esdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
-  esdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
-  esdTrackCuts->SetRequireTPCRefit(kTRUE);
-
-  esdTrackCuts->SetMinNsigmaToVertex(3);
-  esdTrackCuts->SetRequireSigmaToVertex(kTRUE);
-  esdTrackCuts->SetAcceptKingDaughters(kFALSE);
-
-  return esdTrackCuts;
-}
diff --git a/PWG0/esdTrackCuts/draw.C b/PWG0/esdTrackCuts/draw.C
new file mode 100644 (file)
index 0000000..d458a72
--- /dev/null
@@ -0,0 +1,68 @@
+void draw()
+{
+  const char* files[] = { "trackCuts_normal.root", "trackCuts_increased.root", "trackCuts_decreased.root" };
+  const char* titles[] = { "default geometry", "+ 10% material", "- 10% material" };
+  Int_t colors[] = { 1, 2, 4 };
+
+  TCanvas* c = new TCanvas;
+
+  TLegend* legend = new TLegend(0.7, 0.7, 0.9, 0.9);
+
+  for (Int_t i=0; i<3; i++) {
+    TFile::Open(files[i]);
+    
+    TH1* ptPrim = gFile->Get("fTrackCutsPrimaries/before_cuts/pt");
+    TH1* ptPrimCut = gFile->Get("fTrackCutsPrimaries/after_cuts/pt_cut");
+
+    TH1* ptSec = gFile->Get("fTrackCutsSecondaries/before_cuts/pt");
+    TH1* ptSecCut = gFile->Get("fTrackCutsSecondaries/after_cuts/pt_cut");
+
+    ptPrim->Add(ptSec);
+    ptPrimCut->Add(ptSecCut);
+
+    ptPrim->Sumw2();
+    ptPrimCut->Sumw2();
+    ptSec->Sumw2();
+    ptSecCut->Sumw2();
+    
+    Printf("%s", titles[i]);
+    Printf("%.2f %% secondaries before cuts", 100.0 * ptSec->GetEntries() / ptPrim->GetEntries());
+    Printf("%.2f %% secondaries after cuts", 100.0 * ptSecCut->GetEntries() / ptPrimCut->GetEntries());
+    Printf("");
+
+    ptSec->Divide(ptSec, ptPrim, 1, 1, "B");
+    ptSecCut->Divide(ptSecCut, ptPrimCut, 1, 1, "B");
+
+    ptSec->SetLineColor(colors[i]);
+    ptSecCut->SetLineColor(colors[i]);
+    ptSec->SetStats(kFALSE);
+
+    ptSec->GetXaxis()->SetRangeUser(0, 2);
+    ptSec->GetYaxis()->SetRangeUser(0, 1);
+
+    ptSec->SetTitle("");
+    ptSec->GetYaxis()->SetTitle("N_{Secondaries} / N_{All}");
+    
+    ptSec->DrawCopy((i > 0) ? "SAME" : "");
+    ptSecCut->DrawCopy("SAME");
+
+    legend->AddEntry(ptSec, titles[i]);
+  }
+
+  legend->Draw();
+
+  c->SaveAs("secondaries_changedmaterial.gif");
+}
+
+void drawStats(const char* fileName = "trackCuts.root")
+{
+  TFile::Open(fileName);
+
+  TH1* stat1 = gFile->Get("fTrackCutsPrimaries/cut_statistics");
+  TH1* stat2 = gFile->Get("fTrackCutsSecondaries/cut_statistics");
+  
+  new TCanvas;
+  stat1->Draw();
+  stat2->SetLineColor(2);
+  stat2->Draw("SAME");
+}
index f4b201a..70cef46 100644 (file)
@@ -34,15 +34,6 @@ void run(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug = kFALSE, B
   // Create the analysis manager
   mgr = new AliAnalysisManager("testAnalysis");
 
-  // selection of esd tracks
-  gROOT->ProcessLine(".L CreateCuts.C");
-  AliESDtrackCuts* esdTrackCuts = CreateTrackCuts();
-  if (!esdTrackCuts)
-  {
-    printf("ERROR: esdTrackCuts could not be created\n");
-    return;
-  }
-
   TString taskName("AliCutTask.cxx+");
   if (aDebug)
     taskName += "+g";
@@ -54,7 +45,24 @@ void run(Char_t* data, Int_t nRuns=20, Int_t offset=0, Bool_t aDebug = kFALSE, B
     gROOT->Macro(taskName);
 
   task = new AliCutTask;
-  task->SetTrackCuts(esdTrackCuts);
+
+  AliPWG0Helper::AnalysisMode analysisMode = AliPWG0Helper::kTPC;
+  task->SetAnalysisMode(analysisMode);
+
+  if (analysisMode != AliPWG0Helper::kSPD)
+  {
+    // selection of esd tracks
+    gROOT->ProcessLine(".L ../CreateStandardCuts.C");
+    AliESDtrackCuts* esdTrackCuts = CreateTrackCuts(analysisMode);
+    if (!esdTrackCuts)
+    {
+      printf("ERROR: esdTrackCuts could not be created\n");
+      return;
+    }
+
+    task->SetTrackCuts(esdTrackCuts);
+  }
+
   mgr->AddTask(task);
 
   // Add ESD handler