i) Adding the centrality train for the MC (taking the impact parameter from the heade...
authorpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 23 Oct 2011 22:27:10 +0000 (22:27 +0000)
committerpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 23 Oct 2011 22:27:10 +0000 (22:27 +0000)
PWG2/EBYE/AliAnalysisTaskBF.cxx
PWG2/EBYE/AliAnalysisTaskBF.h
PWG2/EBYE/macros/AddTaskBalanceMCCentralityTrain.C [new file with mode: 0644]
PWG2/EBYE/macros/runBalanceFunctionMC.C [new file with mode: 0755]

index 9bce442..717aa30 100755 (executable)
@@ -5,6 +5,7 @@
 #include "TGraphErrors.h"\r
 #include "TH1F.h"\r
 #include "TH2F.h"\r
+#include "TArrayF.h"\r
 \r
 #include "AliAnalysisTaskSE.h"\r
 #include "AliAnalysisManager.h"\r
@@ -15,6 +16,8 @@
 #include "AliAODEvent.h"\r
 #include "AliAODTrack.h"\r
 #include "AliAODInputHandler.h"\r
+#include "AliGenEventHeader.h"\r
+#include "AliGenHijingEventHeader.h"\r
 #include "AliMCEventHandler.h"\r
 #include "AliMCEvent.h"\r
 #include "AliStack.h"\r
@@ -25,7 +28,7 @@
 \r
 \r
 // Analysis task for the BF code\r
-// Authors: Panos Cristakoglou@cern.ch\r
+// Authors: Panos.Christakoglou@nikhef.nl\r
 \r
 ClassImp(AliAnalysisTaskBF)\r
 \r
@@ -58,6 +61,11 @@ AliAnalysisTaskBF::AliAnalysisTaskBF(const char *name)
   fUseCentrality(kFALSE),\r
   fCentralityPercentileMin(0.), \r
   fCentralityPercentileMax(5.),\r
+  fImpactParameterMin(0.),\r
+  fImpactParameterMax(20.),\r
+  fUseMultiplicity(kFALSE),\r
+  fNumberOfAcceptedTracksMin(0),\r
+  fNumberOfAcceptedTracksMax(10000),\r
   fUseOfflineTrigger(kFALSE),\r
   fVxMax(0.3),\r
   fVyMax(0.3),\r
@@ -247,6 +255,7 @@ void AliAnalysisTaskBF::UserExec(Option_t *) {
   TObjArray *array         = new TObjArray();\r
   AliESDtrack *track_TPC   = NULL;\r
 \r
+  Int_t gNumberOfAcceptedTracks = 0;\r
 \r
   // vector holding the charges of all tracks\r
   vector<Int_t> chargeVectorShuffle;   // this will be shuffled\r
@@ -375,7 +384,7 @@ void AliAnalysisTaskBF::UserExec(Option_t *) {
 \r
     // store offline trigger bits\r
     fHistTriggerStats->Fill(aodHeader->GetOfflineTrigger());\r
-\r
+    \r
     // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
     fHistEventStats->Fill(1); //all events\r
     Bool_t isSelected = kTRUE;\r
@@ -496,6 +505,7 @@ void AliAnalysisTaskBF::UserExec(Option_t *) {
                    fHistEta->Fill(eta);\r
                    fHistPhi->Fill(aodTrack->Phi()*TMath::RadToDeg());\r
                    \r
+                   gNumberOfAcceptedTracks += 1;\r
                    // fill BF array\r
                    array->Add(aodTrack);\r
                    \r
@@ -514,37 +524,203 @@ void AliAnalysisTaskBF::UserExec(Option_t *) {
     }//triggered event \r
   }//AOD analysis\r
 \r
-  //MC analysis\r
-  else if(gAnalysisLevel == "MC") {\r
-    \r
+  //MC-ESD analysis\r
+  if(gAnalysisLevel == "MCESD") {\r
     AliMCEvent*  mcEvent = MCEvent(); \r
     if (!mcEvent) {\r
       Printf("ERROR: mcEvent not available");\r
       return;\r
     }\r
-    \r
-    Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());\r
-    for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {\r
-      AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));\r
-      if (!track) {\r
-       Printf("ERROR: Could not receive particle %d", iTracks);\r
-       continue;\r
-      }\r
 \r
-      if( track->Pt() < fPtMin || track->Pt() > fPtMax)      continue;\r
-      if( track->Eta() < fEtaMin || track->Eta() > fEtaMax)  continue;\r
+    AliESDEvent* gESD = dynamic_cast<AliESDEvent*>(InputEvent()); // from TaskSE\r
+    if (!gESD) {\r
+      Printf("ERROR: gESD not available");\r
+      return;\r
+    }\r
+\r
+    // store offline trigger bits\r
+    fHistTriggerStats->Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());\r
 \r
-      array->Add(track);\r
+    // event selection done in AliAnalysisTaskSE::Exec() --> this is not used\r
+    fHistEventStats->Fill(1); //all events\r
+    Bool_t isSelected = kTRUE;\r
+    if(fUseOfflineTrigger)\r
+      isSelected = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();\r
+    if(isSelected) {\r
+      fHistEventStats->Fill(2); //triggered events\r
 \r
-      // fill charge vector\r
-      chargeVector.push_back(track->Charge());\r
-      if(fRunShuffling) {\r
-       chargeVectorShuffle.push_back(track->Charge());\r
+      if(fUseCentrality) {\r
+       //Centrality stuff\r
+       AliCentrality *centrality = gESD->GetCentrality();\r
+       //Int_t nCentrality = 0;\r
+       //nCentrality = centrality->GetCentralityClass5(fCentralityEstimator.Data());\r
+       //cout<<nCentrality<<" "<<centrality->IsEventInCentralityClass(fCentralityPercentileMin,fCentralityPercentileMax,fCentralityEstimator.Data())<<endl;\r
+       \r
+       // take only events inside centrality class\r
+       if(!centrality->IsEventInCentralityClass(fCentralityPercentileMin,\r
+                                                fCentralityPercentileMax,\r
+                                                fCentralityEstimator.Data()))\r
+         return;\r
+       \r
+       // centrality QA (V0M)\r
+       fHistV0M->Fill(gESD->GetVZEROData()->GetMTotV0A(), gESD->GetVZEROData()->GetMTotV0C());\r
       }\r
+       \r
+      const AliESDVertex *vertex = gESD->GetPrimaryVertex();\r
+      if(vertex) {\r
+       if(vertex->GetNContributors() > 0) {\r
+         if(vertex->GetZRes() != 0) {\r
+           fHistEventStats->Fill(3); //events with a proper vertex\r
+           if(TMath::Abs(vertex->GetXv()) < fVxMax) {\r
+             if(TMath::Abs(vertex->GetYv()) < fVyMax) {\r
+               if(TMath::Abs(vertex->GetZv()) < fVzMax) {\r
+                 fHistEventStats->Fill(4); //analayzed events\r
+                 fHistVx->Fill(vertex->GetXv());\r
+                 fHistVy->Fill(vertex->GetYv());\r
+                 fHistVz->Fill(vertex->GetZv());\r
+                 \r
+                 //Printf("There are %d tracks in this event", gESD->GetNumberOfTracks());\r
+                 for (Int_t iTracks = 0; iTracks < gESD->GetNumberOfTracks(); iTracks++) {\r
+                   AliESDtrack* track = dynamic_cast<AliESDtrack *>(gESD->GetTrack(iTracks));\r
+                   if (!track) {\r
+                     Printf("ERROR: Could not receive track %d", iTracks);\r
+                     continue;\r
+                   }   \r
+                   \r
+                   Int_t label = TMath::Abs(track->GetLabel());\r
+                   if(label > mcEvent->GetNumberOfTracks()) continue;\r
+                   if(label > mcEvent->GetNumberOfPrimaries()) continue;\r
+                   \r
+                   AliMCParticle* mcTrack = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(label));\r
+                   if(!mcTrack) continue;\r
+\r
+                   // take only TPC only tracks\r
+                   track_TPC   = new AliESDtrack();\r
+                   if(!track->FillTPCOnlyTrack(*track_TPC)) continue;\r
+                   \r
+                   //ESD track cuts\r
+                   if(fESDtrackCuts) \r
+                     if(!fESDtrackCuts->AcceptTrack(track_TPC)) continue;\r
+                   \r
+                   // fill QA histograms\r
+                   Float_t b[2];\r
+                   Float_t bCov[3];\r
+                   track_TPC->GetImpactParameters(b,bCov);\r
+                   if (bCov[0]<=0 || bCov[2]<=0) {\r
+                     AliDebug(1, "Estimated b resolution lower or equal zero!");\r
+                     bCov[0]=0; bCov[2]=0;\r
+                   }\r
+                   \r
+                   Int_t nClustersTPC = -1;\r
+                   nClustersTPC = track_TPC->GetTPCNclsIter1();   // TPC standalone\r
+                   //nClustersTPC = track->GetTPCclusters(0);   // global track\r
+                   Float_t chi2PerClusterTPC = -1;\r
+                   if (nClustersTPC!=0) {\r
+                     chi2PerClusterTPC = track_TPC->GetTPCchi2Iter1()/Float_t(nClustersTPC);      // TPC standalone\r
+                     //chi2PerClusterTPC = track->GetTPCchi2()/Float_t(nClustersTPC);     // global track\r
+                   }\r
+                   \r
+                   fHistClus->Fill(track_TPC->GetITSclusters(0),nClustersTPC);\r
+                   fHistDCA->Fill(b[1],b[0]);\r
+                   fHistChi2->Fill(chi2PerClusterTPC);\r
+                   fHistPt->Fill(mcTrack->Pt());\r
+                   fHistEta->Fill(mcTrack->Eta());\r
+                   fHistPhi->Fill(mcTrack->Phi()*TMath::RadToDeg());\r
+                   \r
+                   // fill BF array\r
+                   array->Add(track_TPC);\r
+                   \r
+                   // fill charge vector\r
+                   chargeVector.push_back(track_TPC->Charge());\r
+                   if(fRunShuffling){\r
+                     chargeVectorShuffle.push_back(track_TPC->Charge());\r
+                   }\r
+                   \r
+                   delete track_TPC;\r
+                   \r
+                 } //track loop\r
+               }//Vz cut\r
+             }//Vy cut\r
+           }//Vx cut\r
+         }//proper vertex resolution\r
+       }//proper number of contributors\r
+      }//vertex object valid\r
+    }//triggered event \r
+  }//MC-ESD analysis\r
 \r
-    } //track loop\r
+  //MC analysis\r
+  else if(gAnalysisLevel == "MC") {\r
+    AliMCEvent*  mcEvent = MCEvent(); \r
+    if (!mcEvent) {\r
+      Printf("ERROR: mcEvent not available");\r
+      return;\r
+    }\r
+\r
+    Double_t gReactionPlane = 0., gImpactParameter = 0.;\r
+    if(fUseCentrality) {\r
+      //Get the MC header\r
+      AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());\r
+      if (headerH) {\r
+       //Printf("=====================================================");\r
+       //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());\r
+       //Printf("Impact parameter: %lf",headerH->ImpactParameter());\r
+       //Printf("=====================================================");\r
+       gReactionPlane = headerH->ReactionPlaneAngle();\r
+       gImpactParameter = headerH->ImpactParameter();\r
+      }\r
+      // take only events inside centrality class\r
+      if((fImpactParameterMin > gImpactParameter) || (fImpactParameterMax < gImpactParameter))\r
+       return;\r
+    }\r
+    \r
+    AliGenEventHeader *header = mcEvent->GenEventHeader();\r
+    if(!header) return;\r
+    \r
+    TArrayF gVertexArray;\r
+    header->PrimaryVertex(gVertexArray);\r
+    //Printf("Vertex: %lf (x) - %lf (y) - %lf (z)",\r
+    //gVertexArray.At(0),\r
+    //gVertexArray.At(1),\r
+    //gVertexArray.At(2));\r
+    fHistEventStats->Fill(3); //events with a proper vertex\r
+    if(TMath::Abs(gVertexArray.At(0)) < fVxMax) {\r
+      if(TMath::Abs(gVertexArray.At(1)) < fVyMax) {\r
+       if(TMath::Abs(gVertexArray.At(2)) < fVzMax) {\r
+         fHistEventStats->Fill(4); //analayzed events\r
+         fHistVx->Fill(gVertexArray.At(0));\r
+         fHistVy->Fill(gVertexArray.At(1));\r
+         fHistVz->Fill(gVertexArray.At(2));\r
+         \r
+         Printf("There are %d tracks in this event", mcEvent->GetNumberOfPrimaries());\r
+         for (Int_t iTracks = 0; iTracks < mcEvent->GetNumberOfPrimaries(); iTracks++) {\r
+           AliMCParticle* track = dynamic_cast<AliMCParticle *>(mcEvent->GetTrack(iTracks));\r
+           if (!track) {\r
+             Printf("ERROR: Could not receive particle %d", iTracks);\r
+             continue;\r
+           }\r
+           \r
+           if( track->Pt() < fPtMin || track->Pt() > fPtMax)      \r
+             continue;\r
+           if( track->Eta() < fEtaMin || track->Eta() > fEtaMax)  \r
+             continue;\r
+           \r
+           array->Add(track);\r
+           \r
+           // fill charge vector\r
+           chargeVector.push_back(track->Charge());\r
+           if(fRunShuffling) {\r
+             chargeVectorShuffle.push_back(track->Charge());\r
+           }\r
+         } //track loop\r
+       }//Vz cut\r
+      }//Vy cut\r
+    }//Vx cut\r
   }//MC analysis\r
   \r
