]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
First version of an analysis code for Spectra analysis on AODs. The code was develope...
authormfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 30 Sep 2011 09:49:00 +0000 (09:49 +0000)
committermfloris <mfloris@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 30 Sep 2011 09:49:00 +0000 (09:49 +0000)
16 files changed:
PWG2/SPECTRA/SpectraAOD/AddTaskSpectraAOD.C [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/AliAnalysisTaskSpectraAOD.cxx [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/AliAnalysisTaskSpectraAOD.h [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/AliSpectraAODEventCuts.cxx [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/AliSpectraAODEventCuts.h [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/AliSpectraAODHistoManager.cxx [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/AliSpectraAODHistoManager.h [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/AliSpectraAODTrackCuts.cxx [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/AliSpectraAODTrackCuts.h [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/README [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/analysis_macro.C [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/run.C [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/runAOD.C [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/runAODProof.C [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/runAODProof_init.C [new file with mode: 0644]
PWG2/SPECTRA/SpectraAOD/runProof.C [new file with mode: 0644]

diff --git a/PWG2/SPECTRA/SpectraAOD/AddTaskSpectraAOD.C b/PWG2/SPECTRA/SpectraAOD/AddTaskSpectraAOD.C
new file mode 100644 (file)
index 0000000..dab181f
--- /dev/null
@@ -0,0 +1,53 @@
+AliAnalysisTaskSpectraAOD * AddTaskSpectraAOD(const char * outfilename)
+{
+  // TODO: add some parameters to set the centrality for this task, and maybe the name of the task
+  // TODO: shall I use the same file and different dirs for the different centralities?
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskSpectraAOD", "No analysis manager to connect to.");
+    return NULL;
+  }  
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler()) {
+    ::Error("AddTaskSpectraAOD", "This task requires an input event handler");
+    return NULL;
+  }
+  TString inputDataType = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+  
+  if (inputDataType != "AOD") {
+    Printf("ERROR! This task can only run on AODs!");
+  }
+
+  // Configure analysis
+  //===========================================================================
+    
+   
+  // Set I/O
+  AliAnalysisDataContainer *cinput0 = mgr->GetCommonInputContainer();
+  AliAnalysisTaskSpectraAOD *task = new AliAnalysisTaskSpectraAOD("TaskAODSpectra");
+  mgr->AddTask(task);
+  // Set the cuts
+  AliSpectraAODVertexCuts * vcuts = new AliSpectraAODVertexCuts("VertexCuts");
+  AliSpectraAODTrackCuts  * tcuts = new AliSpectraAODTrackCuts ("TracksCuts");
+  tcuts->SetTrackType(6);
+  tcuts->SetEta(1.);
+  task->SetVertexCuts(vcuts);
+  task->SetTrackCuts (tcuts);
+
+
+  AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+  AliAnalysisDataContainer *coutputpt1 = mgr->CreateContainer("chistpt", AliSpectraAODHistoManager::Class(),  AliAnalysisManager::kOutputContainer, outfilename);
+  AliAnalysisDataContainer *coutputpt2 = mgr->CreateContainer("cvcutpt", AliSpectraAODVertexCuts::Class(),    AliAnalysisManager::kOutputContainer, outfilename);
+  AliAnalysisDataContainer *coutputpt3 = mgr->CreateContainer("ctcutpt", AliSpectraAODTrackCuts::Class(),     AliAnalysisManager::kOutputContainer, outfilename);
+
+  mgr->ConnectInput(task, 0, cinput);
+  mgr->ConnectOutput(task, 1, coutputpt1);
+  mgr->ConnectOutput(task, 2, coutputpt2);
+  mgr->ConnectOutput(task, 3, coutputpt3);
+
+  return task;
+}   
diff --git a/PWG2/SPECTRA/SpectraAOD/AliAnalysisTaskSpectraAOD.cxx b/PWG2/SPECTRA/SpectraAOD/AliAnalysisTaskSpectraAOD.cxx
new file mode 100644 (file)
index 0000000..816bc4f
--- /dev/null
@@ -0,0 +1,249 @@
+/**************************************************************************\r
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
+ *                                                                        *\r
+ * Author: The ALICE Off-line Project.                                    *\r
+ * Contributors are mentioned in the code where appropriate.              *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+\r
+//-----------------------------------------------------------------\r
+//         AliAnalysisTaskSpectraAOD class\r
+//-----------------------------------------------------------------\r
+\r
+#include "TChain.h"\r
+#include "TTree.h"\r
+#include "TLegend.h"\r
+#include "TH1F.h"\r
+#include "TH2F.h"\r
+#include "TCanvas.h"\r
+#include "AliAnalysisTask.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliAODTrack.h"\r
+#include "AliAODMCParticle.h"\r
+#include "AliAODEvent.h"\r
+#include "AliAODInputHandler.h"\r
+#include "AliAnalysisTaskSpectraAOD.h"\r
+#include "AliAnalysisTaskESDfilter.h"\r
+#include "AliAnalysisDataContainer.h"\r
+#include "AliSpectraAODHistoManager.h"\r
+#include "AliSpectraAODTrackCuts.h"\r
+#include "AliSpectraAODEventCuts.h"\r
+#include "AliCentrality.h"\r
+#include "TProof.h"\r
+#include "AliPID.h"\r
+#include "AliVEvent.h"\r
+#include "AliPIDResponse.h"\r
+#include <iostream>\r
+\r
+using namespace AliSpectraNameSpace;\r
+using namespace std;\r
+\r
+ClassImp(AliAnalysisTaskSpectraAOD) // EX1 This stuff tells root to implement the streamer, inspector methods etc (we discussed about it today)\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskSpectraAOD::AliAnalysisTaskSpectraAOD(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0), fIsMC(0), fPIDResponse(0), fNSigmaPID(0), fYCut(0)\r
+{\r
+   // Default constructor\r
+   fNSigmaPID = 3.;\r
+   fYCut = .5;\r
+   DefineInput(0, TChain::Class());\r
+   DefineOutput(1, AliSpectraAODHistoManager::Class());\r
+   DefineOutput(2, AliSpectraAODEventCuts::Class());\r
+   DefineOutput(3, AliSpectraAODTrackCuts::Class());\r
+\r
+}\r
+//________________________________________________________________________\r
+Bool_t AliAnalysisTaskSpectraAOD::CheckYCut(AODParticleSpecies_t species, AliAODTrack* track) const\r
+{\r
+    // check if the rapidity is within the set range\r
+    // note: masses are hardcoded for now. we could look them up in the pdg database, but that would mean accecing it 100k+ times per run ...\r
+    Double_t y;\r
+    if (species == kProton) { y = track->Y(9.38271999999999995e-01); }\r
+    if ( species == kKaon ) { y = track->Y(4.93676999999999977e-01); }\r
+    if ( species == kPion)  { y = track->Y(1.39570000000000000e-01); }\r
+    if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;\r
+    return kTRUE;\r
+}\r
+//____________________________________________________________________________\r
+Bool_t AliAnalysisTaskSpectraAOD::CheckYCut(AliAODMCParticle* particle) const\r
+{\r
+    // check if the rapidity is within the set range\r
+    Double_t y = particle->Y();\r
+    if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;\r
+    return kTRUE;\r
+}\r
+//________________________________________________________________________\r
+void AliAnalysisTaskSpectraAOD::UserCreateOutputObjects()\r
+{\r
+  // create output objects\r
+    fHistMan = new AliSpectraAODHistoManager("SpectraHistos");\r
+\r
+   if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
+   if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
+\r
+   AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();\r
+   AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());\r
+   fPIDResponse = inputHandler->GetPIDResponse();\r
+\r
+   PostData(1, fHistMan);\r
+   PostData(2, fEventCuts);\r
+   PostData(3, fTrackCuts);\r
+\r
+}\r
+//________________________________________________________________________\r
+void AliAnalysisTaskSpectraAOD::UserExec(Option_t *)\r
+{\r
+  // main event loop\r
+   fAOD = dynamic_cast<AliAODEvent*>(fInputEvent);\r
+   if (strcmp(fAOD->ClassName(), "AliAODEvent"))\r
+   {\r
+      AliFatal("Not processing AODs");\r
+   }\r
+\r
+   if(!fPIDResponse) {\r
+     AliError("Cannot get pid response");\r
+     return;\r
+   }\r
+\r
+   if (!fEventCuts->IsSelected(fAOD)) return;\r
+\r
+   //AliCentrality fAliCentral*;\r
+//   if ((fAOD->GetCentrality())->GetCentralityPercentile("V0M") > 5.) return;\r
+   \r
+   // First do MC to fill up the MC particle array, such that we can use it later\r
+   TClonesArray *arrayMC = 0;\r
+   if (fIsMC)\r
+   {\r
+      arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
+      if (!arrayMC) {\r
+         AliFatal("Error: MC particles branch not found!\n");\r
+      }\r
+      Int_t nMC = arrayMC->GetEntries();\r
+      for (Int_t iMC = 0; iMC < nMC; iMC++)\r
+      {\r
+         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
+         fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt());\r
+         \r
+        if(TMath::Abs(partMC->Eta()) > fTrackCuts->GetEta()) continue;\r
+         // check for true PID + and fill P_t histos \r
+         if (partMC->IsPhysicalPrimary() && CheckYCut(partMC) ) {// only primary vertices and y cut satisfied\r
+            if ( partMC->PdgCode() == 2212) { fHistMan->GetHistogram(kHistPtGenTruePrimaryProtonPlus)->Fill(partMC->Pt()); } \r
+            if ( partMC->PdgCode() == -2212) { fHistMan->GetHistogram(kHistPtGenTruePrimaryProtonMinus)->Fill(partMC->Pt()); } \r
+            if ( partMC->PdgCode() == 321)  { fHistMan->GetHistogram(kHistPtGenTruePrimaryKaonPlus)->Fill(partMC->Pt()); } \r
+            if ( partMC->PdgCode() == -321) { fHistMan->GetHistogram(kHistPtGenTruePrimaryKaonMinus)->Fill(partMC->Pt()); } \r
+            if ( partMC->PdgCode() == 211)  { fHistMan->GetHistogram(kHistPtGenTruePrimaryPionPlus)->Fill(partMC->Pt()); } \r
+            if ( partMC->PdgCode() == -211) { fHistMan->GetHistogram(kHistPtGenTruePrimaryPionMinus)->Fill(partMC->Pt()); } \r
+         }\r
+      }\r
+   }\r
+   \r
+   for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++)\r
+   {\r
+\r
+      AliAODTrack* track = fAOD->GetTrack(iTracks);\r
+      if (!fTrackCuts->IsSelected(track)) continue;\r
+\r
+      fHistMan->GetHistogram(kHistPtRec)->Fill(track->Pt());  // PT histo\r
+      fHistMan->GetPIDHistogram(kHistPIDTPC)->Fill(track->P(), track->GetTPCsignal()); // PID histo\r
+      fHistMan->GetPIDHistogram(kHistPIDTPCPt)->Fill(track->Pt(), track->GetTPCsignal());\r
+      \r
+      // Sigma identify particles (for both MC and regular data):\r
+      // note: as of now, very crude identification. how many sigmas are allowed? what to do with high pt?\r
+      AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(track);\r
+      Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)); \r
+      Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton));\r
+      Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)); \r
+      if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton )) { \r
+       \r
+          if ((nsigmaTPCkKaon > fNSigmaPID) || (!CheckYCut(kKaon, track) ) ) continue;\r
+          if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaKaonPlus)->Fill(track->Pt()); } \r
+          else { fHistMan->GetHistogram(kHistPtRecSigmaKaonMinus)->Fill(track->Pt()); } \r
+      }\r
+      if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
+          if ( nsigmaTPCkProton > fNSigmaPID || (!CheckYCut(kProton, track) ) )  continue;\r
+          if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaProtonPlus)->Fill(track->Pt()); }\r
+          else { fHistMan->GetHistogram(kHistPtRecSigmaProtonMinus)->Fill(track->Pt()); }\r
+      }\r
+      if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
+          if (nsigmaTPCkPion > fNSigmaPID || (!CheckYCut(kPion, track) ) ) continue;\r
+          if ( track->Charge() > 0 )  { fHistMan->GetHistogram(kHistPtRecSigmaPionPlus)->Fill(track->Pt()); }\r
+          else  { fHistMan->GetHistogram(kHistPtRecSigmaPionMinus)->Fill(track->Pt()); }\r
+      }\r
+      /* MC Part */\r
+      if (arrayMC) {\r
+         AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
+         if (!partMC) { \r
+             AliError("Cannot get MC particle");\r
+             continue; }\r
+         if (CheckYCut(partMC)) {\r
+         // primaries, true pid\r
+         if ( partMC->PdgCode() == 2212) { fHistMan->GetHistogram(kHistPtRecTrueProtonPlus)->Fill(partMC->Pt()); \r
+            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryProtonPlus)->Fill(partMC->Pt()); }}\r
+         if ( partMC->PdgCode() == -2212) { fHistMan->GetHistogram(kHistPtRecTrueProtonMinus)->Fill(partMC->Pt()); \r
+            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryProtonMinus)->Fill(partMC->Pt()); }}\r
+         if ( partMC->PdgCode() == 321) { fHistMan->GetHistogram(kHistPtRecTrueKaonPlus)->Fill(partMC->Pt()); \r
+            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryKaonPlus)->Fill(partMC->Pt()); }}\r
+         if ( partMC->PdgCode() == -321) { fHistMan->GetHistogram(kHistPtRecTrueKaonMinus)->Fill(partMC->Pt()); \r
+            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryKaonMinus)->Fill(partMC->Pt()); }}\r
+         if ( partMC->PdgCode() == 211) { fHistMan->GetHistogram(kHistPtRecTruePionPlus)->Fill(partMC->Pt()); \r
+            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryPionPlus)->Fill(partMC->Pt()); }}\r
+         if ( partMC->PdgCode() == -211) { fHistMan->GetHistogram(kHistPtRecTruePionMinus)->Fill(partMC->Pt());  \r
+            if (partMC->IsPhysicalPrimary()) {fHistMan->GetHistogram(kHistPtRecTruePrimaryPionMinus)->Fill(partMC->Pt()); }}\r
+         }\r
+\r
+         // primaries, sigma pid\r
+         if (partMC->IsPhysicalPrimary()) { \r
+            if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) {\r
+               if ( (nsigmaTPCkKaon > fNSigmaPID ) || (!CheckYCut(kKaon, track) ) ) continue; \r
+               if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryKaonPlus)->Fill(track->Pt()); } \r
+               else { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryKaonMinus)->Fill(track->Pt()); } \r
+            }\r
+            if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
+               if ( (nsigmaTPCkProton > fNSigmaPID ) || (!CheckYCut(kProton, track) ) ) continue;\r
+               if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryProtonPlus)->Fill(track->Pt()); }\r
+               else { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryProtonMinus)->Fill(track->Pt()); }\r
+            }\r
+            if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
+               if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
+               if ( track->Charge() > 0 )  { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryPionPlus)->Fill(track->Pt()); }\r
+               else  { fHistMan->GetHistogram(kHistPtRecSigmaPrimaryPionMinus)->Fill(track->Pt()); }\r
+            }\r
+         }\r
+         //secondaries\r
+         else {\r
+            if( ( nsigmaTPCkKaon < nsigmaTPCkPion ) && ( nsigmaTPCkKaon < nsigmaTPCkProton ) ) { \r
+               if ( (nsigmaTPCkKaon > fNSigmaPID )  || (!CheckYCut(kKaon, track) ) ) continue;\r
+               if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryKaonPlus)->Fill(track->Pt()); } \r
+               else { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryKaonMinus)->Fill(track->Pt()); } \r
+            }\r
+            if( ( nsigmaTPCkProton < nsigmaTPCkKaon ) && ( nsigmaTPCkProton < nsigmaTPCkPion ) ) {\r
+               if ( (nsigmaTPCkProton > fNSigmaPID )  || (!CheckYCut(kProton, track) ) ) continue;\r
+               if ( track->Charge() > 0 ) { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryProtonPlus)->Fill(track->Pt()); }\r
+               else { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryProtonMinus)->Fill(track->Pt()); }\r
+            }\r
+            if( (nsigmaTPCkPion < nsigmaTPCkProton ) && ( nsigmaTPCkPion < nsigmaTPCkKaon ) ) {\r
+               if ( ( nsigmaTPCkPion > fNSigmaPID )  || (!CheckYCut(kPion, track) ) ) continue;\r
+               if ( track->Charge() > 0 )  { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryPionPlus)->Fill(track->Pt()); }\r
+               else  { fHistMan->GetHistogram(kHistPtRecSigmaSecondaryPionMinus)->Fill(track->Pt()); }\r
+            }\r
+         }\r
+      }\r
+   }\r
+\r
+   PostData(1, fHistMan);\r
+   PostData(2, fEventCuts);\r
+   PostData(3, fTrackCuts);\r
+}\r
+  \r
+//_________________________________________________________________\r
+void   AliAnalysisTaskSpectraAOD::Terminate(Option_t *)\r
+{\r
+   // Terminate\r
+}\r
diff --git a/PWG2/SPECTRA/SpectraAOD/AliAnalysisTaskSpectraAOD.h b/PWG2/SPECTRA/SpectraAOD/AliAnalysisTaskSpectraAOD.h
new file mode 100644 (file)
index 0000000..e04139a
--- /dev/null
@@ -0,0 +1,65 @@
+#ifndef ALIANALYSISTASKSPECTRAAOD_H\r
+#define ALIANALYSISTASKSPECTRAAOD_H\r
+\r
+/*  See cxx source for full Copyright notice */\r
+\r
+//-------------------------------------------------------------------------\r
+//                      AliAnalysisTaskSpectraAOD\r
+//\r
+//\r
+//\r
+//\r
+// Author: Michele Floris, CERN\r
+//-------------------------------------------------------------------------\r
+\r
+class TH1F;\r
+class TH2F;\r
+class AliAODEvent;\r
+class AliSpectraAODHistoManager;\r
+class AliSpectraAODTrackCuts;\r
+class AliSpectraAODEventCuts;\r
+#include "AliSpectraAODHistoManager.h"\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+\r
+class AliAnalysisTaskSpectraAOD : public AliAnalysisTaskSE\r
+{\r
+public:\r
+\r
+   // constructors\r
+   AliAnalysisTaskSpectraAOD() : AliAnalysisTaskSE(), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0), fIsMC(0), fPIDResponse(0), fNSigmaPID(0), fYCut(0) {}\r
+   AliAnalysisTaskSpectraAOD(const char *name);\r
+   virtual ~AliAnalysisTaskSpectraAOD() {}\r
+\r
+   void SetIsMC(Bool_t isMC = kFALSE)    {fIsMC = isMC; };\r
+   Bool_t GetIsMC()           const           { return fIsMC;};\r
+   void SetNSigmaForIdentification (Double_t sigma ) { fNSigmaPID = sigma; }\r
+   Double_t GetNSigmaForIdentification () const {return fNSigmaPID; }\r
+   void SetYCut (Double_t y ) { fYCut = y; }\r
+   Double_t GetYCut () const {return fYCut; }\r
+\r
+   virtual void   UserCreateOutputObjects();\r
+   Bool_t         CheckYCut(AliSpectraNameSpace::AODParticleSpecies_t species, AliAODTrack* track) const;\r
+   Bool_t         CheckYCut(AliAODMCParticle* particle) const;\r
+   virtual void   UserExec(Option_t *option);\r
+   virtual void   Terminate(Option_t *);\r
+   void SetTrackCuts(AliSpectraAODTrackCuts * tc)   {      fTrackCuts = tc;   }\r
+   void SetEventCuts(AliSpectraAODEventCuts * vc)   {      fEventCuts = vc;   }\r
+\r
+private:\r
+\r
+   AliAODEvent           * fAOD;         //! AOD object\r
+   AliSpectraAODHistoManager      * fHistMan;       // Histogram Manager\r
+   AliSpectraAODTrackCuts      * fTrackCuts;     // Track Cuts\r
+   AliSpectraAODEventCuts      * fEventCuts;     // Event Cuts\r
+   Bool_t          fIsMC;// true if processing MC\r
+   AliPIDResponse                        *fPIDResponse;     // ! PID response object\r
+   Double_t        fNSigmaPID; // Maximum number of sigmas allowed in particle identification\r
+   Double_t        fYCut; // Maximum rapidity - calculated from identified particle's mass\r
+   AliAnalysisTaskSpectraAOD(const AliAnalysisTaskSpectraAOD&);\r
+   AliAnalysisTaskSpectraAOD& operator=(const AliAnalysisTaskSpectraAOD&);\r
+\r
+   ClassDef(AliAnalysisTaskSpectraAOD, 1);\r
+};\r
+\r
+#endif\r
diff --git a/PWG2/SPECTRA/SpectraAOD/AliSpectraAODEventCuts.cxx b/PWG2/SPECTRA/SpectraAOD/AliSpectraAODEventCuts.cxx
new file mode 100644 (file)
index 0000000..d050f58
--- /dev/null
@@ -0,0 +1,138 @@
+
+/**************************************************************************
+ * Copyright(c) 1998-2009, 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.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//         AliSpectraAODEventCuts class
+//-----------------------------------------------------------------
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TLegend.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAnalysisTaskESDfilter.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliSpectraAODEventCuts.h"
+#include "AliSpectraAODHistoManager.h"
+#include <iostream>
+
+using namespace std;
+
+ClassImp(AliSpectraAODEventCuts)
+
+AliSpectraAODEventCuts::AliSpectraAODEventCuts(const char *name) : TNamed(name, "AOD Event Cuts"), fAOD(0), fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fHistoCuts(0)
+{
+   // Constructor
+   fHistoCuts = new TH1I("fEventCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
+   fCentralityCutMin = 0.0;      // default value of centrality cut minimum, 0 ~ no cut
+   fCentralityCutMax = 10000.0;  // default value of centrality cut maximum,  ~ no cut
+
+}
+
+//______________________________________________________
+Bool_t AliSpectraAODEventCuts::IsSelected(AliAODEvent * aod)
+{
+// Returns true if Event Cuts are selected and applied
+   fAOD = aod;
+   fIsSelected = (CheckVtxRange() && CheckCentralityCut());
+   if(fIsSelected)  fHistoCuts->Fill(kAcceptedEvents);
+   return fIsSelected;
+}
+
+//______________________________________________________
+Bool_t AliSpectraAODEventCuts::CheckVtxRange()
+{
+  // reject events outside of range
+   AliAODVertex * vertex = fAOD->GetPrimaryVertex();
+   if (!vertex)
+   {
+      fHistoCuts->Fill(kVtxNoEvent);
+      return kFALSE;
+   }
+   if (TMath::Abs(vertex->GetZ()) < 10)
+   {
+      return kTRUE;
+   }
+   fHistoCuts->Fill(kVtxRange);
+   return kFALSE;
+}
+
+//______________________________________________________
+Bool_t AliSpectraAODEventCuts::CheckCentralityCut()
+{
+   // Check centrality cut
+   if ( (fAOD->GetCentrality()->GetCentralityPercentile("V0M") <= fCentralityCutMax)  &&  (fAOD->GetCentrality()->GetCentralityPercentile("V0M") >= fCentralityCutMin) )  return kTRUE;   
+   fHistoCuts->Fill(kVtxCentral);
+   return kFALSE;
+}
+
+//______________________________________________________
+void AliSpectraAODEventCuts::PrintCuts()
+{
+    // print info about event cuts
+    cout << "Event Stats" << endl;
+    cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
+    cout << " > Vertex out of range: " << fHistoCuts->GetBinContent(kVtxRange + 1) << endl;
+    cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
+    cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
+    }
+//______________________________________________________
+
+Long64_t AliSpectraAODEventCuts::Merge(TCollection* list)
+{
+  // Merge a list of AliSpectraAODEventCuts objects with this.
+  // Returns the number of merged objects (including this).
+
+  //  AliInfo("Merging");
+
+  if (!list)
+    return 0;
+
+  if (list->IsEmpty())
+    return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj;
+
+  // collections of all histograms
+  TList collections;
+
+  Int_t count = 0;
+
+  while ((obj = iter->Next())) {
+    AliSpectraAODEventCuts* entry = dynamic_cast<AliSpectraAODEventCuts*> (obj);
+    if (entry == 0) 
+      continue;
+
+    TH1I * histo = entry->GetHistoCuts();      
+    collections.Add(histo);
+    count++;
+  }
+  
+  fHistoCuts->Merge(&collections);
+  
+  delete iter;
+
+  return count+1;
+}
+
diff --git a/PWG2/SPECTRA/SpectraAOD/AliSpectraAODEventCuts.h b/PWG2/SPECTRA/SpectraAOD/AliSpectraAODEventCuts.h
new file mode 100644 (file)
index 0000000..56d25ab
--- /dev/null
@@ -0,0 +1,61 @@
+#ifndef ALISPECTRAAODEVENTCUTS_H
+#define ALISPECTRAAODEVENTCUTS_H
+
+/*  See cxx source for full Copyright notice */
+
+//-------------------------------------------------------------------------
+//                      AliSpectraAODEventCuts
+//
+//
+//
+//
+// Authors: Michele Floris, CERN, Philip Versteeg, UU, Redmer Bertens, UU
+//-------------------------------------------------------------------------
+
+class AliAODEvent;
+class AliSpectraAODHistoManager;
+
+#include "TNamed.h"
+
+class AliSpectraAODEventCuts : public TNamed
+{
+public:
+   enum { kAcceptedEvents = 0, kVtxRange, kVtxCentral, kVtxNoEvent, kNVtxCuts};
+
+   // Constructors
+   AliSpectraAODEventCuts() : TNamed(), fAOD(0), fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fHistoCuts(0) {}
+   AliSpectraAODEventCuts(const char *name);
+   virtual  ~AliSpectraAODEventCuts() {}
+
+   // Methods
+   Bool_t IsSelected(AliAODEvent * aod);
+   Bool_t CheckVtxRange();
+   Bool_t CheckCentralityCut();
+   void  SetCentralityCutMin(Float_t cut)  { fCentralityCutMin = cut; }
+   void  SetCentralityCutMax(Float_t cut)  { fCentralityCutMax = cut; }
+
+   
+   TH1I * GetHistoCuts()         {  return fHistoCuts; }
+   Float_t  GetCentralityMin()  const {  return fCentralityCutMin; }
+   Float_t  GetCentralityMax()  const {  return fCentralityCutMax; }
+   void   PrintCuts();
+   Float_t  NumberOfEvents()     { return fHistoCuts->GetBinContent(kAcceptedEvents+1); }
+
+  Long64_t Merge(TCollection* list);
+
+
+private:
+
+   AliAODEvent     *fAOD;              //! AOD event
+   Bool_t          fIsSelected;        // True if cuts are selected
+   Float_t         fCentralityCutMin;     // minimum centrality percentile
+  Float_t         fCentralityCutMax;     // maximum centrality percentile
+   TH1I            *fHistoCuts;        // Cuts statistics
+   AliSpectraAODEventCuts(const AliSpectraAODEventCuts&);
+   AliSpectraAODEventCuts& operator=(const AliSpectraAODEventCuts&);
+
+   ClassDef(AliSpectraAODEventCuts, 2);
+
+};
+#endif
+
diff --git a/PWG2/SPECTRA/SpectraAOD/AliSpectraAODHistoManager.cxx b/PWG2/SPECTRA/SpectraAOD/AliSpectraAODHistoManager.cxx
new file mode 100644 (file)
index 0000000..fa7faae
--- /dev/null
@@ -0,0 +1,130 @@
+
+/**************************************************************************
+ * Copyright(c) 1998-2009, 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.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//         AliSpectraAODHistoManager class
+//-----------------------------------------------------------------
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TLegend.h"
+#include "TH1.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAnalysisTaskESDfilter.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliSpectraAODHistoManager.h"
+#include <iostream>
+
+using namespace std;
+
+ClassImp(AliSpectraAODHistoManager)
+
+
+using namespace AliSpectraNameSpace;
+
+AliSpectraAODHistoManager::AliSpectraAODHistoManager(const char *name): TNamed(name, "AOD Spectra Histo Manager"), fOutputList(0)
+{
+  // ctor
+   fOutputList = new TList;
+   fOutputList->SetOwner();
+   Bool_t oldStatus = TH1::AddDirectoryStatus();
+   TH1::AddDirectory(kFALSE);
+   for (Int_t ihist  = 0; ihist < kNHist ; ihist++)
+   {
+      if (ihist <= kNPtHist) BookPtHistogram(kHistName[ihist]);  // PT histos
+      if (ihist > kNPtHist) BookPIDHistogram(kHistName[ihist]);  // PID histos
+   }
+   TH1::AddDirectory(oldStatus);
+
+}
+//_______________________________________________________
+TH1F* AliSpectraAODHistoManager::BookPtHistogram(const char * name)
+{
+   // Return a pt histogram with predefined binning, set the ID and add it to the output list
+   AliInfo(Form("Booking pt histogram %s", name));
+
+   TH1F * hist = new TH1F(name, Form("P_{T} distribution (%s)", name), 40, 0, 1.2);
+   hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
+   hist->GetYaxis()->SetTitle("No. of events");
+   hist->SetMarkerStyle(kFullCircle);
+   hist->Sumw2();
+   fOutputList->Add(hist);
+
+   return hist;
+}
+//_____________________________________________________________________________
+
+TH2F* AliSpectraAODHistoManager::BookPIDHistogram(const char * name)
+{
+   // Return a pt histogram with predefined binning, set the ID and add it to the output list
+   AliInfo(Form("Booking pt histogram %s", name));
+
+   TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 100, 0, 1.2, 1000, 0, 1000);
+   hist->GetXaxis()->SetTitle("P (GeV / c)");
+   hist->GetYaxis()->SetTitle("TPC");
+//  hist->Sumw2();
+   fOutputList->Add(hist);
+
+   return hist;
+}
+
+
+Long64_t AliSpectraAODHistoManager::Merge(TCollection* list)
+{
+  // Merge a list of AliSpectraAODHistoManager objects with this.
+  // Returns the number of merged objects (including this).
+
+  //  AliInfo("Merging");
+
+  if (!list)
+    return 0;
+
+  if (list->IsEmpty())
+    return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj;
+
+  // collections of all histograms
+  TList collections;
+
+  Int_t count = 0;
+
+  while ((obj = iter->Next())) {
+    AliSpectraAODHistoManager* entry = dynamic_cast<AliSpectraAODHistoManager*> (obj);
+    if (entry == 0) 
+      continue;
+
+    TList * hlist = entry->GetOutputList();      
+    collections.Add(hlist);
+    count++;
+  }
+  
+  fOutputList->Merge(&collections);
+  
+  delete iter;
+
+  return count+1;
+}
+
diff --git a/PWG2/SPECTRA/SpectraAOD/AliSpectraAODHistoManager.h b/PWG2/SPECTRA/SpectraAOD/AliSpectraAODHistoManager.h
new file mode 100644 (file)
index 0000000..d875ce6
--- /dev/null
@@ -0,0 +1,204 @@
+
+#ifndef ALISPECTRAAODHISTOMANAGER_H
+#define ALISPECTRAAODHISTOMANAGER_H
+
+/*  See cxx source for full Copyright notice */
+
+//-------------------------------------------------------------------------
+//                      AliSpectraAODHistoManager
+//
+//
+//
+//
+// Authors: Michele Floris, CERN, Philip Versteeg, UU, Redmer Bertens, UU
+//-------------------------------------------------------------------------
+
+class AliAODEvent;
+class TH1F;
+class TH2F;
+class TList;
+#include "TNamed.h"
+
+namespace AliSpectraNameSpace
+{
+   enum AODPtHist_t
+   {
+      // 6 Pt Reconstructed Sigma
+      kHistPtRecSigmaProtonPlus = 0,            // Pt histo for protons +, reconstructed tracks, sigma ID
+      kHistPtRecSigmaKaonPlus,                  // Pt histo for kaons +, reconsructed tracks, sigma ID
+      kHistPtRecSigmaPionPlus,                  // Pt histo for pions +, reconstructed tracks, sigma ID
+      kHistPtRecSigmaProtonMinus,               // Pt histo for protons -, reconstructed tracks, sigma ID
+      kHistPtRecSigmaKaonMinus,                 // Pt histo for kaons -, reconstructed tracks, sigma ID
+      kHistPtRecSigmaPionMinus,                 // Pt histo for pions -, reconstructed tracks, sigma ID
+      
+      // 6 Pt Reconstructed True
+      kHistPtRecTrueProtonPlus,                 // Pt histo for protons +, reconstructed tracks, true ID
+      kHistPtRecTrueKaonPlus,                   // Pt histo for kaons +, reconsructed tracks, true ID
+      kHistPtRecTruePionPlus,                   // Pt histo for pions +, reconstructed tracks, true ID
+      kHistPtRecTrueProtonMinus,                // Pt histo for protons -, reconstructed tracks, true ID
+      kHistPtRecTrueKaonMinus,                  // Pt histo for kaons -, reconstructed tracks, true ID
+      kHistPtRecTruePionMinus,                  // Pt histo for pions -, reconstructed tracks, true ID
+            
+      // 6 Pt Reconstructed Sigma Primary
+      kHistPtRecSigmaPrimaryProtonPlus,         // Pt histo for protons +, reconstructed tracks, sigma ID, primary Event
+      kHistPtRecSigmaPrimaryKaonPlus,           // Pt histo for kaons +, reconsructed tracks, sigma ID, primary Event
+      kHistPtRecSigmaPrimaryPionPlus,           // Pt histo for pions +, reconstructed tracks, sigma ID, primary Event
+      kHistPtRecSigmaPrimaryProtonMinus,        // Pt histo for protons -, reconstructed tracks, sigma ID, primary Event
+      kHistPtRecSigmaPrimaryKaonMinus,          // Pt histo for kaons -, reconstructed tracks, sigma ID, primary Event
+      kHistPtRecSigmaPrimaryPionMinus,          // Pt histo for pions -, reconstructed tracks, sigma ID, primary Event
+            
+      // 6 Pt Reconstructed Sigma Secondary
+      kHistPtRecSigmaSecondaryProtonPlus,       // Pt histo for protons +, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryKaonPlus,         // Pt histo for kaons +, reconsructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryPionPlus,         // Pt histo for pions +, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryProtonMinus,      // Pt histo for protons -, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryKaonMinus,        // Pt histo for kaons -, reconstructed tracks, sigma ID, secondary Event
+      kHistPtRecSigmaSecondaryPionMinus,        // Pt histo for pions -, reconstructed tracks, sigma ID, secondary Event
+
+      // 6 Pt Generated True Primary
+      kHistPtGenTruePrimaryProtonPlus,          // Pt histo for protons +, generated tracks, true ID, primary Event
+      kHistPtGenTruePrimaryKaonPlus,            // Pt histo for kaons +, generated tracks, true ID, primary Event
+      kHistPtGenTruePrimaryPionPlus,            // Pt histo for pions +, generated tracks, true ID, primary Event
+      kHistPtGenTruePrimaryProtonMinus,         // Pt histo for protons -, generated tracks, true ID, primary Event
+      kHistPtGenTruePrimaryKaonMinus,           // Pt histo for kaons -, generated tracks, true ID, primary Event
+      kHistPtGenTruePrimaryPionMinus,           // Pt histo for pions -, generated tracks, true ID, primary Event
+      kNPtSpecies = kHistPtGenTruePrimaryPionMinus,
+      
+      // 6 Pt Reconstructed True Primary
+      kHistPtRecTruePrimaryProtonPlus,          // Pt histo for protons +, reconstructed tracks, true ID, primary event
+      kHistPtRecTruePrimaryKaonPlus,            // Pt histo for kaons +, reconsructed tracks, true ID, primary event
+      kHistPtRecTruePrimaryPionPlus,            // Pt histo for pions +, reconstructed tracks, true ID, primary event
+      kHistPtRecTruePrimaryProtonMinus,         // Pt histo for protons -, reconstructed tracks, true ID, primary event
+      kHistPtRecTruePrimaryKaonMinus,           // Pt histo for kaons -, reconstructed tracks, true ID, primary event
+      kHistPtRecTruePrimaryPionMinus,           // Pt histo for pions -, reconstructed tracks, true ID, primary event
+      
+      // Rest
+      kHistPtRec,                               // Pt histo for all particles, reconstructed tracks
+      kHistPtGen,                               // Pt histo for all particles, generated tracks
+      kNPtHist = kHistPtGen,                    // Number of pt-likehistos histos
+      kHistPIDTPC,                              // Particle Identification histo
+      kHistPIDTPCPt,
+      kNHist,                                   // Total number of histos
+   };  // Type of events plotted in Pt Histogram
+
+   const char * kHistName[] =
+   {
+      // 6 Pt Reconstructed Sigma
+      "histPtRecSigmaProtonPlus",               // Pt histo for protons +, reconstructed tracks, sigma ID
+      "histPtRecSigmaKaonPlus",                 // Pt histo for kaons +, reconsructed tracks, sigma ID
+      "histPtRecSigmaPionPlus",                 // Pt histo for pions +, reconstructed tracks, sigma ID
+      "histPtRecSigmaProtonMinus",              // Pt histo for protons -, reconstructed tracks, sigma ID
+      "histPtRecSigmaKaonMinus",                // Pt histo for kaons -, reconstructed tracks, sigma ID
+      "histPtRecSigmaPionMinus",                // Pt histo for pions -, reconstructed tracks, sigma ID
+      
+      // 6 Pt Reconstructed True
+      "histPtRecTrueProtonPlus",                // Pt histo for protons +, reconstructed tracks, true ID
+      "histPtRecTrueKaonPlus",                  // Pt histo for kaons +, reconsructed tracks, true ID
+      "histPtRecTruePionPlus",                  // Pt histo for pions +, reconstructed tracks, true ID
+      "histPtRecTrueProtonMinus",               // Pt histo for protons -, reconstructed tracks, true ID
+      "histPtRecTrueKaonMinus",                 // Pt histo for kaons -, reconstructed tracks, true ID
+      "histPtRecTruePionMinus",                 // Pt histo for pions -, reconstructed tracks, true ID
+
+      // 6 Pt Reconstructed Sigma Primary
+      "histPtRecSigmaPrimaryProtonPlus",        // Pt histo for protons +, reconstructed tracks, sigma ID, primary Event
+      "histPtRecSigmaPrimaryKaonPlus",          // Pt histo for kaons +, reconsructed tracks, sigma ID, primary Event
+      "histPtRecSigmaPrimaryPionPlus",          // Pt histo for pions +, reconstructed tracks, sigma ID, primary Event
+      "histPtRecSigmaPrimaryProtonMinus",       // Pt histo for protons -, reconstructed tracks, sigma ID, primary Event
+      "histPtRecSigmaPrimaryKaonMinus",         // Pt histo for kaons -, reconstructed tracks, sigma ID, primary Event
+      "histPtRecSigmaPrimaryPionMinus",         // Pt histo for pions -, reconstructed tracks, sigma ID, primary Event
+      
+      // 6 Pt Reconstructed Sigma Seconday
+      "histPtRecSigmaSecondaryProtonPlus",      // Pt histo for protons +, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryKaonPlus",        // Pt histo for kaons +, reconsructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryPionPlus",        // Pt histo for pions +, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryProtonMinus",     // Pt histo for protons -, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryKaonMinus",       // Pt histo for kaons -, reconstructed tracks, sigma ID, secondary Event
+      "histPtRecSigmaSecondaryPionMinus",       // Pt histo for pions -, reconstructed tracks, sigma ID, secondary Event
+      
+      // 6 Pt Reconstructed Sigma Primary
+      "histPtGenTruePrimaryProtonPlus",         // Pt histo for protons +, generated tracks, sigma ID, primary Event
+      "histPtGenTruePrimaryKaonPlus",           // Pt histo for kaons +, generated tracks, sigma ID, primary Event
+      "histPtGenTruePrimaryPionPlus",           // Pt histo for pions +, generated tracks, sigma ID, primary Event
+      "histPtGenTruePrimaryProtonMinus",          // Pt histo for protons -, generated tracks, sigma ID, primary Event
+      "histPtGenTruePrimaryKaonMinus",            // Pt histo for kaons -, generated tracks, sigma ID, primary Event
+      "histPtGenTruePrimaryPionMinus",            // Pt histo for pions -, generated tracks, sigma ID, primary Event
+
+      // 6 Pt Reconstructed True
+      "histPtRecTruePrimaryProtonPlus",         // Pt histo for protons +, reconstructed tracks, true ID, primary event
+      "histPtRecTruePrimaryKaonPlus",           // Pt histo for kaons +, reconsructed tracks, true ID, primary event
+      "histPtRecTruePrimaryPionPlus",           // Pt histo for pions +, reconstructed tracks, true ID, primary event
+      "histPtRecTruePrimaryProtonMinus",        // Pt histo for protons -, reconstructed tracks, true ID, primary event
+      "histPtRecTruePrimaryKaonMinus",          // Pt histo for kaons -, reconstructed tracks, true ID, primary event
+      "histPtRecTruePrimaryPionMinus",          // Pt histo for pions -, reconstructed tracks, true ID, primary event
+      
+      // Rest
+      "histPtRec",                              // Pt histo for all particles, reconstructed tracks
+      "histPtGen",                              // Pt histo for all particles, generated tracks
+      "histPIDTPC",                             // Particle Identification histo
+      "histPIDTPCPt",
+   };
+
+   enum AODParticleSpecies_t
+   {
+       kProton = 0,
+       kKaon,
+       kPion,
+       kNSpecies = kPion,
+   }; // Particle species used in plotting
+
+   const char * kParticleSpecies[] =
+   {
+       "ProtonPlus",
+       "ProtonMinus",
+       "KaonPlus",
+       "KaonMinus",
+       "PionPlus",
+       "PionMinus",
+   };
+   enum AODHistoType_t
+   {
+       kHistSpectraRec = 0,
+       kHistSpectraGen,
+       kHNHistoTypes = kHistSpectraGen,
+   }; // Types of histos
+
+}
+
+using namespace AliSpectraNameSpace;
+
+class AliSpectraAODHistoManager : public TNamed
+{
+public:
+   AliSpectraAODHistoManager() :  TNamed(), fOutputList(0) {}
+   AliSpectraAODHistoManager(const char *name);
+   virtual  ~AliSpectraAODHistoManager() {}
+
+
+   TH1F*   BookPtHistogram(const char * name);
+   TH2F*   BookPIDHistogram(const char * name);
+   TH1*     GetHistogram(UInt_t id)      {      return (TH1*) fOutputList->At(id);   }
+   TH1*     GetPtHistogram(UInt_t id)    {      return (TH1*) fOutputList->At(id);   }
+   TH1*     GetPtHistogram(const char * name)   {      return (TH1*) fOutputList->FindObject(name);   }
+   TH1*     GetPtHistogramByName(UInt_t id)     {      return (TH1*) fOutputList->FindObject(kHistName[id]); }  // Use this if you want to read a file saved with a different histo list   
+   TH2*     GetPIDHistogram(UInt_t id)   {      return (TH2*) fOutputList->At(id);   }
+   TH2*     GetPIDHistogram(const char * name)  {      return (TH2*) fOutputList->FindObject(name);   }
+   TH2*     GetPIDHistogramByName(UInt_t id)    {      return (TH2*) fOutputList->FindObject(kHistName[id]);  }// Use this if you want to read a file saved with a different histo list
+   TH1F*   GetTH1F(UInt_t id)            {      return (TH1F*) GetPtHistogram(id);   }
+   TH2F*   GetTH2F(UInt_t id)            {      return (TH2F*) GetPIDHistogram(id);   }
+
+  TList * GetOutputList() {return fOutputList;}
+
+  Long64_t Merge(TCollection* list);
+
+
+private:
+   TList     *fOutputList;  // List of Pt Histo's
+
+   AliSpectraAODHistoManager(const AliSpectraAODHistoManager&);
+   AliSpectraAODHistoManager& operator=(const AliSpectraAODHistoManager&);
+
+   ClassDef(AliSpectraAODHistoManager, 1);
+
+};
+#endif
+
diff --git a/PWG2/SPECTRA/SpectraAOD/AliSpectraAODTrackCuts.cxx b/PWG2/SPECTRA/SpectraAOD/AliSpectraAODTrackCuts.cxx
new file mode 100644 (file)
index 0000000..de04d4d
--- /dev/null
@@ -0,0 +1,167 @@
+
+/**************************************************************************
+ * Copyright(c) 1998-2009, 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.                  *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+//         AliSpectraAODTrackCuts class
+//-----------------------------------------------------------------
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TLegend.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAnalysisTaskESDfilter.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliSpectraAODTrackCuts.h"
+#include "AliSpectraAODHistoManager.h"
+#include <iostream>
+
+using namespace std;
+
+ClassImp(AliSpectraAODTrackCuts)
+
+
+AliSpectraAODTrackCuts::AliSpectraAODTrackCuts(const char *name) : TNamed(name, "AOD Track Cuts"), fIsSelected(0), fTrackBits(0), fEtaCut(0), fDCACut(0), fPCut(0), fPtCut(0), fHistoCuts(0), fTrack(0)
+
+{
+   // Constructor
+   fHistoCuts = new TH1I("fTrkCuts", "Track Cuts", kNTrkCuts, -0.5, kNTrkCuts - 0.5);
+   fEtaCut = 100000.0; // default value of eta cut ~ no cut
+   fDCACut = 100000.0; // default value of dca cut ~ no cut
+   fPCut = 100000.0; // default value of p cut ~ no cut
+   fPtCut = 100000.0; // default value of pt cut ~ no cut 
+   
+}
+
+//_______________________________________________________
+Bool_t AliSpectraAODTrackCuts::IsSelected(AliAODTrack * track)
+{
+// Returns true if Track Cuts are selected and applied
+   if (!track)
+   {
+      printf("ERROR: Could not receive track");
+      return kFALSE;
+   }
+   fTrack = track;
+   fIsSelected = (CheckTrackType() && CheckEtaCut() && CheckDCACut() && CheckPCut() && CheckPtCut());
+   return fIsSelected ;
+}
+//_________________________________________________________
+
+Bool_t AliSpectraAODTrackCuts::CheckTrackType()
+{
+   // Check track cuts
+   if (fTrack->TestFilterBit(fTrackBits)) return kTRUE;
+   fHistoCuts->Fill(kTrkBit);
+   return kFALSE;
+}
+//________________________________________________________
+Bool_t AliSpectraAODTrackCuts::CheckEtaCut()
+{
+   // Check eta cut
+   if (fTrack->Eta() < fEtaCut && fTrack->Eta() > - fEtaCut) return kTRUE;
+   fHistoCuts->Fill(kTrkEta);
+   return kFALSE;
+}
+//_______________________________________________________
+Bool_t AliSpectraAODTrackCuts::CheckDCACut()
+{
+   // Check DCA cut
+   if (fTrack->DCA() < fDCACut) return kTRUE;
+   fHistoCuts->Fill(kTrkDCA);
+   return kFALSE;
+}
+//________________________________________________________
+Bool_t AliSpectraAODTrackCuts::CheckPCut()
+{
+   // Check P cut
+   if (fTrack->P() < fPCut) return kTRUE;
+   fHistoCuts->Fill(kTrkP);
+   return kFALSE;
+}
+//_______________________________________________________
+Bool_t AliSpectraAODTrackCuts::CheckPtCut()
+{
+    // check Pt cut
+//    if ((fTrack->Pt() < fPtCut) && (fTrack->Pt() > 0.3 )) return kTRUE;
+   if (fTrack->Pt() < fPtCut) return kTRUE;
+    fHistoCuts->Fill(kTrkPt);
+    return kFALSE;
+}
+//_______________________________________________________
+void AliSpectraAODTrackCuts::PrintCuts() const
+{
+  // Print cuts
+    cout << "Track Cuts" << endl;
+    cout << " > TrackBit\t" << fTrackBits << endl;
+    cout << " > Eta cut\t" << fEtaCut << endl;
+    cout << " > DCA cut\t" << fDCACut << endl;
+    cout << " > P cut\t" << fPCut << endl;
+    cout << " > Pt cut \t" << fPtCut << endl;
+}
+//_______________________________________________________
+void AliSpectraAODTrackCuts::SetTrackType(UInt_t bit)
+{
+   // Set the type of track to be used. The argument should be the bit number. The mask is produced automatically.
+   fTrackBits = (0x1 << (bit - 1));
+}
+//_______________________________________________________
+
+Long64_t AliSpectraAODTrackCuts::Merge(TCollection* list)
+{
+  // Merge a list of AliSpectraAODTrackCuts objects with this.
+  // Returns the number of merged objects (including this).
+
+  //  AliInfo("Merging");
+
+  if (!list)
+    return 0;
+
+  if (list->IsEmpty())
+    return 1;
+
+  TIterator* iter = list->MakeIterator();
+  TObject* obj;
+
+  // collections of all histograms
+  TList collections;
+
+  Int_t count = 0;
+
+  while ((obj = iter->Next())) {
+    AliSpectraAODTrackCuts* entry = dynamic_cast<AliSpectraAODTrackCuts*> (obj);
+    if (entry == 0) 
+      continue;
+
+    TH1I * histo = entry->GetHistoCuts();      
+    collections.Add(histo);
+    count++;
+  }
+  
+  fHistoCuts->Merge(&collections);
+  
+  delete iter;
+
+  return count+1;
+}
+
diff --git a/PWG2/SPECTRA/SpectraAOD/AliSpectraAODTrackCuts.h b/PWG2/SPECTRA/SpectraAOD/AliSpectraAODTrackCuts.h
new file mode 100644 (file)
index 0000000..b00f92a
--- /dev/null
@@ -0,0 +1,74 @@
+#ifndef ALISPECTRAAODTRACKCUTS_H
+#define ALISPECTRAAODTRACKCUTS_H
+
+/*  See cxx source for full Copyright notice */
+
+//-------------------------------------------------------------------------
+//                      AliSpectraAODTrackCuts
+//
+//
+//
+//
+// Authors: Michele Floris, CERN, Philip Versteeg, UU, Redmer Bertens, UU
+//-------------------------------------------------------------------------
+
+class AliAODEvent;
+class AliSpectraAODHistoManager;
+
+#include "TNamed.h"
+
+class AliSpectraAODTrackCuts : public TNamed
+{
+public:
+
+   enum { kTrkBit = 0, kTrkEta, kTrkDCA, kTrkP, kTrkPt, kNTrkCuts};
+
+
+   AliSpectraAODTrackCuts() : TNamed(), fIsSelected(0), fTrackBits(0), fEtaCut(0), fPCut(0), fPtCut(0), fHistoCuts(0), fTrack(0) {}
+
+   AliSpectraAODTrackCuts(const char *name);
+   virtual  ~AliSpectraAODTrackCuts() {} // To be implemented
+
+   Bool_t IsSelected(AliAODTrack * track);
+
+   void SetTrackType(UInt_t bit);
+   Bool_t CheckTrackType();
+   Bool_t CheckEtaCut();
+   Bool_t CheckDCACut();
+   Bool_t CheckPCut();
+   Bool_t CheckPtCut();
+   void PrintCuts() const;
+
+   UInt_t GetTrackType()  const    { return fTrackBits; }
+   TH1I * GetHistoCuts()      { return fHistoCuts; }
+   void SetEta(Float_t eta)   { fEtaCut = eta; }
+   void SetDCA(Float_t dca)   { fDCACut = dca; }
+   void SetP(Float_t p)       { fPCut = p; }
+   void SetPt(Float_t pt)     { fPtCut = pt; }
+   Float_t GetEta()       const    { return fEtaCut; }
+   Float_t GetDCA()       const    { return fDCACut; }
+   Float_t GetP()         const    { return fPCut; }
+   Float_t GetPt()        const    { return fPtCut; }
+
+  Long64_t Merge(TCollection* list);
+
+
+private:
+
+   Bool_t         fIsSelected;      // True if cuts are selected
+   UInt_t         fTrackBits;       // Type of track to be used
+   Float_t        fEtaCut;          // Allowed absolute maximum value of Eta
+   Float_t        fDCACut;          // Maximum value of DCA
+   Float_t        fPCut;            // Maximum value of P
+   Float_t        fPtCut;           // Maximum value of Pt
+
+   TH1I *         fHistoCuts;       // Cuts statistics
+   AliAODTrack *  fTrack;           //! Track pointer
+   
+   AliSpectraAODTrackCuts(const AliSpectraAODTrackCuts&);
+   AliSpectraAODTrackCuts& operator=(const AliSpectraAODTrackCuts&);
+
+   ClassDef(AliSpectraAODTrackCuts, 1);
+};
+#endif
+
diff --git a/PWG2/SPECTRA/SpectraAOD/README b/PWG2/SPECTRA/SpectraAOD/README
new file mode 100644 (file)
index 0000000..91849fc
--- /dev/null
@@ -0,0 +1 @@
+This code was developed in Summer 2011 by 2 summer students (Philip Versteeg, UU, Redmer Bertens, UU) under the supervision of Michele Floris. It can serve as a starting point for AOD analysis.
diff --git a/PWG2/SPECTRA/SpectraAOD/analysis_macro.C b/PWG2/SPECTRA/SpectraAOD/analysis_macro.C
new file mode 100644 (file)
index 0000000..438fbe2
--- /dev/null
@@ -0,0 +1,356 @@
+#include "AliAODTrack.h"
+#include <iostream>
+#include "TFile.h"
+#include "TH1I.h"
+#include "TH2.h"
+#include "AliSpectraAODHistoManager.h"
+#include "AliSpectraAODEventCuts.h"
+#include "AliSpectraAODTrackCuts.h"
+#include "TCanvas.h"
+#include "TLegend.h"
+
+using namespace std;
+
+void analysis_macro(const char * dataFile = "Pt.AOD.1._data_ptcut.root", const char * mcFile ="Pt.AOD.1._MC.root") {
+   
+  // gSystem->Load("libTree.so");
+  // gSystem->Load("libGeom.so");
+  // gSystem->Load("libVMC.so");
+  // gSystem->Load("libPhysics.so");
+  // gSystem->Load("libSTEERBase.so");
+  // gSystem->Load("libESD.so");
+  // gSystem->Load("libAOD.so");
+  // gSystem->Load("libANALYSIS.so");
+  // gSystem->Load("libANALYSISalice.so");
+  // gSystem->Load("libANALYSIS");
+  // gSystem->Load("libANALYSISalice");
+  // gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
+  // gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+  // gStyle->SetPalette(1);
+  // gStyle->SetFillColor(kWhite);
+
+  // gROOT->LoadMacro("AliSpectraAODTrackCuts.cxx+g");
+  // gROOT->LoadMacro("AliSpectraAODEventCuts.cxx+g");
+  // gROOT->LoadMacro("AliSpectraAODHistoManager.cxx+g");
+  // gROOT->LoadMacro("AliAnalysisTaskSpectraAOD.cxx+g");
+  const char * histoname[] =
+    {
+      // names of histos we want to draw. please observe order of histos, elsewise name will not make sense
+      "#frac{MC_{#sigma, rec}}{Data}, P^{+}",
+      "#frac{MC_{#sigma, rec}}{Data}, K^{+}",
+      "#frac{MC_{#sigma, rec}}{Data}, #pi^{+}",
+      "#frac{MC_{#sigma, rec}}{Data}, P^{-}",
+      "#frac{MC_{#sigma, rec}}{Data}, K^{-}",
+      "#frac{MC_{#sigma, rec}}{Data}, #pi^{-}",      
+          
+      "#frac{MC_{#sigma, rec, prim}}{Data}, P^{+}",
+      "#frac{MC_{#sigma, rec, prim}}{Data}, K^{+}",
+      "#frac{MC_{#sigma, rec, prim}}{Data}, #pi^{+}",
+      "#frac{MC_{#sigma, rec, prim}}{Data}, P^{-}",
+      "#frac{MC_{#sigma, rec, prim}}{Data}, K^{-}",
+      "#frac{MC_{#sigma, rec, prim}}{Data}, #pi^{-}",     
+         
+      "#frac{MC_{#sigma, rec, sec}}{Data}, P^{+}",
+      "#frac{MC_{#sigma, rec, sec}}{Data}, K^{+}",
+      "#frac{MC_{#sigma, rec, sec}}{Data}, #pi^{+}",
+      "#frac{MC_{#sigma, rec, sec}}{Data}, P^{-}",
+      "#frac{MC_{#sigma, rec, sec}}{Data}, K^{-}",
+      "#frac{MC_{#sigma, rec, sec}}{Data}, #pi^{-}",     
+         
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, P^{+}",
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, K^{+}",
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, #pi^{+}",
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, P^{-}",
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, K^{-}",
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}}, #pi^{-}",     
+              
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, P^{+}",
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, K^{+}",
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, #pi^{+}",
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, P^{-}",
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, K^{-}",
+      "#frac{MC_{true, gen}}{MC_{#sigma, rec}} x Data, #pi^{-}",     
+
+      "Data, P^{+}",
+      "Data, K^{+}",
+      "Data, #pi^{+}",
+      "Data, P^{-}",
+      "Data, K^{-}",
+      "Data, #pi^{-}",      
+
+      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, P^{+}",
+      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, K^{+}",
+      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, #pi^{+}",
+      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, P^{-}",
+      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, K^{-}",
+      "#frac{MC_{rec #sigma}}{MC_{true, rec}}, #pi^{-}",     
+
+      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, P^{+}",
+      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, K^{+}",
+      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, #pi^{+}",
+      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, P^{-}",
+      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, K^{-}",
+      "#frac{MC_{rec #sigma sec}}{MC_{#sigma, rec}}, #pi^{-}",     
+    };
+
+  const char * histotitle[] =
+    {
+      // titles of histos we want to draw. please observe order of histos, elsewise name will not make sense
+      "MC #sigma, rec / Data, P^{+}",
+      "MC #sigma, rec / Data, K^{+}",
+      "MC #sigma, rec / Data, #pi^{+}",
+      "MC #sigma, rec / Data, P^{-}",
+      "MC #sigma, rec / Data, K^{-}",
+      "MC #sigma, rec / Data, #pi^{-}",     
+   
+      "MC #sigma, rec, prim / Data, P^{+}",
+      "MC #sigma, rec, prim  / Data, K^{+}",
+      "MC #sigma, rec, prim  / Data, #pi^{+}",
+      "MC #sigma, rec, prim  / Data, P^{-}",
+      "MC #sigma, rec, prim  / Data, K^{-}",
+      "MC #sigma, rec, prim  / Data, #pi^{-}",      
+    
+      "MC #sigma, rec, sec / Data, P^{+}",
+      "MC #sigma, rec, sec / Data, K^{+}",
+      "MC #sigma, rec, sec / Data, #pi^{+}",
+      "MC #sigma, rec, sec / Data, P^{-}",
+      "MC #sigma, rec, sec / Data, K^{-}",
+      "MC #sigma, rec, sec / Data, #pi^{-}",      
+          
+      "Correction factor, P^{+}",
+      "Correction factor, K^{+}",
+      "Correction factor, #pi^{+}",
+      "Correction factor, P^{-}",
+      "Correction factor, K^{-}",
+      "Correction factor, #pi^{-}",      
+
+      "True Data, P^{+}",
+      "True Data, K^{+}",
+      "True Data, #pi^{+}",
+      "True Data, P^{-}",
+      "True Data, K^{-}",
+      "True Data, #pi^{-}",      
+
+      "Raw Data, P^{+}",
+      "Raw Data, K^{+}",
+      "Raw Data, #pi^{+}",
+      "Raw Data, P^{-}",
+      "Raw Data, K^{-}",
+      "Raw Data, #pi^{-}",      
+
+      "Identificaiton quality, P^{+}",
+      "Identificaiton quality, K^{+}",
+      "Identificaiton quality, #pi^{+}",
+      "Identificaiton quality, P^{-}",
+      "Identificaiton quality, K^{-}",
+      "Identificaiton quality, #pi^{-}",     
+
+      "Contribution of secondaries, P^{+}",
+      "Contribution of secondaries, K^{+}",
+      "Contribution of secondaries, #pi^{+}",
+      "Contribution of secondaries, P^{-}",
+      "Contribution of secondaries, K^{-}",
+      "Contribution of secondaries, #pi^{-}",     
+    };
+
+  // Open root MC file and get classes
+  cout << "Analysis Macro" << endl;
+  cout << "  > Reading MC data" << endl;
+  TFile *_mc = TFile::Open(mcFile);
+  AliSpectraAODHistoManager* hman_mc = (AliSpectraAODHistoManager*) _mc->Get("SpectraHistos");
+  AliSpectraAODEventCuts* ecuts_mc = (AliSpectraAODEventCuts*) _mc->Get("Event Cuts");
+  AliSpectraAODTrackCuts* tcuts_mc = (AliSpectraAODTrackCuts*) _mc->Get("Track Cuts");
+  // print info about mc track and Event cuts
+  cout << " -- Info about MC -- "<< endl;
+  ecuts_mc->PrintCuts();
+  tcuts_mc->PrintCuts();
+  // get the mc histo's necessary to write the contamination histos from the rootfile
+  TH1F* spectrahistos_mc[30];
+  cout << " -- Reading and normalizing histograms -- " << endl;
+  for( Int_t i = 0 ; i <= AliSpectraNameSpace::kNPtSpecies ; i++ ) { 
+    spectrahistos_mc[i] = dynamic_cast<TH1F*>(hman_mc->GetHistogram((AliSpectraNameSpace::AODPtHist_t)i)); 
+    Double_t events_mc =  1. / ecuts_mc->NumberOfEvents();
+    if(!spectrahistos_mc[i]) {
+      continue;
+    }    
+    spectrahistos_mc[i]->Scale(events_mc, "width");
+  }
+  cout << "  > Reading real data " << endl;
+  // proceed likewise for data
+  TFile *_data = TFile::Open(dataFile);
+  AliSpectraAODHistoManager* hman_data = (AliSpectraAODHistoManager*) _data->Get("SpectraHistos");
+  AliSpectraAODEventCuts* ecuts_data = (AliSpectraAODEventCuts*) _data->Get("Event Cuts");
+  AliSpectraAODTrackCuts* tcuts_data = (AliSpectraAODTrackCuts*) _data->Get("Track Cuts");
+  // print info about track and Event cuts
+  cout << " -- Info about data -- " << endl;
+  ecuts_data->PrintCuts();
+  tcuts_data->PrintCuts();
+  // get the histo's necessary to write the contamination histos from the rootfile
+  TH1F* spectrahistos_data[30];
+  cout << " -- Reading and normalizing histograms -- " << endl;
+  for( Int_t i = 0 ; i <= AliSpectraNameSpace::kNPtSpecies ; i++ ) { 
+    spectrahistos_data[i] = dynamic_cast<TH1F*>(hman_data->GetHistogram((AliSpectraNameSpace::AODPtHist_t)i));
+    if(!spectrahistos_data[i]) {
+      continue;
+    }    
+    Double_t events_data =  1. / ecuts_data->NumberOfEvents();
+    //normalize the histos    
+    spectrahistos_data[i]->Scale(events_data, "width");
+  }
+  //create output file and select it as the active file
+  TFile * output = new TFile("analysis_output.root", "RECREATE");
+  output->cd();
+  // Write the scaled data and MC histos in the output file for convenience
+  output->mkdir("Data");
+  output->cd("Data");
+  for(Int_t ihist = 0; ihist <= AliSpectraNameSpace::kNPtSpecies; ihist++){
+      spectrahistos_data[ihist]->Write();
+  }
+  output->mkdir("MC");
+  output->cd("MC");
+  for(Int_t ihist = 0; ihist <= AliSpectraNameSpace::kNPtSpecies; ihist++){
+      spectrahistos_mc[ihist]->Write();
+  }
+  output->mkdir("ANALYSIS");
+  output->cd("ANALYSIS");
+  // perform some analysis!
+  //  next loop: MC_Sigma_Rec over Rec_Data
+  TCanvas* c1 = new TCanvas("MC reconstructed over Data","MC reconstructed over Data");
+  c1->Divide(3,2);
+  TH1F* kRatioPtMCSigmaRecSigmaRec[6];
+  for (Int_t i = 0 ; i < 6 ; i ++ )
+    {
+      kRatioPtMCSigmaRecSigmaRec[i] = (TH1F*)(spectrahistos_mc[i])->Clone(histotitle[i]);
+      kRatioPtMCSigmaRecSigmaRec[i]->Divide(spectrahistos_data[i]);
+      c1->cd(1+i);
+      kRatioPtMCSigmaRecSigmaRec[i]->SetTitle(histoname[i]);
+      kRatioPtMCSigmaRecSigmaRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
+      kRatioPtMCSigmaRecSigmaRec[i]->Draw();
+      kRatioPtMCSigmaRecSigmaRec[i]->Write();
+    }
+  // next loop: MC_sigma_primaries over Rec_data
+  TCanvas* c2 = new TCanvas("MC reconstructed primaries over Data", "MC reconstructed primaries over Data");
+  c2->Divide(3,2);
+  TH1F* kRatioPtMCSigmaRecPrimarySigmaRec[6];
+  for (Int_t i = 0 ; i < 6 ; i ++ )
+    {
+      kRatioPtMCSigmaRecPrimarySigmaRec[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaPrimaryProtonPlus])->Clone(histotitle[i+6]);
+      kRatioPtMCSigmaRecPrimarySigmaRec[i]->Divide(spectrahistos_data[i]);
+      c2->cd(1+i);
+      kRatioPtMCSigmaRecPrimarySigmaRec[i]->SetTitle(histoname[i+6]);
+      kRatioPtMCSigmaRecPrimarySigmaRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
+      kRatioPtMCSigmaRecPrimarySigmaRec[i]->Draw();
+      kRatioPtMCSigmaRecPrimarySigmaRec[i]->Write();
+    }
+  // next loop: MC_Sigma_secondaries over Rec_data
+  TCanvas* c3 = new TCanvas("MC reconstructed secondaries over Data","MC reconstructed secondaries over Data");
+  c3->Divide(3,2);
+  TH1F* kRatioPtMCSigmaRecSigmaRecSecondary[6];
+  for (Int_t i = 0 ; i < 6 ; i ++ )
+    {
+      kRatioPtMCSigmaRecSigmaRecSecondary[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+12]);
+      kRatioPtMCSigmaRecSigmaRecSecondary[i]->Divide(spectrahistos_data[i]);
+      c3->cd(1+i);
+      kRatioPtMCSigmaRecSigmaRecSecondary[i]->SetTitle(histoname[12+i]);
+      kRatioPtMCSigmaRecSigmaRecSecondary[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
+      kRatioPtMCSigmaRecSigmaRecSecondary[i]->Draw();
+      kRatioPtMCSigmaRecSigmaRecSecondary[i]->Write();
+    }
+  //
+  TCanvas* c4 = new TCanvas("MC true generated over MC true reconstructed","MC true generated over MC true reconstructed");
+  c4->Divide(3,2);
+  TH1F* kRatioPtMCTrueMCRec[6];
+  for (Int_t i = 0 ; i < 6 ; i ++ )
+    {
+      kRatioPtMCTrueMCRec[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus])->Clone(histotitle[i+18]);
+      kRatioPtMCTrueMCRec[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaPrimaryProtonPlus]);
+      c4->cd(1+i);
+      kRatioPtMCTrueMCRec[i]->SetTitle(histoname[18+i]);
+      kRatioPtMCTrueMCRec[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
+      kRatioPtMCTrueMCRec[i]->Draw();
+      kRatioPtMCTrueMCRec[i]->Write();
+    }
+  
+  TCanvas* c5 = new TCanvas("True Data","True Data");
+  c5->Divide(3,2);
+  TH1F* kRatioPtTrueData[6];
+  for (Int_t i = 0 ; i < 6 ; i ++ )
+    {
+      kRatioPtTrueData[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus])->Clone(histotitle[i+24]);
+      kRatioPtTrueData[i]->Divide(spectrahistos_mc[i]); // Michele: Since we are correcting all contributions at once, this should be exaclty the same as in the data!
+      kRatioPtTrueData[i]->Multiply(spectrahistos_data[i]);
+      c5->cd(1+i);
+      kRatioPtTrueData[i]->SetTitle(histoname[24+i]);
+      kRatioPtTrueData[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
+      kRatioPtTrueData[i]->Draw();
+      kRatioPtTrueData[i]->Write();
+    }
+  TCanvas* c6 = new TCanvas("Raw data", "Raw data");
+  c6->Divide(3,2);
+  TH1F* kRawData[6];
+  for (Int_t i = 0 ; i < 6 ; i ++ )
+    {
+      kRawData[i] = (TH1F*)(spectrahistos_data[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus])->Clone(histotitle[i+30]);
+      c6->cd(1+i);
+      kRawData[i]->SetTitle(histoname[30+i]);
+      kRawData[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
+      kRawData[i]->Draw();
+      kRawData[i]->Write();
+    }
+  
+  TCanvas* c7 = new TCanvas("contamination", "contamination");
+  c7->Divide(3,2);
+  TH1F* kContamination[6];
+  for (Int_t i = 0 ; i < 6 ; i ++ )
+    {
+      kContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus])->Clone(histotitle[i+36]);
+      kContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecTrueProtonPlus]);
+      c7->cd(1+i);
+      kContamination[i]->SetTitle(histoname[36+i]);
+      kContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
+      kContamination[i]->Draw();
+      kContamination[i]->Write();
+    }
+  
+  TCanvas* c8 = new TCanvas("contamination secondary", "contamination secondary");
+  c8->Divide(3,2);
+  TH1F* kSecContamination[6];
+  for (Int_t i = 0 ; i < 6 ; i ++ )
+    {
+      kSecContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+42]);
+      kSecContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus]);
+      c8->cd(1+i);
+      kSecContamination[i]->SetTitle(histoname[42+i]);
+      kSecContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
+      kSecContamination[i]->Draw();
+      kSecContamination[i]->Write();
+    }
+   
+  TCanvas* c9 = new TCanvas("true data and generated", "true data and generated");
+  c9->Divide(3,2);
+  TH1F* ktrueandgenerated[6];
+  for (Int_t i = 0 ; i < 6 ; i ++ )
+    {
+      // kSecContamination[i] = (TH1F*)(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaSecondaryProtonPlus])->Clone(histotitle[i+42]);
+      // kSecContamination[i]->Divide(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtRecSigmaProtonPlus]);
+      c9->cd(1+i);
+      // kSecContamination[i]->SetTitle(histoname[42+i]);
+      // kSecContamination[i]->GetYaxis()->SetTitle("#frac{d^{2} N}{dy dp_{T}} (c / GeV)");
+      kRatioPtTrueData[i]->Draw();
+      spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetMarkerColor(kRed);
+      spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetMarkerStyle(kOpenCircle);
+      spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->SetLineColor(kRed);
+      spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus]->Draw("same");
+
+      // legend
+      TLegend * leg = new TLegend(0.6,0.64,0.89,0.79);  //coordinates are fractions
+      leg->AddEntry(spectrahistos_mc[i+AliSpectraNameSpace::kHistPtGenTruePrimaryProtonPlus],"Generated","p"); 
+      leg->AddEntry(kRatioPtTrueData[i],"Corrected","p"); 
+      leg->Draw();
+    }
+   
+  cout << "Analysis complete" << endl;
+  cout << " - Please exit root to close and save output file " << endl;
+}
diff --git a/PWG2/SPECTRA/SpectraAOD/run.C b/PWG2/SPECTRA/SpectraAOD/run.C
new file mode 100644 (file)
index 0000000..4b48f85
--- /dev/null
@@ -0,0 +1,32 @@
+void run()
+{
+
+   gSystem->Load("libTree.so");
+   gSystem->Load("libGeom.so");
+   gSystem->Load("libVMC.so");
+   gSystem->Load("libPhysics.so");
+   gSystem->Load("libSTEERBase.so");
+   gSystem->Load("libESD.so");
+   gSystem->Load("libAOD.so");
+   gSystem->Load("libANALYSIS.so");
+   gSystem->Load("libANALYSISalice.so");
+   gSystem->Load("libANALYSIS");
+   gSystem->Load("libANALYSISalice");
+   gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
+   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+   gStyle->SetPalette(1);
+   gStyle->SetFillColor(kWhite);
+
+   gROOT->LoadMacro("AliSpectraAODTrackCuts.cxx+g");
+   gROOT->LoadMacro("AliSpectraAODEventCuts.cxx+g");
+   gROOT->LoadMacro("AliSpectraAODHistoManager.cxx+g");
+   gROOT->LoadMacro("AliAnalysisTaskSpectraAOD.cxx+g");
+
+
+
+
+
+
+
+
+}
diff --git a/PWG2/SPECTRA/SpectraAOD/runAOD.C b/PWG2/SPECTRA/SpectraAOD/runAOD.C
new file mode 100644 (file)
index 0000000..8815a9e
--- /dev/null
@@ -0,0 +1,78 @@
+void runAOD(Bool_t isMC = 0)
+{
+
+//gEnv->SetValue("XSec.GSI.DelegProxy", "2");
+//TProof::Open("pverstee@alice-caf.cern.ch");
+
+   gSystem->Load("libTree.so");
+   gSystem->Load("libGeom.so");
+   gSystem->Load("libVMC.so");
+   gSystem->Load("libPhysics.so");
+   gSystem->Load("libSTEERBase.so");
+   gSystem->Load("libESD.so");
+   gSystem->Load("libAOD.so");
+   gSystem->Load("libANALYSIS.so");
+   gSystem->Load("libANALYSISalice.so");
+   gSystem->Load("libANALYSIS");
+   gSystem->Load("libANALYSISalice");
+   gROOT->LoadMacro("$ALICE_ROOT/PWG0/CreateESDChain.C");
+
+   ifstream flist;
+   flist.open("filelist.txt");
+   TChain * chain = new TChain("aodTree");
+   TString line;
+   while (line.ReadLine(flist))
+   {
+      gSystem->ExpandPathName(line);
+      chain->AddFile(line.Data());
+      cout << "Adding file " << line.Data() << endl;
+   }
+//  gProof->UploadPackage("AF-v4-19-04-AN");
+//  gProof->EnablePackage("AF-v4-19-04-AN");
+   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+   gStyle->SetPalette(1);
+
+
+
+   AliAnalysisManager *mgr = new AliAnalysisManager("testAnalysis");
+   AliAODInputHandler* aodH = new AliAODInputHandler();
+   mgr->SetInputEventHandler(aodH);
+
+   gROOT->LoadMacro("AliSpectraAODTrackCuts.cxx+g");
+   gROOT->LoadMacro("AliSpectraAODVertexCuts.cxx+g");
+   gROOT->LoadMacro("AliSpectraAODHistoManager.cxx+g");
+   gROOT->LoadMacro("AliAnalysisTaskSpectraAOD.cxx+g");
+
+   gROOT->LoadMacro("$(ALICE_ROOT)/ANALYSIS/macros/AddTaskPIDResponse.C");
+   AliAnalysisTask * taskPID = AddTaskPIDResponse();
+   mgr->AddTask(taskPID);
+
+
+   AliAnalysisTaskSpectraAOD *task = new AliAnalysisTaskSpectraAOD("TaskAODExercise");
+   task->SetIsMC(isMC);
+   mgr->AddTask(task);
+
+   // Set the cuts
+   AliSpectraAODVertexCuts * vcuts = new AliSpectraAODVertexCuts("Vertex Cuts");
+   AliSpectraAODTrackCuts  * tcuts = new AliSpectraAODTrackCuts("Tracks Cuts");
+   tcuts->SetTrackType(6);
+   tcuts->SetEta(1.);
+   task->SetVertexCuts(vcuts);
+   task->SetTrackCuts(tcuts);
+
+
+   AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+   AliAnalysisDataContainer *coutputpt1 = mgr->CreateContainer("chistpt", AliSpectraAODHistoManager::Class(),  AliAnalysisManager::kOutputContainer, "Pt.AOD.1.root");
+   AliAnalysisDataContainer *coutputpt2 = mgr->CreateContainer("cvcutpt", AliSpectraAODVertexCuts::Class(),    AliAnalysisManager::kOutputContainer, "Pt.AOD.1.root");
+   AliAnalysisDataContainer *coutputpt3 = mgr->CreateContainer("ctcutpt", AliSpectraAODTrackCuts::Class(),     AliAnalysisManager::kOutputContainer, "Pt.AOD.1.root");
+
+   mgr->ConnectInput(task, 0, cinput);
+   mgr->ConnectOutput(task, 1, coutputpt1);
+   mgr->ConnectOutput(task, 2, coutputpt2);
+   mgr->ConnectOutput(task, 3, coutputpt3);
+   mgr->SetDebugLevel(2);
+
+   if (!mgr->InitAnalysis()) return;
+   mgr->PrintStatus();
+   mgr->StartAnalysis("local", chain, 100);
+}
diff --git a/PWG2/SPECTRA/SpectraAOD/runAODProof.C b/PWG2/SPECTRA/SpectraAOD/runAODProof.C
new file mode 100644 (file)
index 0000000..acd4f8c
--- /dev/null
@@ -0,0 +1,84 @@
+void runAODProof(const char * proofMode = "full")
+{
+
+   gEnv->SetValue("XSec.GSI.DelegProxy", "2");
+   //  TProof::Open("rbertens@alice-caf.cern.ch");
+
+   gSystem->Load("libTree.so");
+   gSystem->Load("libGeom.so");
+   gSystem->Load("libVMC.so");
+   gSystem->Load("libPhysics.so");
+   gSystem->Load("libSTEERBase.so");
+   gSystem->Load("libESD.so");
+   gSystem->Load("libAOD.so");
+   gSystem->Load("libANALYSIS.so");
+   gSystem->Load("libOADB.so");
+   gSystem->Load("libANALYSISalice.so");
+   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+
+
+   AliAnalysisAlien * handler = new AliAnalysisAlien("test");
+   handler->SetOverwriteMode();
+   handler->SetRunMode(proofMode);
+   handler->SetProofReset(0);
+   handler->SetAliROOTVersion("v4-21-29-AN");
+   handler->SetProofCluster(Form("%s@alice-caf.cern.ch", gSystem->Getenv("CAFUSER")));
+   //handler->SetProofCluster("rbertens@skaf.saske.sk");
+   //   handler->SetProofDataSet("/alice/data/LHC10h_000138653_p2_AOD049#aodTree");
+   handler->SetProofDataSet("/alice/sim/LHC11a10a_000138653_AOD048#aodTree");
+   handler->SetNproofWorkersPerSlave(1);
+   handler->SetAliRootMode("default");
+   handler->SetAdditionalLibs("AliSpectraAODHistoManager.cxx AliSpectraAODHistoManager.h AliSpectraAODEventCuts.cxx AliSpectraAODEventCuts.h AliSpectraAODTrackCuts.cxx AliSpectraAODTrackCuts.h AliAnalysisTaskSpectraAOD.cxx AliAnalysisTaskSpectraAOD.h");
+   handler->SetAnalysisSource("AliSpectraAODHistoManager.cxx+ AliSpectraAODEventCuts.cxx+ AliSpectraAODTrackCuts.cxx+ AliAnalysisTaskSpectraAOD.cxx+");
+   handler->SetFileForTestMode("filelist.txt"); // list of local files for testing
+
+   //  handler->SetAliRootMode("");
+   handler->SetClearPackages();
+
+
+   AliAnalysisManager *mgr = new AliAnalysisManager("testAnalysis");
+   mgr->SetGridHandler(handler);
+   AliAODInputHandler* aodH = new AliAODInputHandler();
+   mgr->SetInputEventHandler(aodH);
+
+   gROOT->LoadMacro("AliSpectraAODTrackCuts.cxx+g");
+   gROOT->LoadMacro("AliSpectraAODEventCuts.cxx+g");
+   gROOT->LoadMacro("AliSpectraAODHistoManager.cxx+g");
+   gROOT->LoadMacro("AliAnalysisTaskSpectraAOD.cxx+g");
+
+   // Add PID task
+   gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
+   AliAnalysisTask * taskPID = AddTaskPIDResponse();
+   mgr->AddTask(taskPID);
+
+   AliAnalysisTaskSpectraAOD *task = new AliAnalysisTaskSpectraAOD("TaskAODExercise");
+   task->SetIsMC(1);
+   mgr->AddTask(task);
+
+   // Set the cuts
+   AliSpectraAODEventCuts * vcuts = new AliSpectraAODEventCuts("Event Cuts");
+   AliSpectraAODTrackCuts  * tcuts = new AliSpectraAODTrackCuts("Track Cuts");
+   tcuts->SetTrackType(6);
+   tcuts->SetEta(1.);
+   tcuts->SetP(.7);
+   // vcuts->SetCentralityCutMin(0.0) // default
+   // vcuts->SetCentralityCutMax(0.0) // default
+   task->SetEventCuts(vcuts);
+   task->SetTrackCuts(tcuts);
+   task->SetNSigmaForIdentification(3.);
+
+   AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+   AliAnalysisDataContainer *coutputpt1 = mgr->CreateContainer("chistpt", AliSpectraAODHistoManager::Class(),  AliAnalysisManager::kOutputContainer, "Pt.AOD.1.root");
+   AliAnalysisDataContainer *coutputpt2 = mgr->CreateContainer("cvcutpt", AliSpectraAODEventCuts::Class(),    AliAnalysisManager::kOutputContainer, "Pt.AOD.1.root");
+   AliAnalysisDataContainer *coutputpt3 = mgr->CreateContainer("ctcutpt", AliSpectraAODTrackCuts::Class(),     AliAnalysisManager::kOutputContainer, "Pt.AOD.1.root");
+
+   mgr->ConnectInput(task, 0, cinput);
+   mgr->ConnectOutput(task, 1, coutputpt1);
+   mgr->ConnectOutput(task, 2, coutputpt2);
+   mgr->ConnectOutput(task, 3, coutputpt3);
+   mgr->SetDebugLevel(2);
+
+   if (!mgr->InitAnalysis()) return;
+   mgr->PrintStatus();
+   mgr->StartAnalysis("proof");
+}
diff --git a/PWG2/SPECTRA/SpectraAOD/runAODProof_init.C b/PWG2/SPECTRA/SpectraAOD/runAODProof_init.C
new file mode 100644 (file)
index 0000000..7f39992
--- /dev/null
@@ -0,0 +1,117 @@
+void runAODProof(Int_t c, const char * proofMode = "full")
+{
+
+   gEnv->SetValue("XSec.GSI.DelegProxy", "2");
+
+   gSystem->Load("libTree.so");
+   gSystem->Load("libGeom.so");
+   gSystem->Load("libVMC.so");
+   gSystem->Load("libPhysics.so");
+   gSystem->Load("libSTEERBase.so");
+   gSystem->Load("libESD.so");
+   gSystem->Load("libAOD.so");
+   gSystem->Load("libANALYSIS.so");
+   gSystem->Load("libOADB.so");
+   gSystem->Load("libANALYSISalice.so");
+   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+
+
+   AliAnalysisAlien * handler = new AliAnalysisAlien("test");
+   handler->SetOverwriteMode();
+   handler->SetRunMode(proofMode);
+   handler->SetProofReset(0);
+   handler->SetAliROOTVersion("v4-21-29-AN");
+
+   handler->SetProofCluster(Form("%s@alice-caf.cern.ch", gSystem->Getenv("CAFUSER")));
+//   handler->SetProofCluster(Form("%s@skaf.saske.sk",gSystem->Getenv("CAFUSER")));
+
+   // Set handler for Real DATA:
+   if (c == 1)
+   {
+     handler->SetProofDataSet("/alice/data/LHC10h_000138653_p2_AOD049#aodTree");
+ //     handler->SetProofDataSet("/alice/sim/LHC11a10a_000139107_AOD048#aodTree|alice/sim/LHC11a10a_000138653_AOD048#aodTree");
+   }
+   if (c == 2)
+   {
+      handler->SetProofDataSet("/alice/sim/LHC11a10a_000139107_AOD048#aodTree|/alice/sim/LHC11a10a_000138653_AOD048#aodTree|/alice/sim/LHC11a10a_000139110_AOD048#aodTree|/alice/sim/LHC11a10a_000138662_AOD048#aodTree|/alice/sim/LHC11a10a_000138666_AOD048#aodTree|/alice/sim/LHC11a10a_000138795_AOD048#aodTree");      
+   }
+   handler->SetNproofWorkersPerSlave(1);
+   handler->SetAliRootMode("default");
+   handler->SetAdditionalLibs("AliSpectraAODHistoManager.cxx AliSpectraAODHistoManager.h AliSpectraAODEventCuts.cxx AliSpectraAODEventCuts.h AliSpectraAODTrackCuts.cxx AliSpectraAODTrackCuts.h AliAnalysisTaskSpectraAOD.cxx AliAnalysisTaskSpectraAOD.h");
+   handler->SetAnalysisSource("AliSpectraAODHistoManager.cxx+ AliSpectraAODEventCuts.cxx+ AliSpectraAODTrackCuts.cxx+ AliAnalysisTaskSpectraAOD.cxx+");
+   handler->SetFileForTestMode("filelist.txt"); // list of local files for testing
+   //  handler->SetAliRootMode("");
+   handler->SetClearPackages();
+
+
+   AliAnalysisManager *mgr = new AliAnalysisManager("testAnalysis");
+   mgr->SetGridHandler(handler);
+   AliAODInputHandler* aodH = new AliAODInputHandler();
+   mgr->SetInputEventHandler(aodH);
+
+   gROOT->LoadMacro("AliSpectraAODTrackCuts.cxx+g");
+   gROOT->LoadMacro("AliSpectraAODEventCuts.cxx+g");
+   gROOT->LoadMacro("AliSpectraAODHistoManager.cxx+g");
+   gROOT->LoadMacro("AliAnalysisTaskSpectraAOD.cxx+g");
+
+
+   // Add PID task
+   gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
+   Bool_t isMC = kFALSE;
+   if (c == 2 || c == 3) isMC = kTRUE;   
+   AliAnalysisTask * taskPID = AddTaskPIDResponse(isMC);
+   mgr->AddTask(taskPID);
+
+
+   AliAnalysisTaskSpectraAOD *task = new AliAnalysisTaskSpectraAOD("TaskAODExercise");
+   mgr->AddTask(task);
+
+   // Set the cuts
+   AliSpectraAODEventCuts * vcuts = new AliSpectraAODEventCuts("Event Cuts");
+   AliSpectraAODTrackCuts  * tcuts = new AliSpectraAODTrackCuts("Track Cuts");
+   tcuts->SetTrackType(6);
+   tcuts->SetEta(.8);
+   //   tcuts->SetDCA(.1);
+   tcuts->SetPt(1.2);   
+   vcuts->SetCentralityCutMax(5.);  // example min max cuts
+   vcuts->SetCentralityCutMin(0.);
+   task->SetEventCuts(vcuts);
+   task->SetTrackCuts(tcuts);
+   task->SetNSigmaForIdentification(3.); // FIXME
+   task->SetYCut(.5);
+   // check for MC or real data
+   if (c == 2 || c == 3)
+   {
+      task->SetIsMC(kTRUE);
+
+
+      AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+      AliAnalysisDataContainer *coutputpt1 = mgr->CreateContainer("chistpt", AliSpectraAODHistoManager::Class(),  AliAnalysisManager::kOutputContainer, "Pt.AOD.1._MC.root");
+      AliAnalysisDataContainer *coutputpt2 = mgr->CreateContainer("cvcutpt", AliSpectraAODEventCuts::Class(),    AliAnalysisManager::kOutputContainer, "Pt.AOD.1._MC.root");
+      AliAnalysisDataContainer *coutputpt3 = mgr->CreateContainer("ctcutpt", AliSpectraAODTrackCuts::Class(),     AliAnalysisManager::kOutputContainer, "Pt.AOD.1._MC.root");
+   }
+   if (c == 1)
+   {
+      AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
+      AliAnalysisDataContainer *coutputpt1 = mgr->CreateContainer("chistpt", AliSpectraAODHistoManager::Class(),  AliAnalysisManager::kOutputContainer, "Pt.AOD.1._data_ptcut.root");
+      AliAnalysisDataContainer *coutputpt2 = mgr->CreateContainer("cvcutpt", AliSpectraAODEventCuts::Class(),    AliAnalysisManager::kOutputContainer, "Pt.AOD.1._data_ptcut.root");
+      AliAnalysisDataContainer *coutputpt3 = mgr->CreateContainer("ctcutpt", AliSpectraAODTrackCuts::Class(),     AliAnalysisManager::kOutputContainer, "Pt.AOD.1._data_ptcut.root");
+
+   }
+   mgr->ConnectInput(task, 0, cinput);
+   mgr->ConnectOutput(task, 1, coutputpt1);
+   mgr->ConnectOutput(task, 2, coutputpt2);
+   mgr->ConnectOutput(task, 3, coutputpt3);
+   mgr->SetDebugLevel(2);
+
+   if (!mgr->InitAnalysis()) return;
+   mgr->PrintStatus();
+   if (c == 3)
+   {
+      mgr->StartAnalysis("local");
+   }
+   if (c != 3)
+   {
+      mgr->StartAnalysis("proof");
+   }
+}
diff --git a/PWG2/SPECTRA/SpectraAOD/runProof.C b/PWG2/SPECTRA/SpectraAOD/runProof.C
new file mode 100644 (file)
index 0000000..a2f0e30
--- /dev/null
@@ -0,0 +1,261 @@
+// TODO:
+// 1. run with Many centrality bins at once
+#include <string.h>
+
+enum { kMyRunModeLocal = 0, kMyRunModeCAF, kMyRunModeGRID};
+
+TList * listToLoad = new TList();
+
+TChain * GetAnalysisChain(const char * incollection);
+
+void runProof(Char_t* data = "/alice/data/LHC10h_000138653_p2_AOD049", Long64_t nev = -1, Long64_t offset = 0, Bool_t debug = kFALSE, Int_t runMode = 1, Bool_t isMC = 0,
+              const char* option = "SAVE", Int_t workers = -1)
+{
+   // runMode:
+   //
+   // 0 local
+   // 1 proof
+   // 2 grid
+
+   if (nev < 0)
+      nev = 1234567890;
+
+   InitAndLoadLibs(runMode, workers, debug);
+
+   // Create the analysis manager
+   mgr = new AliAnalysisManager;
+
+   // Add AOD handler
+   AliAODInputHandler* AODH = new AliAODInputHandler;
+   mgr->SetInputEventHandler(AODH);
+   mgr->SetDebugLevel(10);
+
+   // Is this needed for AOD?
+   // if(isMC) {
+   //   AliMCEventHandler* handler = new AliMCEventHandler;
+   //   handler->SetPreReadMode(AliMCEventHandler::kLmPreRead);
+   //   mgr->SetMCtruthEventHandler(handler);
+   // }
+
+
+   // If we are running on grid, we need the alien handler
+   if (runMode == kMyRunModeGRID)
+   {
+      // Create and configure the alien handler plugin
+      gROOT->LoadMacro("CreateAlienHandler.C");
+      AliAnalysisGrid *alienHandler = CreateAlienHandler(data, listToLoad, "full", isMC);
+      if (!alienHandler)
+      {
+         cout << "Cannot create alien handler" << endl;
+         exit(1);
+      }
+      mgr->SetGridHandler(alienHandler);
+   }
+
+   // Parse option strings
+   TString optionStr(option);
+
+   // remove SAVE option if set
+   // This  is copied from a macro by Jan. The reason I kept it is that I may want to pass textual options to the new task at some point
+   Bool_t doSave = kFALSE;
+   TString optionStr(option);
+   if (optionStr.Contains("SAVE"))
+   {
+      optionStr = optionStr(0, optionStr.Index("SAVE")) + optionStr(optionStr.Index("SAVE") + 4, optionStr.Length());
+      doSave = kTRUE;
+   }
+
+   //  AliLog::SetClassDebugLevel("AliESDtrackCuts", AliLog::kDebug);// FIXME
+
+   // load my task
+   gROOT->ProcessLine(".L AddTaskSpectraAOD.C");
+   AliAnalysisTaskSpectraAOD * task = AddTaskSpectraAOD("SpectraAOD.root");
+   task->SetIsMC(isMC);
+   // Init and run the analy
+   if (!mgr->InitAnalysis()) return;
+
+   mgr->PrintStatus();
+
+   if (runMode == kMyRunModeLocal)
+   {
+      // If running in local mode, create chain of ESD files
+      cout << "RUNNING LOCAL, CHAIN" << endl;
+      TChain * chain = GetAnalysisChain(data);
+      //    chain->Print();
+      mgr->StartAnalysis("local", chain, nev);
+   }
+   else if (runMode == kMyRunModeCAF)
+   {
+      mgr->StartAnalysis("proof", TString(data) + "#aodTree", nev);
+   }
+   else if (runMode == kMyRunModeGRID)
+   {
+      mgr->StartAnalysis("grid");
+   }
+   else
+   {
+      cout << "ERROR: unknown run mode" << endl;
+   }
+
+   TString pathsuffix = "";
+   if (doSave) MoveOutput(data, pathsuffix.Data());
+
+}
+
+
+void MoveOutput(const char * data, const char * suffix = "")
+{
+
+   TString path("output/");
+   path = path + TString(data).Tokenize("/")->Last()->GetName() + suffix;
+
+   TString fileName = "SpectraAOD.root";
+   gSystem->mkdir(path, kTRUE);
+   gSystem->Rename(fileName, path + "/" + fileName);
+   // for(Int_t ibin = 0; ibin < 20; ibin++){
+   //   TString fileBin = fileName;
+   //   fileBin.ReplaceAll(".root",Form("_%2.2d.root",ibin));
+   //   gSystem->Rename(fileBin, path + "/" + fileBin);
+   // }
+   Printf(">>>>> Moved files to %s", path.Data());
+}
+
+
+
+TChain * GetAnalysisChain(const char * incollection)
+{
+   // Builds a chain of esd files
+   // incollection can be
+   // - a single root file
+   // - an xml collection of files on alien
+   // - a ASCII containing a list of local root files
+   TChain* analysisChain = 0;
+   // chain
+   analysisChain = new TChain("aodTree");
+   if (TString(incollection).Contains(".root"))
+   {
+      analysisChain->Add(incollection);
+   }
+   else if (TString(incollection).Contains("xml"))
+   {
+      TGrid::Connect("alien://");
+      TAlienCollection * coll = TAlienCollection::Open(incollection);
+      while (coll->Next())
+      {
+         analysisChain->Add(TString("alien://") + coll->GetLFN());
+      }
+   }
+   else
+   {
+      ifstream file_collect(incollection);
+      TString line;
+      while (line.ReadLine(file_collect))
+      {
+         analysisChain->Add(line.Data());
+      }
+   }
+   analysisChain->GetListOfFiles()->Print();
+
+   return analysisChain;
+}
+
+
+void InitAndLoadLibs(Int_t runMode = kMyRunModeLocal, Int_t workers = 0, Bool_t debug = 0)
+{
+   // Loads libs and par files + custom task and classes
+
+   // Custom stuff to be loaded
+
+   listToLoad->Add(new TObjString("AliSpectraAODHistoManager.cxx+"));
+   listToLoad->Add(new TObjString("AliSpectraAODEventCuts.cxx+"));
+   listToLoad->Add(new TObjString("AliSpectraAODTrackCuts.cxx+"));
+   listToLoad->Add(new TObjString("AliAnalysisTaskSpectraAOD.cxx+"));
+
+   if (runMode == kMyRunModeCAF)
+   {
+      cout << "Init in CAF mode" << endl;
+
+      gEnv->SetValue("XSec.GSI.DelegProxy", "2");
+      TProof * p = TProof::Open("alice-caf.cern.ch", workers > 0 ? Form("workers=%d", workers) : "1x");
+      //TProof * p = TProof::Open("skaf.saske.sk", workers>0 ? Form("workers=%d",workers) : "");
+      p->Exec("TObject *o = gEnv->GetTable()->FindObject(\"Proof.UseMergers\"); gEnv->GetTable()->Remove(o);", kTRUE);
+
+      gProof->EnablePackage("VO_ALICE@AliRoot::v4-21-29-AN");
+      gSystem->Load("libCore.so");
+      gSystem->Load("libTree.so");
+      gSystem->Load("libGeom.so");
+      gSystem->Load("libVMC.so");
+      gSystem->Load("libPhysics.so");
+      gSystem->Load("libSTEERBase");
+      gSystem->Load("libESD");
+      gSystem->Load("libAOD");
+      gSystem->Load("libANALYSIS");
+      gSystem->Load("libOADB");
+      gSystem->Load("libANALYSISalice");
+
+      // Enable the needed package
+      // gProof->UploadPackage("$ALICE_ROOT/obj/STEERBase");
+      // gProof->EnablePackage("$ALICE_ROOT/obj/STEERBase");
+      // gProof->UploadPackage("$ALICE_ROOT/obj/ESD");
+      // gProof->EnablePackage("$ALICE_ROOT/obj/ESD");
+      // gProof->UploadPackage("$ALICE_ROOT/obj/AOD");
+      // gProof->EnablePackage("$ALICE_ROOT/obj/AOD");
+      // gProof->UploadPackage("$ALICE_ROOT/obj/ANALYSIS");
+      // gProof->EnablePackage("$ALICE_ROOT/obj/ANALYSIS");
+      // gProof->UploadPackage("$ALICE_ROOT/obj/OADB");
+      // gProof->EnablePackage("$ALICE_ROOT/obj/OADB");
+      // gProof->UploadPackage("$ALICE_ROOT/obj/ANALYSISalice");
+      // gProof->EnablePackage("$ALICE_ROOT/obj/ANALYSISalice");
+      // gProof->UploadPackage("$ALICE_ROOT/obj/PWG0base");
+      // gProof->EnablePackage("$ALICE_ROOT/obj/PWG0base");
+      // gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG0/multPb"));
+      // gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background"));
+   }
+   else
+   {
+      cout << "Init in Local or Grid mode" << endl;
+      gSystem->Load("libCore.so");
+      gSystem->Load("libTree.so");
+      gSystem->Load("libGeom.so");
+      gSystem->Load("libVMC.so");
+      gSystem->Load("libPhysics.so");
+      gSystem->Load("libSTEERBase");
+      gSystem->Load("libESD");
+      gSystem->Load("libAOD");
+      gSystem->Load("libANALYSIS");
+      gSystem->Load("libOADB");
+      gSystem->Load("libANALYSISalice");
+      // Use AliRoot includes to compile our task
+      gROOT->ProcessLine(".include $ALICE_ROOT/include");
+
+      // gSystem->Load("libVMC");
+      // gSystem->Load("libTree");
+      // gSystem->Load("libSTEERBase");
+      // gSystem->Load("libESD");
+      // gSystem->Load("libAOD");
+      // gSystem->Load("libANALYSIS");
+      // gSystem->Load("libANALYSISalice");
+      // gSystem->Load("libPWG0base");
+
+      // gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG0/multPb"));
+      // gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background"));
+      //    gROOT->ProcessLine(gSystem->ExpandPathName(".include $ALICE_ROOT/PWG1/background/"));
+   }
+   // Load helper classes
+   TIterator * iter = listToLoad->MakeIterator();
+   TObjString * name = 0;
+   while (name = (TObjString *)iter->Next())
+   {
+      gSystem->ExpandPathName(name->String());
+      cout << name->String().Data();
+      if (runMode == kMyRunModeCAF)
+      {
+         gProof->Load(name->String() + (debug ? "+g" : ""));
+      }
+      else
+      {
+         gROOT->LoadMacro(name->String() + (debug ? "+g" : ""));
+      }
+   }
+
+}