+  if(fUseMultiplicity) \r
+    if((gNumberOfAcceptedTracks < fNumberOfAcceptedTracksMin)||(gNumberOfAcceptedTracks > fNumberOfAcceptedTracksMax))\r
+      return;\r
+  \r
   // calculate balance function\r
   fBalance->CalculateBalance(array,chargeVector);\r
 \r
@@ -555,8 +731,6 @@ void AliAnalysisTaskBF::UserExec(Option_t *) {
   }\r
   \r
   delete array;\r
-  \r
-\r
 }      \r
 \r
 //________________________________________________________________________\r
index df68359..171bbc4 100755 (executable)
@@ -72,8 +72,19 @@ class AliAnalysisTaskBF : public AliAnalysisTaskSE {
     fCentralityPercentileMin=min;\r
     fCentralityPercentileMax=max;\r
   }\r
+  void SetImpactParameterRange(Double_t min, Double_t max) { \r
+    fUseCentrality = kTRUE;\r
+    fImpactParameterMin=min;\r
+    fImpactParameterMax=max;\r
+  }\r
 \r
-  void UseOfflineTrigger() {fUseOfflineTrigger = kTRUE;}\r
+  //multiplicity\r
+  void SetCentralityPercentileRange(Int_t min, Int_t max) {\r
+    fUseMultiplicity = kTRUE;\r
+    fNumberOfAcceptedTracksMin = min;\r
+    fNumberOfAcceptedTracksMax = max;}\r
+  \r
+    void UseOfflineTrigger() {fUseOfflineTrigger = kTRUE;}\r
 \r
  private:\r
   AliBalance *fBalance; //BF object\r
@@ -104,8 +115,14 @@ class AliAnalysisTaskBF : public AliAnalysisTaskSE {
 \r
   TString fCentralityEstimator;      //"V0M","TRK","TKL","ZDC","FMD"\r
   Bool_t fUseCentrality;//use the centrality (PbPb) or not (pp)\r
-  Double_t fCentralityPercentileMin;\r
-  Double_t fCentralityPercentileMax;\r
+  Double_t fCentralityPercentileMin;//centrality percentile min\r
+  Double_t fCentralityPercentileMax;//centrality percentile max\r
+  Double_t fImpactParameterMin;//impact parameter min (used for MC)\r
+  Double_t fImpactParameterMax;//impact parameter max (used for MC)\r
+\r
+  Bool_t fUseMultiplicity;//use the multiplicity cuts\r
+  Int_t fNumberOfAcceptedTracksMin;//min. number of number of accepted tracks (used for the multiplicity dependence study - pp)\r
+  Int_t fNumberOfAcceptedTracksMax;//max. number of number of accepted tracks (used for the multiplicity dependence study - pp)\r
 \r
   Bool_t fUseOfflineTrigger;//Usage of the offline trigger selection\r
 \r
diff --git a/PWG2/EBYE/macros/AddTaskBalanceMCCentralityTrain.C b/PWG2/EBYE/macros/AddTaskBalanceMCCentralityTrain.C
new file mode 100644 (file)
index 0000000..efc5b6b
--- /dev/null
@@ -0,0 +1,134 @@
+// now in options\r
+//=============================================//\r
+const char* centralityEstimator = "V0M";\r
+//const char* centralityEstimator = "CL1";\r
+//const char* centralityEstimator = "TRK";\r
+//=============================================//\r
+//Bool_t gRunShuffling = kFALSE;\r
+//Bool_t gRunShuffling = kTRUE;\r
+//=============================================//\r
+//_________________________________________________________//\r
+AliAnalysisTaskBF *AddTaskBalanceMCCentralityTrain(Double_t centrMin=0.,\r
+                                                  Double_t centrMax=100.,\r
+                                                  Double_t impactParameterMin=0.,\r
+                                                  Double_t impactParameterMax=20.,\r
+                                                  Bool_t gRunShuffling=kFALSE,\r
+                                                  Double_t vertexZ=10.,\r
+                                                  Double_t ptMin=0.3,\r
+                                                  Double_t ptMax=1.5,\r
+                                                  Double_t etaMin=-0.8,\r
+                                                  Double_t etaMax=0.8,\r
+                                                  TString fileNameBase="AnalysisResults") {\r
+\r
+  // Creates a balance function analysis task and adds it to the analysis manager.\r
+  // Get the pointer to the existing analysis manager via the static access method.\r
+  TString centralityName("");\r
+  centralityName+=Form("%.0f",centrMin);\r
+  centralityName+="-";\r
+  centralityName+=Form("%.0f",centrMax);\r
+\r
+  TString outputFileName(fileNameBase);\r
+  outputFileName.Append(".root");\r
+\r
+  //===========================================================================\r
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();\r
+  if (!mgr) {\r
+    ::Error("AddTaskBF", "No analysis manager to connect to.");\r
+    return NULL;\r
+  }\r
+\r
+  // Check the analysis type using the event handlers connected to the analysis manager.\r
+  //===========================================================================\r
+  if (!mgr->GetInputEventHandler()) {\r
+    ::Error("AddTaskBF", "This task requires an input event handler");\r
+    return NULL;\r
+  }\r
+  TString analysisType = "MC";\r
+\r
+  // for local changed BF configuration\r
+  //gROOT->LoadMacro("./configBalanceFunctionAnalysis.C");\r
+  gROOT->LoadMacro("$ALICE_ROOT/PWG2/EBYE/macros/configBalanceFunctionAnalysis.C");\r
+  AliBalance *bf  = 0;  // Balance Function object\r
+  AliBalance *bfs = 0;  // shuffled Balance function object\r
+\r
+  if (analysisType=="ESD"){\r
+    bf  = GetBalanceFunctionObject("ESD");\r
+    if(gRunShuffling) bfs = GetBalanceFunctionObject("ESD",kTRUE);\r
+  }\r
+  else if (analysisType=="AOD"){\r
+    bf  = GetBalanceFunctionObject("AOD");\r
+    if(gRunShuffling) bfs = GetBalanceFunctionObject("AOD",kTRUE);\r
+  }\r
+  else if (analysisType=="MC"){\r
+    bf  = GetBalanceFunctionObject("MC");\r
+    if(gRunShuffling) bfs = GetBalanceFunctionObject("MC",kTRUE);\r
+  }\r
+  else{\r
+    ::Error("AddTaskBF", "analysis type NOT known.");\r
+    return NULL;\r
+  }\r
+\r
+  // Create the task, add it to manager and configure it.\r
+  //===========================================================================\r
+  AliAnalysisTaskBF *taskBF = new AliAnalysisTaskBF("TaskBF");\r
+  taskBF->SetAnalysisObject(bf);\r
+  if(gRunShuffling) taskBF->SetShufflingObject(bfs);\r
+\r
+  if(analysisType == "ESD") {\r
+    AliESDtrackCuts *trackCuts = GetTrackCutsObject(ptMin,ptMax,etaMin,etaMax,maxTPCchi2,DCAxy,DCAz,minNClustersTPC);\r
+    taskBF->SetAnalysisCutObject(trackCuts);\r
+    // centrality estimator (default = V0M)\r
+    taskBF->SetCentralityEstimator(centralityEstimator);\r
+    taskBF->SetCentralityPercentileRange(impactParameterMin,\r
+                                        impactParameterMax);\r
+  }\r
+  else if(analysisType == "AOD") {\r
+    // pt and eta cut (pt_min, pt_max, eta_min, eta_max)\r
+    taskBF->SetAODtrackCutBit(128);\r
+    taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax);\r
+    taskBF->SetCentralityEstimator(centralityEstimator);\r
+    taskBF->SetCentralityPercentileRange(impactParameterMin,\r
+                                        impactParameterMax);\r
+\r
+    // set extra DCA cuts (-1 no extra cut)\r
+    taskBF->SetExtraDCACutsAOD(DCAxy,DCAz);\r
+\r
+    // set extra TPC chi2 / nr of clusters cut\r
+    taskBF->SetExtraTPCCutsAOD(maxTPCchi2, minNClustersTPC);\r
+    taskBF->SetCentralityEstimator(centralityEstimator);    \r
+  }\r
+  else if(analysisType == "MC") {\r
+    taskBF->SetKinematicsCutsAOD(ptMin,ptMax,etaMin,etaMax); \r
+    taskBF->SetImpactParameterRange(impactParameterMin,\r
+                                   impactParameterMax);\r
+  }\r
+\r
+  // offline trigger selection (AliVEvent.h)\r
+  // taskBF->UseOfflineTrigger(); // NOT used (selection is done with the AliAnalysisTaskSE::SelectCollisionCandidates()) \r
+  // with this only selected events are analyzed (first 2 bins in event QA histogram are the same))\r
+  // documentation in https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PWG1EvSelDocumentation\r
+  //taskBF->SelectCollisionCandidates(AliVEvent::kMB);\r
+    \r
+  // vertex cut (x,y,z)\r
+  taskBF->SetVertexDiamond(.3,.3,vertexZ);\r
+  \r
+\r
+\r
+  //bf->PrintAnalysisSettings();\r
+  mgr->AddTask(taskBF);\r
+  \r
+  // Create ONLY the output containers for the data produced by the task.\r
+  // Get and connect other common input/output containers via the manager as below\r
+  //==============================================================================\r
+  TString outputFileName = AliAnalysisManager::GetCommonFileName();\r
+  outputFileName += ":PWG2EbyE.outputBalanceFunctionAnalysis";\r
+  AliAnalysisDataContainer *coutQA = mgr->CreateContainer(Form("listQA_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  AliAnalysisDataContainer *coutBF = mgr->CreateContainer(Form("listBF_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  if(gRunShuffling) AliAnalysisDataContainer *coutBFS= mgr->CreateContainer(Form("listBFShuffled_%s",centralityName.Data()), TList::Class(),AliAnalysisManager::kOutputContainer,outputFileName.Data());\r
+  mgr->ConnectInput(taskBF, 0, mgr->GetCommonInputContainer());\r
+  mgr->ConnectOutput(taskBF, 1, coutQA);\r
+  mgr->ConnectOutput(taskBF, 2, coutBF);\r
+  if(gRunShuffling) mgr->ConnectOutput(taskBF, 3, coutBFS);\r
+\r
+  return taskBF;\r
+}\r
diff --git a/PWG2/EBYE/macros/runBalanceFunctionMC.C b/PWG2/EBYE/macros/runBalanceFunctionMC.C
new file mode 100755 (executable)
index 0000000..f3e5e50
--- /dev/null
@@ -0,0 +1,246 @@
+enum analysisModes {mLocal,mLocalPAR,mPROOF,mGrid,mGridPAR};\r
+enum analysisTypes {mESD,mAOD,mMC,mMCESD};\r
+\r
+//\r
+class AliAnalysisGrid;\r
+class AliAnalysisTaskBF;\r
+class AliBalance;\r
+\r
+//Centrality stuff\r
+Int_t binfirst = 0;  //where do we start numbering bins\r
+Int_t binlast = 8;  //where do we stop numbering bins\r
+const Int_t numberOfCentralityBins = 9;\r
+Float_t centralityArray[numberOfCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.}; // in centrality percentile\r
+Float_t impactParameterArray[numberOfCentralityBins+1] = {0.0,3.79,5.30,7.41,9.04,10.40,11.61,12.68,13.67,14.63}; // in fm (impact parametger taken from the MC header)\r
+\r
+//________________________________________________________________________//\r
+void runBalanceFunctionMC(Int_t mode = mLocal, \r
+                         Int_t type = mMC,\r
+                         Bool_t DATA = kFALSE) {\r
+  // Time:\r
+  TStopwatch timer;\r
+  timer.Start();\r
+  \r
+  //Check analysis mode\r
+  if((mode < 0) || (mode > 4)) {\r
+    Printf("Analysis mode not recognized!");\r
+    Printf("You can select out of 0: local, 1: local with par files, 2: proof, 3: grid, 4: grid with par files");\r
+    return;\r
+  }\r
+  \r
+  //Check analysis type\r
+  if((type < 0) || (type > 3)) {\r
+    Printf("Analysis type not recognized!");\r
+    Printf("You can select out of 0: ESD, 1: AOD, 2: MC (stack), 3: MC (from the ESD)");\r
+    return;\r
+  }\r
+  \r
+  // Load needed libraries:\r
+  LoadLibraries(mode);\r
+  \r
+\r
+  // Create and configure the AliEn plug-in:\r
+  if(mode == mGrid || mode == mGridPAR) {\r
+    gROOT->LoadMacro("CreateAlienHandler.C");\r
+    AliAnalysisGrid *alienHandler = CreateAlienHandler(runListFileName);  \r
+    if (!alienHandler) return;\r
+  }\r
+  // Chains:   \r
+  if(mode == mLocal || mode == mLocalPAR) {\r
+    TChain* chain = 0x0;\r
+    if((type == mESD)||(type == mMCESD))  \r
+      chain = new TChain("esdTree");\r
+    else if(type == mAOD)\r
+      chain = new TChain("aodTree");\r
+    else if(type == mMC)\r
+      chain = new TChain("TE");\r
+\r
+    TString filename;\r
+    for(Int_t i = 1; i < 10; i++) {\r
+      filename = "/data/alice2/pchrist/HIJING/Full/2.76TeV/";\r
+      filename += "Set"; filename += i; \r
+      if((type == mESD)||(type == mMCESD))  \r
+       filename += "/AliESDs.root";\r
+      else if(type == mAOD)\r
+       filename += "/AliAOD.root";\r
+      else if(type == mMC)\r
+       filename += "/galice.root";\r
+\r
+      chain->Add(filename.Data());\r
+    }\r
+  }\r
+  //Proof\r
+  if(mode == mPROOF) {\r
+    gROOT->ProcessLine(Form(".include %s/include", gSystem->ExpandPathName("$ALICE_ROOT")));\r
+  }\r
+  \r
+  // analysis manager\r
+  AliAnalysisManager* mgr = new AliAnalysisManager("balanceFunctionManager");\r
+  if(mode == mGrid || mode == mGridPAR)\r
+    mgr->SetGridHandler(alienHandler);\r
+    \r
+  // input handler (ESD or AOD)\r
+  AliVEventHandler* inputH = NULL;\r
+  if((type == mESD)||(type == mMCESD)||(type == mMC))  \r
+    inputH = new AliESDInputHandler();\r
+  else if(type == mAOD)\r
+    inputH = new AliAODInputHandler();\r
+  mgr->SetInputEventHandler(inputH);\r
+    \r
+  // mc event handler\r
+  if((type == mMC) || (type == mMCESD)) {\r
+    AliMCEventHandler* mchandler = new AliMCEventHandler();\r
+    // Not reading track references\r
+    mchandler->SetReadTR(kFALSE);\r
+    mgr->SetMCtruthEventHandler(mchandler);\r
+  }   \r
+\r
+  if(mAOD){\r
+    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C");\r
+    AliCentralitySelectionTask *taskCentrality = AddTaskCentrality();\r
+\r
+    // Add physics selection task (NOT needed for AODs)\r
+    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");\r
+    AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(DATA);\r
+  }\r
+\r
+  //Add the BF task (all centralities)\r
+  gROOT->LoadMacro("AddTaskBalanceMCCentralityTrain.C"); \r
+  for (Int_t i=binfirst; i<binlast+1; i++) {\r
+    Float_t lowCentralityBinEdge = centralityArray[i];\r
+    Float_t highCentralityBinEdge = centralityArray[i+1];\r
+    Printf("\nWagon for centrality bin %i: %.0f-%.0f",i,lowCentralityBinEdge,highCentralityBinEdge);\r
+    AddTaskBalanceMCCentralityTrain(lowCentralityBinEdge,\r
+                                   highCentralityBinEdge,\r
+                                   impactParameterArray[i],\r
+                                   impactParameterArray[i+1],\r
+                                   kTRUE,\r
+                                   10.,0.3,1.5,-0.8,0.8);\r
+  }\r
+\r
+  // enable debug printouts\r
+  mgr->SetDebugLevel(2);\r
+  mgr->SetUseProgressBar(1,100);\r
+  if (!mgr->InitAnalysis()) return;\r
+  mgr->PrintStatus();\r
+  \r
+  // start analysis\r
+  if(mode == mLocal || mode == mLocalPAR) \r
+    mgr->StartAnalysis("local",chain);\r
+  else if(mode == mPROOF) \r
+    mgr->StartAnalysis("proof",dataDir,nRuns,offset);\r
+  else if(mode == mGrid || mode == mGridPAR) \r
+    mgr->StartAnalysis("grid");\r
+  \r
+  // Print real and CPU time used for analysis:  \r
+  timer.Stop();\r
+  timer.Print();\r
+}\r
+\r
+//=============================================================//\r
+void LoadLibraries(const analysisModes mode) {  \r
+  //--------------------------------------\r
+  // Load the needed libraries most of them already loaded by aliroot\r
+  //--------------------------------------\r
+  gSystem->Load("libCore.so");        \r
+  gSystem->Load("libGeom.so");\r
+  gSystem->Load("libVMC.so");\r
+  gSystem->Load("libPhysics.so");\r
+  gSystem->Load("libTree.so");\r
+\r
+  //----------------------------------------------------------\r
+  // >>>>>>>>>>> Local mode <<<<<<<<<<<<<< \r
+  //----------------------------------------------------------\r
+  if (mode==mLocal || mode==mGrid || mode == mGridPAR) {\r
+    //--------------------------------------------------------\r
+    // If you want to use already compiled libraries \r
+    // in the aliroot distribution\r
+    //--------------------------------------------------------\r
+    gSystem->Load("libSTEERBase.so");\r
+    gSystem->Load("libESD.so");\r
+    gSystem->Load("libAOD.so");\r
+    gSystem->Load("libANALYSIS.so");\r
+    gSystem->Load("libANALYSISalice.so");\r
+    gSystem->Load("libPWG2ebye.so");\r
+    // Use AliRoot includes to compile our task\r
+    gROOT->ProcessLine(".include $ALICE_ROOT/include");\r
+  }\r
+  \r
+  else if (mode == mLocalPAR) {\r
+    //--------------------------------------------------------\r
+    //If you want to use root and par files from aliroot\r
+    //--------------------------------------------------------  \r
+    SetupPar("STEERBase");\r
+    SetupPar("ESD");\r
+    SetupPar("AOD");\r
+    SetupPar("ANALYSIS");\r
+    SetupPar("ANALYSISalice");\r
+    SetupPar("PWG2ebye");\r
+}\r
+  \r
+  //---------------------------------------------------------\r
+  // <<<<<<<<<< PROOF mode >>>>>>>>>>>>\r
+  //---------------------------------------------------------\r
+  else if (mode==mPROOF) {\r
+    // Connect to proof\r
+    printf("*** Connect to PROOF ***\n");\r
+    gEnv->SetValue("XSec.GSI.DelegProxy","2");\r
+    // Put appropriate username here\r
+    TProof::Open("alice-caf.cern.ch");\r
+    //TProof::Open("skaf.saske.sk");\r
+    //TProof::Open("prf000-iep-grid.saske.sk");\r
+\r
+    gProof->EnablePackage("VO_ALICE@AliRoot::v4-21-12-AN");\r
+  }  \r
+  \r
+} // end of void LoadLibraries(const anaModes mode)\r
+\r
+//======================================================================//\r
+void SetupPar(char* pararchivename) {\r
+  //Load par files, create analysis libraries\r
+  //For testing, if par file already decompressed and modified\r
+  //classes then do not decompress.\r
+  \r
+  TString cdir(Form("%s", gSystem->WorkingDirectory() )) ; \r
+  TString parpar(Form("%s.par", pararchivename)) ; \r
+  if ( gSystem->AccessPathName(parpar.Data()) ) {\r
+    gSystem->ChangeDirectory(gSystem->Getenv("ALICE_ROOT")) ;\r
+    TString processline(Form(".! make %s", parpar.Data())) ; \r
+    gROOT->ProcessLine(processline.Data()) ;\r
+    gSystem->ChangeDirectory(cdir) ; \r
+    processline = Form(".! mv /tmp/%s .", parpar.Data()) ;\r
+    gROOT->ProcessLine(processline.Data()) ;\r
+  } \r
+  if ( gSystem->AccessPathName(pararchivename) ) {  \r
+    TString processline = Form(".! tar xvzf %s",parpar.Data()) ;\r
+    gROOT->ProcessLine(processline.Data());\r
+  }\r
+  \r
+  TString ocwd = gSystem->WorkingDirectory();\r
+  gSystem->ChangeDirectory(pararchivename);\r
+  \r
+  // check for BUILD.sh and execute\r
+  if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {\r
+    printf("*******************************\n");\r
+    printf("*** Building PAR archive    ***\n");\r
+    cout<<pararchivename<<endl;\r
+    printf("*******************************\n");\r
+    if (gSystem->Exec("PROOF-INF/BUILD.sh")) {\r
+      Error("runProcess","Cannot Build the PAR Archive! - Abort!");\r
+      return -1;\r
+    }\r
+  }\r
+  // check for SETUP.C and execute\r
+  if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {\r
+    printf("*******************************\n");\r
+    printf("*** Setup PAR archive       ***\n");\r
+    cout<<pararchivename<<endl;\r
+    printf("*******************************\n");\r
+    gROOT->Macro("PROOF-INF/SETUP.C");\r
+  }\r
+  \r
+  gSystem->ChangeDirectory(ocwd.Data());\r
+  printf("Current dir: %s\n", ocwd.Data());\r
+\r
+} // end of void SetupPar(char* pararchivename) \r
+\r