]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Added analysis task for effiency studies (original by Veronica Canoa Roman). Now...
authoresicking <esicking@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 7 Oct 2011 13:00:15 +0000 (13:00 +0000)
committeresicking <esicking@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 7 Oct 2011 13:00:15 +0000 (13:00 +0000)
test/vmctest/scripts/efficiency/AddTaskEfficiency.C [new file with mode: 0644]
test/vmctest/scripts/efficiency/AliAnalysisTaskEfficiency.cxx [new file with mode: 0644]
test/vmctest/scripts/efficiency/AliAnalysisTaskEfficiency.h [new file with mode: 0644]
test/vmctest/scripts/efficiency/CreateAnalysisPlugin.C [new file with mode: 0644]
test/vmctest/scripts/efficiency/SetupAnalysis.C [new file with mode: 0644]
test/vmctest/scripts/efficiency/doubleRatio.C [new file with mode: 0644]
test/vmctest/scripts/efficiency/macrosalida.C [new file with mode: 0644]
test/vmctest/scripts/efficiency/runALICE.C [new file with mode: 0644]

diff --git a/test/vmctest/scripts/efficiency/AddTaskEfficiency.C b/test/vmctest/scripts/efficiency/AddTaskEfficiency.C
new file mode 100644 (file)
index 0000000..c20ade4
--- /dev/null
@@ -0,0 +1,163 @@
+AliAnalysisTaskEfficiency * AddTaskEfficiency(Int_t runNumber)
+
+{
+  // Creates a QA task exploiting simple symmetries phi, eta +/-, charge ...
+  
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskEfficiency", "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("AddTaskEfficiency", "This task requires an input event handler");
+    return NULL;
+  }
+   
+  TString inputDataType = mgr->GetInputEventHandler()->GetDataType(); 
+  // can be "ESD" or "AOD"
+  
+  // Configure analysis
+  //===========================================================================
+  
+  //TASKS
+  //Task for global tracks
+  AliAnalysisTaskEfficiency *task0 = 
+    new AliAnalysisTaskEfficiency("AliAnalysisTaskEfficiency_Global");
+  task0->SelectCollisionCandidates(); // default setting: kMB = min bias trigger
+  task0->SetTrackType(0);// global tracks
+
+  // Task for ITS tracks
+  AliAnalysisTaskEfficiency *task1 = 
+    new AliAnalysisTaskEfficiency("AliAnalysisTaskEfficiency_ITS");
+  task1->SelectCollisionCandidates();
+  task1->SetTrackType(1);// ITS tracks
+
+  //Task for ITS_SA tracks
+  AliAnalysisTaskEfficiency *task1sa = 
+    new AliAnalysisTaskEfficiency("AliAnalysisTaskEfficiency_ITS_SA");
+  task1sa->SelectCollisionCandidates();
+  task1sa->SetTrackType(1);// ITS tracks
+   
+  //Task for TPC tracks 
+  AliAnalysisTaskEfficiency *task2 = 
+    new AliAnalysisTaskEfficiency("AliAnalysisTaskEfficiency_TPC");
+  task2->SelectCollisionCandidates();
+  task2->SetTrackType(2);// TPC only tracks
+  //CUTS
+  //cuts for global tracks
+  AliESDtrackCuts* esdTrackCutsL0 = new AliESDtrackCuts("AliESDtrackCuts0","Global");
+  esdTrackCutsL0->SetMinNClustersTPC(70);
+  esdTrackCutsL0->SetRequireTPCRefit(kTRUE);
+  esdTrackCutsL0->SetRequireITSRefit(kTRUE);
+  esdTrackCutsL0->SetMaxDCAToVertexXY(3.);
+  esdTrackCutsL0->SetMaxDCAToVertexZ(3.);
+  esdTrackCutsL0->SetAcceptKinkDaughters(kFALSE);
+   
+  //cuts for ITS tracks
+  AliESDtrackCuts* esdTrackCutsL1 = new AliESDtrackCuts("AliESDtrackCuts1","ITS");
+  esdTrackCutsL1->SetMaxDCAToVertexXY(3.);
+  esdTrackCutsL1->SetMaxDCAToVertexZ(3.);
+  esdTrackCutsL1->SetAcceptKinkDaughters(kFALSE);
+  esdTrackCutsL1->SetRequireITSRefit(kTRUE);
+  esdTrackCutsL1->SetRequireITSStandAlone(kTRUE); 
+
+  //cuts for ITS tracks SA
+  AliESDtrackCuts* esdTrackCutsL1sa = new AliESDtrackCuts("AliESDtrackCuts1sa","ITS_SA");
+  esdTrackCutsL1sa->SetMaxDCAToVertexXY(3.);
+  esdTrackCutsL1sa->SetMaxDCAToVertexZ(3.);
+  esdTrackCutsL1sa->SetAcceptKinkDaughters(kFALSE);
+  esdTrackCutsL1sa->SetRequireITSRefit(kTRUE);
+  esdTrackCutsL1sa->SetRequireITSPureStandAlone(kTRUE);
+   
+  //cuts for TPC tracks
+  AliESDtrackCuts* esdTrackCutsL2 = new AliESDtrackCuts("AliESDtrackCuts2","TPC");
+  esdTrackCutsL2->SetRequireTPCRefit(kFALSE);
+  esdTrackCutsL2->SetAcceptKinkDaughters(kFALSE);
+  //jacek's cuts:
+  esdTrackCutsL2->SetMinNClustersTPC(70);
+  // cut on max ncl=160 in Task
+  esdTrackCutsL2->SetMaxDCAToVertexXY(3.);
+  esdTrackCutsL2->SetMaxDCAToVertexZ(3.);
+  esdTrackCutsL2->SetMaxChi2PerClusterTPC(3.999);
+  //cut minChi=0 in task
+  //esdTrackCutsL2->SetPRange(0.15,16); // not needed for QA
+  //esdTrackCutsL2->SetEtaRange(-0.8, 0.7999); // not needed for QA
+  //cuts for ITS tracks
+  
+  //add cuts to tasks
+  task0->SetCuts(esdTrackCutsL0);
+  task1->SetCuts(esdTrackCutsL1);
+  task1sa->SetCuts(esdTrackCutsL1sa);
+  task2->SetCuts(esdTrackCutsL2);
+
+  // add tasks to manager
+  mgr->AddTask(task0);
+  mgr->AddTask(task1);
+  mgr->AddTask(task1sa);
+  mgr->AddTask(task2);
+     
+  //output container for tasks
+  AliAnalysisDataContainer *cout0    = 0;
+  AliAnalysisDataContainer *cout1    = 0;
+  AliAnalysisDataContainer *cout1sa  = 0;
+  AliAnalysisDataContainer *cout2    = 0;
+
+   
+  if(runNumber>0){ 
+    cout0   =  mgr->CreateContainer("QAHists_Global",TList::Class(),
+                                   AliAnalysisManager::kOutputContainer, 
+                                   Form("run%d.root",runNumber));
+    cout1   =  mgr->CreateContainer("QAHists_ITS",TList::Class(),
+                                   AliAnalysisManager::kOutputContainer,
+                                   Form("run%d.root",runNumber));
+    cout1sa =  mgr->CreateContainer("QAHists_ITS_SA",TList::Class(),
+                                   AliAnalysisManager::kOutputContainer, 
+                                   Form("run%d.root",runNumber));
+    cout2   =  mgr->CreateContainer("QAHists_TPC",TList::Class(),
+                                   AliAnalysisManager::kOutputContainer, 
+                                   Form("run%d.root",runNumber));
+  }
+   
+  else{
+    cout0   = mgr->CreateContainer("QAHists_Global",TList::Class(),
+                                  AliAnalysisManager::kOutputContainer, 
+                                  Form("%s:QAHists",AliAnalysisManager::
+                                       GetCommonFileName()));
+    cout1   = mgr->CreateContainer("QAHists_ITS",TList::Class(),
+                                  AliAnalysisManager::kOutputContainer, 
+                                  Form("%s:QAHists",AliAnalysisManager::
+                                       GetCommonFileName()));
+    cout1sa = mgr->CreateContainer("QAHists_ITS_SA",TList::Class(),
+                                  AliAnalysisManager::kOutputContainer, 
+                                  Form("%s:QAHists",AliAnalysisManager::
+                                       GetCommonFileName()));
+    cout2   = mgr->CreateContainer("QAHists_TPC",TList::Class(),
+                                  AliAnalysisManager::kOutputContainer, 
+                                  Form("%s:QAHists",AliAnalysisManager::
+                                       GetCommonFileName()));
+  }
+  
+  //connect input to task
+  mgr->ConnectInput  (task0,   0, mgr->GetCommonInputContainer());
+  mgr->ConnectInput  (task1,   0, mgr->GetCommonInputContainer());
+  mgr->ConnectInput  (task1sa, 0, mgr->GetCommonInputContainer());
+  mgr->ConnectInput  (task2,   0, mgr->GetCommonInputContainer());
+
+  //connect output to task
+  mgr->ConnectOutput (task0,   1, cout0);
+  mgr->ConnectOutput (task1,   1, cout1);
+  mgr->ConnectOutput (task1sa, 1, cout1sa);
+  mgr->ConnectOutput (task2,   1, cout2);
+
+  
+  return task0;
+
+}
+
+
diff --git a/test/vmctest/scripts/efficiency/AliAnalysisTaskEfficiency.cxx b/test/vmctest/scripts/efficiency/AliAnalysisTaskEfficiency.cxx
new file mode 100644 (file)
index 0000000..037f023
--- /dev/null
@@ -0,0 +1,1085 @@
+#include "TChain.h"\r
+#include "TTree.h"\r
+#include "TH1F.h"\r
+#include "TH2F.h"\r
+#include "TProfile.h"\r
+\r
+#include "TCanvas.h"\r
+#include "TList.h"\r
+#include "TParticle.h"\r
+#include "TParticlePDG.h"\r
+#include "TProfile.h"\r
+\r
+#include "AliAnalysisTask.h"\r
+#include "AliAnalysisManager.h"\r
+\r
+#include "AliESDEvent.h"\r
+#include "AliStack.h"\r
+#include "AliESDVertex.h"\r
+#include "AliESDInputHandler.h"\r
+#include "AliESDtrackCuts.h"\r
+#include "AliMultiplicity.h"\r
+\r
+#include "AliMCParticle.h"\r
+#include "AliMCEvent.h"\r
+#include "AliAnalysisTaskEfficiency.h"\r
+#include "AliExternalTrackParam.h"\r
+#include "AliTrackReference.h"\r
+#include "AliStack.h"\r
+#include "AliHeader.h"\r
+#include "AliGenEventHeader.h"\r
+#include "AliGenCocktailEventHeader.h"\r
+#include "AliGenPythiaEventHeader.h"\r
+#include "AliGenDPMjetEventHeader.h"\r
+\r
+// Analysis Task for basic QA on the ESD\r
+// Authors: Jan Fiete Grosse-Oetringhaus, Christian Klein-Boesing, Veronica Canoa, AM, Eva Sicking\r
+\r
+ClassImp(AliAnalysisTaskEfficiency)\r
+\r
+//________________________________________________________________________\r
+  AliAnalysisTaskEfficiency::AliAnalysisTaskEfficiency(const char *name) \r
+    : AliAnalysisTaskSE(name) \r
+    ,fFieldOn(kTRUE)\r
+    ,fHists(0)\r
+    ,fHistRECpt(0)\r
+    ,fHistMCpt(0)\r
+    ,fHistMCNRpt(0)\r
+    ,fHistFAKEpt(0)\r
+    ,fHistMULTpt(0)\r
+    ,fHistRECeta(0)\r
+    ,fHistMCeta(0)\r
+    ,fHistMCNReta(0)\r
+    ,fHistFAKEeta(0)\r
+    ,fHistMULTeta(0)\r
+    ,fHistRECphi(0)\r
+    ,fHistMCphi(0)\r
+    ,fHistMCNRphi(0)\r
+    ,fHistFAKEphi(0)\r
+    ,fHistMULTphi(0)\r
+    ,fHistRECHPTeta(0)\r
+    ,fHistMCHPTeta(0)\r
+    ,fHistMCNRHPTeta(0)\r
+    ,fHistFAKEHPTeta(0)\r
+    ,fHistRECHPTphi(0)\r
+    ,fHistMCHPTphi(0)\r
+    ,fHistMCNRHPTphi(0)\r
+    ,fHistFAKEHPTphi(0)\r
+    ,fHistRecMult(0)\r
+    ,fHistNCluster(0)  \r
+    ,fh1VtxEff(0) \r
+    ,fh2VtxTpcSpd(0)\r
+    ,fh2EtaCorrelation(0)      \r
+    ,fh2EtaCorrelationShift(0)  \r
+    ,fh2PhiCorrelation(0)      \r
+    ,fh2PhiCorrelationShift(0)  \r
+    ,fh2PtCorrelation(0)      \r
+    ,fh2PtCorrelationShift(0)\r
+    ,fHistDeltaphiprimaries(0)\r
+    ,fHistDeltaphisecondaries(0)           \r
+    ,fHistDeltaphireject(0) \r
+    ,fHistTPCITSdeltax(0)\r
+    ,fHistTPCITSdeltay(0)\r
+    ,fHistTPCITSdeltaz(0)\r
+    ,fHistTPCITSdeltar(0)\r
+    ,fHistTPCTRDdeltax(0)\r
+    ,fHistTPCTRDdeltay(0)\r
+    ,fHistTPCTRDdeltaz(0)\r
+    ,fHistTPCTRDdeltar(0)\r
+    ,fHistTPCTRDdeltaphi(0)\r
+    ,fPtPullsPos(0)\r
+    ,fPtPullsNeg(0)\r
+    ,fPtShiftPos(0)\r
+    ,fPtShiftNeg(0)\r
+    ,fTrackType(0)\r
+    ,fCuts(0)\r
+    ,fVertexRvsZ(0)\r
+    ,fVertexRvsZC(0)\r
+    ,fEtaMultiFluc(0)\r
+    ,fEtaMulti(0)\r
+    ,fEtaMultiH(0)\r
+    ,fPhiMultiFluc(0)\r
+    ,fPhiMulti(0)\r
+    ,fPhiMultiH(0)\r
+{\r
+  //\r
+  // Constructor\r
+  //\r
+\r
+  for(int i = 0;i< 2;++i){\r
+    fh2MultSpdChips[i] = 0;\r
+  }\r
+\r
+  for(int i = 0;i < kTotalAC;++i){\r
+    fh2PhiPadRow[i] = fh2PhiLayer[i] = 0;\r
+  }\r
+\r
+  for(int i = 0; i < kTotalVtx; ++i){ \r
+    fh2VertexCorrelation[i] = \r
+      fh2VertexCorrelationShift[i] = 0;\r
+    fh1VertexShift[i] =         \r
+      fh1VertexShiftNorm[i] = 0;      \r
+  }\r
+\r
+  for(int i = 0;i< 8;++i){\r
+    fHistRECptCharge[i]=0;\r
+    fHistMCptCharge[i]=0;\r
+    fHistMCNRptCharge[i]=0;\r
+    fHistFAKEptCharge[i]=0;\r
+\r
+    fHistRECetaCharge[i]=0;\r
+    fHistMCetaCharge[i]=0;\r
+    fHistMCNRetaCharge[i]=0;\r
+    fHistFAKEetaCharge[i]=0;\r
+\r
+    fHistRECphiCharge[i]=0;\r
+    fHistMCphiCharge[i]=0;\r
+    fHistMCNRphiCharge[i]=0;\r
+    fHistFAKEphiCharge[i]=0;\r
+\r
+  }\r
+\r
+  DefineOutput(1,  TList::Class());\r
+}\r
+\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskEfficiency::UserCreateOutputObjects()\r
+{\r
+  // Create histograms\r
+\r
+\r
+  fHists = new TList();\r
+  fHistRECpt   = new TH1F("fHistRECpt",  " p_{T} distribution: all reconstructed",        100, 0., 20.);\r
+  fHistMCpt    = new TH1F("fHistMCpt",   " p_{T} distribution: all MC",                   100, 0., 20.);\r
+  fHistMCNRpt  = new TH1F("fHistMCNRpt", " p_{T} distribution: all not-reconstructed MC", 100, 0., 20.);\r
+  fHistFAKEpt  = new TH1F("fHistFAKEpt", " p_{T} distribution: all Fake",                 100, 0., 20.);\r
+  fHistMULTpt  = new TH1F("fHistMULTpt", " p_{T} distribution: multiply rec.",            100, 0., 20.);\r
+    \r
+  fHistRECeta   = new TH1F("fHistRECeta",  " #eta distribution: all reconstructed",        100,-1.0,1.0);\r
+  fHistMCeta    = new TH1F("fHistMCeta",   " #eta distribution: all MC",                   100,-1.0,1.0);\r
+  fHistMCNReta  = new TH1F("fHistMCNReta", " #eta distribution: all not-reconstructed MC", 100,-1.0,1.0);\r
+  fHistFAKEeta  = new TH1F("fHistFAKEeta", " #eta distribution: all Fake",                 100,-1.0,1.0);\r
+  fHistMULTeta  = new TH1F("fHistMULTeta", " #eta distribution: multiply rec.",            100,-1.0,1.0);\r
+\r
+  fHistRECphi   = new TH1F("fHistRECphi",  " #phi distribution: all reconstructed",        314, 0., 6.28);\r
+  fHistMCphi    = new TH1F("fHistMCphi",   " #phi distribution: all MC",                   314, 0., 6.28);\r
+  fHistMCNRphi  = new TH1F("fHistMCNRphi", " #phi distribution: all not-reconstructed MC", 314, 0., 6.28);\r
+  fHistFAKEphi  = new TH1F("fHistFAKEphi", " #phi distribution: all Fake",                 314, 0., 6.28);\r
+  fHistMULTphi  = new TH1F("fHistMULTphi", " #phi distribution: multipli rec.",            314, 0., 6.28);\r
+\r
+  fHistRECHPTeta   = new TH1F("fHistRECHPTeta",  " #eta distribution: all reconstructed",        100,-1.0,1.0);\r
+  fHistMCHPTeta    = new TH1F("fHistMCHPTeta",   " #eta distribution: all MC",                   100,-1.0,1.0);\r
+  fHistMCNRHPTeta  = new TH1F("fHistMCNRHPTeta", " #eta distribution: all not-reconstructed MC", 100,-1.0,1.0);\r
+  fHistFAKEHPTeta  = new TH1F("fHistFAKEHPTeta", " #eta distribution: all Fake",                 100,-1.0,1.0);\r
+\r
+  fHistRECHPTphi   = new TH1F("fHistRECHPTphi",  " #phi distribution: all reconstructed",        314, 0., 6.28);\r
+  fHistMCHPTphi    = new TH1F("fHistMCHPTphi",   " #phi distribution: all MC",                   314, 0., 6.28);\r
+  fHistMCNRHPTphi  = new TH1F("fHistMCNRHPTphi", " #phi distribution: all not-reconstructed MC", 314, 0., 6.28);\r
+  fHistFAKEHPTphi  = new TH1F("fHistFAKEHPTphi", " #phi distribution: all Fake",                 314, 0., 6.28);\r
+\r
+  fHistRecMult     = new TH1F("fHistRecMult",  "Multiple reconstructed tracks", 50, 0., 50.);\r
+  fHistNCluster    = new TH1F("fHistNCluster", "Number of clusters for suspicious tracks", 300, 0., 300.);\r
+\r
+  // CKB\r
+\r
+  fh1VtxEff = new TH1F("fh1VtxEff","VtxEff TPC",4,-0.5,3.5);\r
+  fh1VtxEff->GetXaxis()->SetBinLabel(1,"NO TPC VTX");\r
+  fh1VtxEff->GetXaxis()->SetBinLabel(2,"TPC VTX");\r
+  fh1VtxEff->GetXaxis()->SetBinLabel(3,"NO SPD VTX");\r
+  fh1VtxEff->GetXaxis()->SetBinLabel(4,"SPD VTX");\r
+\r
+  TString labels[kTotalAC];\r
+  labels[kNegA]="NegA";\r
+  labels[kPosA]="PosA";\r
+  labels[kNegC]="NegC";\r
+  labels[kPosC]="PosC";\r
+\r
+  for(int i = 0;i< kTotalAC;++i){\r
+    fh2PhiPadRow[i] = new TH2F(Form("fh2PhiPadRow_%d",i),Form("Padrow vs phi AC %s;phi;padrow",labels[i].Data()),360,0.,360.,159,-0.5,158.5);\r
+    fh2PhiLayer[i] = new TH2F(Form("fh2PhiLayer_%d",i),Form("layer vs phi AC %s;phi;layer",labels[i].Data()),360,0.,360.,6,-0.5,5.5);\r
+  }\r
+  for(int i = 0;i<2;++i){\r
+    fh2MultSpdChips[i] = new TH2F(Form("fh2ChipsSpdMult_%d",i),"mult SPD vs chips;chips;Mult",201,-1.5,199.5,201,-1.5,199.5);\r
+  }\r
+  fh2VtxTpcSpd  = new TH2F("fh2VtxTpcSpd","SPD vtx vs TPC;tpc-z;spd-z",200,-20.,20.,200,-20,20.);\r
+\r
+  TString labelsVtx[kTotalVtx];\r
+  labelsVtx[kTPC] = "TPC";\r
+  labelsVtx[kSPD] = "SPD";\r
+  for(int i = 0; i < kTotalVtx; ++i){ \r
+    fh2VertexCorrelation[i] = new TH2F(Form("fh2VertexCorrelation_%d",i),Form("VertexCorrelation %s;MC z-vtx;ESD z-vtx",labelsVtx[i].Data()), 120, -30, 30, 120, -30, 30);\r
+    fh2VertexCorrelationShift[i] = new TH2F(Form("fh2VertexCorrelationShift_%d",i), Form("VertexCorrelationShift %s;MC z-vtx;MC z-vtx - ESD z-vtx",labelsVtx[i].Data()), 120, -30, 30, 100, -1, 1); \r
+    fh1VertexShift[i] =  new TH1F(Form("fh1VertexShift_%d",i), Form("VertexShift %s;(MC z-vtx - ESD z-vtx);Entries",labelsVtx[i].Data()), 201, -2, 2);       \r
+    fh1VertexShiftNorm[i] = new TH1F(Form("fh1VertexShiftNorm_%d",i), Form("VertexShiftNorm %s;(MC z-vtx - ESD z-vtx) / #sigma_{ESD z-vtx};Entries",labelsVtx[i].Data()), 200, -100, 100);      \r
+  }\r
+\r
+    \r
+  fh2EtaCorrelation = new TH2F("fh2EtaCorrelation", "EtaCorrelation;MC #eta;ESD #eta", 120, -3, 3, 120, -3, 3);\r
+  fh2EtaCorrelationShift = new TH2F("fh2EtaCorrelationShift", "EtaCorrelationShift;ESD #eta;MC #eta - ESD #eta", 120, -3, 3, 100, -0.1, 0.1);\r
+  fh2PhiCorrelation = new TH2F("fh2PhiCorrelation", "PhiCorrelation;MC #phi;ESD #phi", 120, 0, 360, 120, 0, 360);\r
+  fh2PhiCorrelationShift = new TH2F("fh2PhiCorrelationShift", "PhiCorrelationShift;ESD #phi;MC #phi - ESD #phi", 120, 0, 360, 100, -10, 10);\r
+  fh2PtCorrelation = new TH2F("fh2PtCorrelation", "PtCorrelation;MC p_T;ESD p_T", 200, 0., 200., 200, 0, 200);\r
+  fh2PtCorrelationShift = new TH2F("fh2PtCorrelationShift", "PtCorrelationShift;ESD p_T;MC p_T - ESD #p_T", 120, 0, 10, 100, -2, 2);\r
+  \r
+  fHistDeltaphiprimaries  =new TH1F("fHistDeltaphiprimaries", " #Delta#phi distribution: primaries",     314, -0.2,0.2);               \r
+  fHistDeltaphisecondaries=new TH1F("fHistDeltaphisecondaries", "#Delta#phi distribution: secondaries", 314, -0.2,0.2);              \r
+  fHistDeltaphireject     =new TH1F("fHistDeltaphireject", " #Delta#phi distribution: reject trackled",  314, -0.2,0.2); \r
+  \r
+  fHistTPCITSdeltax =new TH1F("fHistTPCITSdeltax", "TPC-ITS matching:#Delta x", 100,-20,20); \r
+  fHistTPCITSdeltay =new TH1F("fHistTPCITSdeltay", "TPC-ITS matching:#Delta y", 100,-20,20); \r
+  fHistTPCITSdeltaz =new TH1F("fHistTPCITSdeltaz", "TPC-ITS matching:#Delta z", 100,-20,20); \r
+  fHistTPCITSdeltar =new TH1F("fHistTPCITSdeltar", "TPC-ITS matching:#Delta r", 100,0,20); \r
+  \r
+  fHistTPCTRDdeltax   = new TH1F("fHistTPCTRDdeltax", "TPC-TRD matching:#Delta x", 400,-400, 400); \r
+  fHistTPCTRDdeltay   = new TH1F("fHistTPCTRDdeltay", "TPC-TRD matching:#Delta y", 400,-400, 400); \r
+  fHistTPCTRDdeltaz   = new TH1F("fHistTPCTRDdeltaz", "TPC-TRD matching:#Delta z", 400,-400, 400); \r
+  fHistTPCTRDdeltar   = new TH1F("fHistTPCTRDdeltar", "TPC-TRD matching:#Delta r", 100,0,20); \r
+  fHistTPCTRDdeltaphi = new TH1F("fHistTPCTRDdeltarphi", "TPC-TRD matching:#Delta #phi", 128,-3.14 , 3.14); \r
+  // Pulls\r
+  fPtPullsPos = new TProfile("fPtPullsPos", "Pt Pulls for pos. primaries", 50 , 0., 10., -5., 5.,"S");\r
+  fPtPullsNeg = new TProfile("fPtPullsNeg", "Pt Pulls for neg. primaries", 50 , 0., 10., -5., 5.,"S");\r
+  fPtShiftPos = new TProfile("fPtShiftPos", "Pt Shift for pos. primaries", 50 , 0., 10., -0.2, 0.2,"S");\r
+  fPtShiftNeg = new TProfile("fPtShiftNeg", "Pt Shift for neg. primaries", 50 , 0., 10., -0.2, 0.2,"S");\r
+\r
+  fEtaMultiFluc = new TProfile("fEtaMultiFluc", "eta multiplicity fluctuations", 10, -2., 2., 0., 600., "S");\r
+  fEtaMulti     = new TH1F("fEtaMulti",     "eta multiplicity fluctuations", 10, -2., 2.);\r
+  fEtaMultiH    = new TH1F("fEtaMultiH",    "eta multiplicity fluctuations", 600,  0., 600.);\r
+\r
+  fPhiMultiFluc = new TProfile("fPhiMultiFluc", "phi multiplicity fluctuations", 64, 0., 6.4, 0., 100., "S");\r
+  fPhiMulti     = new TH1F("fPhiMulti",         "phi multiplicity fluctuations", 64, 0., 6.4);\r
+  fPhiMultiH    = new TH1F("fPhiMultiH",        "phi multiplicity fluctuations", 100, 0., 100.);\r
+\r
+  // SPD Vertex\r
+  fVertexRvsZ  = new TH2F("fVertexRvsZ", "SPD Vertex R vs Z", 200, -1., 1., 200, -1., 1.);\r
+  fVertexRvsZC = new TH2F("fVertexRvsZC", "SPD Vertex R vs Z", 200, -1., 1., 200, -1., 1.);\r
+  \r
+  TString charge[8]={"Deu", "AntiDeu", "Tri", "AntiTri", "He3", "AntiHe3", "He4", "AntiHe4"};\r
+  for(Int_t i=0;i<8;i++){\r
+    fHistRECptCharge[i]   = new TH1F(Form("fHistRECptCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all reconstructed",\r
+                                    100, 0., 20.);\r
+    fHistMCptCharge[i]    = new TH1F(Form("fHistMCptCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all MC",\r
+                                    100, 0., 20.);\r
+    fHistMCNRptCharge[i]  = new TH1F(Form("fHistMCNRptCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all not-reconstructed MC", \r
+                                    100, 0., 20.);\r
+    fHistFAKEptCharge[i]  = new TH1F( Form("fHistFAKEptCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all Fake",                 \r
+                                    100, 0., 20.);\r
+\r
+    fHistRECetaCharge[i]   = new TH1F(Form("fHistRECetaCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all reconstructed",\r
+                                    100, 0., 20.);\r
+    fHistMCetaCharge[i]    = new TH1F(Form("fHistMCetaCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all MC",\r
+                                    100, 0., 20.);\r
+    fHistMCNRetaCharge[i]  = new TH1F(Form("fHistMCNRetaCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all not-reconstructed MC", \r
+                                    100, 0., 20.);\r
+    fHistFAKEetaCharge[i]  = new TH1F( Form("fHistFAKEetaCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all Fake",                 \r
+                                    100, 0., 20.);\r
+\r
+    fHistRECphiCharge[i]   = new TH1F(Form("fHistRECphiCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all reconstructed",\r
+                                    100, 0., 20.);\r
+    fHistMCphiCharge[i]    = new TH1F(Form("fHistMCphiCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all MC",\r
+                                    100, 0., 20.);\r
+    fHistMCNRphiCharge[i]  = new TH1F(Form("fHistMCNRphiCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all not-reconstructed MC", \r
+                                    100, 0., 20.);\r
+    fHistFAKEphiCharge[i]  = new TH1F( Form("fHistFAKEphiCharge%s",charge[i].Data()), \r
+                                    "p_{T} distribution: all Fake",                 \r
+                                    100, 0., 20.);\r
+\r
+  }\r
+  // BKC\r
+\r
+\r
+  fHists->SetOwner();\r
+\r
+  fHists->Add(fHistRECpt);\r
+  fHists->Add(fHistMCpt);\r
+  fHists->Add(fHistMCNRpt);\r
+  fHists->Add(fHistFAKEpt);\r
+  fHists->Add(fHistMULTpt);\r
+\r
+  fHists->Add(fHistRECeta);\r
+  fHists->Add(fHistMCeta);\r
+  fHists->Add(fHistMCNReta);\r
+  fHists->Add(fHistFAKEeta);\r
+  fHists->Add(fHistMULTeta);\r
+\r
+  fHists->Add(fHistRECphi);\r
+  fHists->Add(fHistMCphi);\r
+  fHists->Add(fHistMCNRphi);\r
+  fHists->Add(fHistFAKEphi);\r
+  fHists->Add(fHistMULTphi);\r
+\r
+  fHists->Add(fHistRECHPTeta);\r
+  fHists->Add(fHistMCHPTeta);\r
+  fHists->Add(fHistMCNRHPTeta);\r
+  fHists->Add(fHistFAKEHPTeta);\r
+\r
+  fHists->Add(fHistRECHPTphi);\r
+  fHists->Add(fHistMCHPTphi);\r
+  fHists->Add(fHistMCNRHPTphi);\r
+  fHists->Add(fHistFAKEHPTphi);\r
+\r
+  // CKB\r
+\r
+  fHists->Add(fh1VtxEff);\r
+  for(int i = 0;i < kTotalAC;++i){\r
+    fHists->Add(fh2PhiPadRow[i]);\r
+    fHists->Add(fh2PhiLayer[i]);\r
+  }\r
+  for(int i = 0;i<2;++i){\r
+    fHists->Add(fh2MultSpdChips[i]);\r
+  }\r
+  fHists->Add(fh2VtxTpcSpd);\r
+    \r
+  for(int i = 0; i < kTotalVtx; ++i){ \r
+    fHists->Add(fh2VertexCorrelation[i]);\r
+    fHists->Add(fh2VertexCorrelationShift[i]);\r
+    fHists->Add(fh1VertexShift[i]);      \r
+    fHists->Add(fh1VertexShiftNorm[i]);\r
+  }\r
+\r
+\r
+  fHists->Add(fh2EtaCorrelation);\r
+  fHists->Add(fh2EtaCorrelationShift);\r
+  fHists->Add(fh2PhiCorrelation);\r
+  fHists->Add(fh2PhiCorrelationShift);\r
+  fHists->Add(fh2PtCorrelation);\r
+  fHists->Add(fh2PtCorrelationShift);\r
+\r
+  fHists->Add(fHistDeltaphiprimaries);              \r
+  fHists->Add(fHistDeltaphisecondaries);\r
+  fHists->Add(fHistDeltaphireject);  \r
+\r
+  fHists->Add(fHistTPCITSdeltax);  \r
+  fHists->Add(fHistTPCITSdeltay);  \r
+  fHists->Add(fHistTPCITSdeltaz);  \r
+  fHists->Add(fHistTPCITSdeltar);  \r
+    \r
+  fHists->Add(fHistTPCTRDdeltax);  \r
+  fHists->Add(fHistTPCTRDdeltay);  \r
+  fHists->Add(fHistTPCTRDdeltaz);  \r
+  fHists->Add(fHistTPCTRDdeltar);  \r
+  fHists->Add(fHistTPCTRDdeltaphi);  \r
+\r
+  fHists->Add(fHistRecMult);\r
+  fHists->Add(fHistNCluster);\r
\r
+\r
+  fHists->Add(fPtPullsPos);\r
+  fHists->Add(fPtPullsNeg);\r
+  fHists->Add(fPtShiftPos);\r
+  fHists->Add(fPtShiftNeg);\r
+\r
+  fHists->Add(fVertexRvsZ);\r
+  fHists->Add(fVertexRvsZC);\r
+  fHists->Add(fEtaMultiFluc);\r
+  fHists->Add(fEtaMulti);\r
+  fHists->Add(fEtaMultiH);\r
+  fHists->Add(fPhiMultiFluc);\r
+  fHists->Add(fPhiMulti);\r
+  fHists->Add(fPhiMultiH);\r
+\r
+  for(Int_t i=0;i<8;i++){\r
+    fHists->Add(fHistRECptCharge[i]);\r
+    fHists->Add(fHistMCptCharge[i]);\r
+    fHists->Add(fHistMCNRptCharge[i]);\r
+    fHists->Add(fHistFAKEptCharge[i]);\r
+\r
+    fHists->Add(fHistRECetaCharge[i]);\r
+    fHists->Add(fHistMCetaCharge[i]);\r
+    fHists->Add(fHistMCNRetaCharge[i]);\r
+    fHists->Add(fHistFAKEetaCharge[i]);\r
+\r
+    fHists->Add(fHistRECphiCharge[i]);\r
+    fHists->Add(fHistMCphiCharge[i]);\r
+    fHists->Add(fHistMCNRphiCharge[i]);\r
+    fHists->Add(fHistFAKEphiCharge[i]);\r
+  }\r
+\r
+  for (Int_t i=0; i<fHists->GetEntries(); ++i) {\r
+    TH1 *h1 = dynamic_cast<TH1*>(fHists->At(i));\r
+    if (h1){\r
+      h1->Sumw2();\r
+    }\r
+  }\r
+  // BKC\r
+\r
+\r
+  // Post output data.\r
+  PostData(1, fHists);\r
+\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskEfficiency::UserExec(Option_t *) \r
+{\r
+  // Event loop\r
+  // Analysis of MC true particles and reconstructed tracks\r
+  // Different track types (Global, TPC, ITS) can be selected\r
+\r
+  if (!fInputEvent) {\r
+    Printf("ERROR: fESD not available");\r
+    return;\r
+  }\r
+  static Int_t nfc = 0;\r
+\r
+  //AliESDInputHandler* esdI = (AliESDInputHandler*) \r
+  //(AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler();\r
+  //AliESDEvent* hltEvent = esdI->GetHLTEvent();\r
+\r
+  AliStack* stack = MCEvent()->Stack();\r
+  AliESDEvent* esdE = (AliESDEvent*) fInputEvent;\r
+    \r
+  // CKB\r
+  // Fetch the SPD vertex:\r
+  Double_t vtxSPD[3];\r
+  const AliESDVertex* spdVertex = esdE->GetPrimaryVertexSPD();\r
+  if (spdVertex) {\r
+    if(spdVertex->GetNContributors()<=0){\r
+      spdVertex = 0;\r
+    }\r
+    else{\r
+      spdVertex->GetXYZ(vtxSPD);\r
+    }\r
+  }\r
+  \r
+  // Fetch the TPC vertex\r
+  Double_t vtxTPC[3];\r
+  const AliESDVertex* tpcVertex = esdE->GetPrimaryVertexTPC();\r
+  if (tpcVertex) {\r
+    if(tpcVertex->GetNContributors()<=0){\r
+      tpcVertex = 0;\r
+    }\r
+    else{\r
+      tpcVertex->GetXYZ(vtxTPC);\r
+    }\r
+  }\r
+  // SPD Vertex\r
+  if (spdVertex) {\r
+    Double_t x,y,z;\r
+    x = spdVertex->GetX() + 0.07;\r
+    y = spdVertex->GetY() - 0.25;\r
+    z = spdVertex->GetZ();\r
+    if (TMath::Abs(z) < 10.) {\r
+      fVertexRvsZ->Fill(x,y);\r
+      if (TMath::Sqrt(x*x + y*y) > 5. * 0.028) \r
+       fVertexRvsZC->Fill(x,y);\r
+    }\r
+  }\r
+  // BKC\r
+  //Printf("%s:%d %5d",(char*)__FILE__,__LINE__, esdE->GetNumberOfTracks());\r
+  TArrayI labels(esdE->GetNumberOfTracks());\r
+  Int_t igood = 0;\r
+  // Track loop to fill a pT spectrum\r
+  //printf("ESD track loop \n");\r
+\r
+  AliESDtrack *tpcP = 0x0;\r
+\r
+  for (Int_t iTracks = 0; iTracks < esdE->GetNumberOfTracks(); iTracks++) {\r
+\r
+    //prevent mem leak for TPConly track\r
+    if(fTrackType==2&&tpcP){\r
+      delete tpcP;\r
+      tpcP = 0;\r
+    }\r
+\r
+    AliVParticle *track = esdE->GetTrack(iTracks);\r
+    AliESDtrack *esdtrack =  dynamic_cast<AliESDtrack*>(track);\r
+    esdtrack->PropagateToDCA(esdE->GetPrimaryVertex(),\r
+                            esdE->GetMagneticField(), 10000.);\r
+\r
+    if (!track) {\r
+      Printf("ERROR: Could not receive track %d", iTracks);\r
+      continue;\r
+    }\r
+    //__________\r
+    // run Task for global tracks or ITS tracks or TPC tracks\r
+    const AliExternalTrackParam *tpcPin = 0x0;\r
+    //Double_t phiIn=0.;\r
+    Float_t phi = 0.;\r
+    Float_t eta = 0.;\r
+    Float_t pt  = 0.;\r
+    Int_t indexAC = GetIndexAC(esdtrack);  \r
+\r
+    //select in the steering macro, which track type should be analysed. \r
+    // Four track types are selectable:\r
+    // 0. Global tracks\r
+    // 1. ITS tracks (SA or Pure SA)\r
+    // 2. TPC only tracks\r
+    // Track selection copied from PWG1/AliAnalysisTaskQAsym\r
+\r
+    if(fTrackType==0){\r
+      //Fill all histograms with global tracks\r
+      tpcP = esdtrack;\r
+      if (!tpcP) continue;\r
+      if (!fCuts->AcceptTrack(tpcP)) continue;\r
+      phi = tpcP->Phi();\r
+      eta = tpcP->Eta();\r
+      pt  = tpcP->Pt();\r
+    }\r
+    else if(fTrackType==1){\r
+      //Fill all histograms with ITS tracks\r
+      tpcP = esdtrack;\r
+      if (!tpcP) continue;\r
+      if (!fCuts->AcceptTrack(tpcP)) continue;\r
+      phi = tpcP->Phi();\r
+      eta = tpcP->Eta();\r
+      pt  = tpcP->Pt();\r
+    }\r
+    else if(fTrackType==2){     \r
+      //Fill all histograms with TPC track information\r
+      tpcPin = esdtrack->GetInnerParam();\r
+      if (!tpcPin) continue;\r
+      tpcP = AliESDtrackCuts::\r
+       GetTPCOnlyTrack(dynamic_cast<AliESDEvent*>(esdE),\r
+                       esdtrack->GetID());\r
+      if (!tpcP) continue;\r
+      if (!fCuts->AcceptTrack(tpcP)) continue;\r
+      if(tpcP->GetNcls(1)>160)continue;//jacek's track cut\r
+      if(tpcP->GetConstrainedChi2TPC()<0)continue; // jacek's track cut\r
+      phi=tpcPin->Phi();\r
+      eta=tpcPin->Eta();\r
+      pt=tpcPin->Pt();\r
+    }\r
+    else{\r
+      Printf("ERROR: wrong track type \n");\r
+      continue;\r
+    }\r
+    //___________\r
+    //\r
+\r
+\r
+    labels.AddAt(-1, iTracks);\r
+      \r
+\r
+      // Selected\r
+      if (TMath::Abs(eta) > 0.9) {\r
+      } else {\r
+\r
+       //Int_t charge=track->Charge();\r
+       // Within acceptance\r
+       Int_t ind = TMath::Abs(esdtrack->GetLabel());\r
+       AliMCParticle* mcPtest = (AliMCParticle*) MCEvent()->GetTrack(ind);\r
+       if (stack->IsPhysicalPrimary(mcPtest->Label())){\r
+         if(TMath::Abs(mcPtest->PdgCode())>=1000020030){//all helium (+-1000020030,+-1000020040)\r
+           pt*=2;// reconstruction takes charge=1 -> for particles with charge=2, pT,rec need to be corrected\r
+         }\r
+         \r
+         Int_t index=ConvertePDG(mcPtest->PdgCode());\r
+         if(fDebug>1)Printf("PDG=%d, index=%d", mcPtest->PdgCode(), index);\r
+         //fill tracks comming from d,t,3He, 4He\r
+         if(index<8){\r
+           fHistRECptCharge [index] ->Fill(pt);\r
+           fHistRECetaCharge[index]->Fill(eta);\r
+           fHistRECphiCharge[index]->Fill(phi);\r
+         }\r
+\r
+       }\r
+       fHistRECpt ->Fill(pt);\r
+       fHistRECeta->Fill(eta);\r
+       fHistRECphi->Fill(phi);\r
+\r
+       if (pt > 2.) {\r
+         fHistRECHPTeta->Fill(eta);\r
+         fHistRECHPTphi->Fill(phi);\r
+       }\r
+         \r
+       //      Int_t ind = TMath::Abs(esdtrack->GetLabel());\r
+       AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(ind);\r
+       if (!(stack->IsPhysicalPrimary(mcP->Label()))) {\r
+         fHistFAKEpt ->Fill(pt);\r
+         fHistFAKEeta->Fill(eta);                \r
+         fHistFAKEphi->Fill(phi);\r
+         \r
+         if (pt > 2.) {\r
+           fHistFAKEHPTeta->Fill(eta);\r
+           fHistFAKEHPTphi->Fill(phi);\r
+         }\r
+           \r
+       }\r
+         \r
+       labels.AddAt(TMath::Abs(esdtrack->GetLabel()), iTracks);\r
+       igood++;\r
+         \r
+         \r
+       Float_t phiDeg = TMath::RadToDeg()*track->Phi();\r
+         \r
+       //TPC-ITS matching\r
+       //TPC-TRD matching  \r
+             \r
+       if (tpcP){\r
+         Int_t labeltpcits = esdtrack->GetLabel();\r
+         AliExternalTrackParam * tpcPCopy = new AliExternalTrackParam(*tpcP);\r
+         Double_t xk = 43.6;  // cm\r
+         Double_t bz = 5.0;   // kG\r
+         if(tpcPCopy->PropagateTo(xk,bz)){\r
+           Double_t alpha=tpcPCopy->GetAlpha();\r
+           // if(tpcPCopy->Rotate(0.)){ \r
+           Float_t xtpc=tpcPCopy->GetX();\r
+           Float_t ytpc=tpcPCopy->GetY();\r
+           Float_t ztpc=tpcPCopy->GetZ();\r
+           Float_t xpr,ypr ; \r
+           xpr=xtpc*TMath::Cos(alpha)-ytpc*TMath::Sin(alpha);\r
+           ypr=xtpc*TMath::Sin(alpha)+ytpc*TMath::Cos(alpha);\r
+               \r
+           AliMCParticle* mcPart = (AliMCParticle*) MCEvent()->GetTrack(abs(labeltpcits));\r
+           Int_t ntref = mcPart->GetNumberOfTrackReferences();\r
+                \r
+           for (Int_t k = 0; k < ntref; k++) {\r
+             AliTrackReference* ref = mcPart->GetTrackReference(k);\r
+             Float_t xhits = ref->X();\r
+             Float_t yhits = ref->Y();\r
+             Float_t radio = TMath::Sqrt(xhits*xhits+yhits*yhits);\r
+               \r
+             if(ref->DetectorId() == 0 && (radio > 42.8 && radio < 43.7)) {\r
+                 \r
+               Float_t xits = ref->X();\r
+               Float_t yits = ref->Y();\r
+               Float_t zits = ref->Z();\r
+                 \r
+               Float_t difx=(xits-xpr);\r
+               fHistTPCITSdeltax->Fill(difx);\r
+               Float_t dify=(yits-ypr);\r
+               fHistTPCITSdeltay->Fill(dify);\r
+               Float_t difz=(zits-ztpc);\r
+               fHistTPCITSdeltaz->Fill(difz);\r
+               Float_t deltar = TMath::Sqrt(difx * difx + dify * dify + difz * difz);\r
+               fHistTPCITSdeltar->Fill(deltar);\r
+               break;\r
+             }\r
+             // }\r
+           } // trackrefs\r
+         } \r
+       } //ITS-TPC maching\r
+         //TPC-TRD maching \r
+              \r
+       const AliExternalTrackParam *trd = esdtrack->GetOuterParam();\r
+         \r
+       if (trd){Int_t labeltpctrd = track->GetLabel();\r
+         \r
+         AliExternalTrackParam * trdCopy = new AliExternalTrackParam(*trd);\r
+         Double_t xktrd=296.0;\r
+         Double_t bztrd=5.0;\r
+         if(trdCopy->PropagateTo(xktrd,bztrd)){\r
+           Float_t xtpc2=trdCopy->GetX();\r
+           Float_t ytpc2=trdCopy->GetY();\r
+           Float_t ztpc2=trdCopy->GetZ();\r
+           Double_t alpha=trdCopy->GetAlpha();\r
+           Float_t xpr,ypr ; \r
+           xpr=xtpc2*TMath::Cos(alpha)-ytpc2*TMath::Sin(alpha);\r
+           ypr=xtpc2*TMath::Sin(alpha)+ytpc2*TMath::Cos(alpha);\r
+           AliMCParticle* mcPart = (AliMCParticle*) MCEvent()->GetTrack(abs(labeltpctrd));\r
+           Int_t ntreftrd = mcPart->GetNumberOfTrackReferences(); \r
+           \r
+           for (Int_t k = 0; k < ntreftrd; k++) {\r
+             AliTrackReference* ref = mcPart->GetTrackReference(k);\r
+             if(ref->DetectorId() == 3 && ref->Label() == labeltpctrd){\r
+               \r
+               Float_t xtrd = ref->X();\r
+               Float_t ytrd = ref->Y();\r
+               Float_t ztrd = ref->Z();\r
+               \r
+               Float_t difx=(xtrd-xpr);\r
+               fHistTPCTRDdeltax->Fill(difx);\r
+               Float_t dify=(ytrd-ypr);\r
+               fHistTPCTRDdeltay->Fill(dify);\r
+               Float_t difz=(ztrd-ztpc2);\r
+               fHistTPCTRDdeltaz->Fill(difz);\r
+               Float_t deltar = TMath::Sqrt(difx * difx + dify * dify + difz * difz);\r
+               fHistTPCTRDdeltar->Fill(deltar);\r
+               Float_t phi_tpc = TMath::ATan2(ypr, xpr);\r
+               Float_t phi_trd = TMath::ATan2(ytrd, xtrd);\r
+               Float_t dphi = phi_trd - phi_tpc;\r
+               if (dphi >   TMath::Pi()) dphi -= 2. * TMath::Pi();\r
+               if (dphi < - TMath::Pi()) dphi += 2. * TMath::Pi();\r
+\r
+               fHistTPCTRDdeltaphi->Fill(dphi);\r
+               break;\r
+             }\r
+             \r
+           }\r
+         }\r
+\r
+       }//TRD-TPC maching\r
+\r
+       // CKB\r
+       if(pt>2.&&fFieldOn){\r
+         TBits bits(esdtrack->GetTPCClusterMap());\r
+         for(unsigned int ip = 0;ip < bits.GetNbits();++ip){\r
+           if(bits[ip]){\r
+             fh2PhiPadRow[indexAC]->Fill(phiDeg,ip);\r
+           }\r
+         }\r
+         for(int i = 0;i < 6;++i){\r
+           if(esdtrack->HasPointOnITSLayer(i))fh2PhiLayer[indexAC]->Fill(phiDeg,i);\r
+         }\r
+       }\r
+       else if(!fFieldOn){ // field not on all momenta are set to MPV \r
+         TBits bits(esdtrack->GetTPCClusterMap());\r
+         for(unsigned int ip = 0;ip < bits.GetNbits();++ip){\r
+           if(bits[ip]){\r
+             fh2PhiPadRow[indexAC]->Fill(phiDeg,ip);\r
+           }\r
+           for(int i = 0;i < 6;++i){\r
+             if(esdtrack->HasPointOnITSLayer(i))fh2PhiLayer[indexAC]->Fill(phiDeg,i);\r
+           }\r
+         }      \r
+\r
+       }\r
+       // Fill the correlation\r
+       if(MCEvent()){\r
+         TParticle *part = MCEvent()->Stack()->Particle(TMath::Abs(esdtrack->GetLabel()));\r
+         if(part){\r
+           Float_t mcEta = part->Eta();\r
+           Float_t mcPhi = TMath::RadToDeg()*part->Phi();\r
+           Float_t mcPt  = part->Pt();\r
+           //             if (pt - mcPt > 20.) {\r
+           fh2EtaCorrelation->Fill(mcEta,eta);\r
+           fh2EtaCorrelationShift->Fill(eta,mcEta-eta);\r
+           fh2PhiCorrelation->Fill(mcPhi,phiDeg);\r
+           fh2PhiCorrelationShift->Fill(phiDeg,mcPhi-phiDeg);\r
+           fh2PtCorrelation->Fill(mcPt,pt);\r
+           fh2PtCorrelationShift->Fill(pt,mcPt-pt);\r
+           //}\r
+           Double_t sigmaPt = TMath::Sqrt(esdtrack->GetSigma1Pt2());\r
+           if (track->Charge() > 0.) {\r
+             fPtPullsPos->Fill(mcPt, (1./pt - 1./mcPt) / sigmaPt); \r
+             fPtShiftPos->Fill(mcPt, (1./pt - 1./mcPt));\r
+           } else {\r
+             fPtPullsNeg->Fill(mcPt, (1./pt - 1./mcPt) / sigmaPt); \r
+             fPtShiftNeg->Fill(mcPt, (1./pt - 1./mcPt));  \r
+           }\r
+         }\r
+         // BKC\r
+       }\r
+      }\r
+   \r
+  } //track loop \r
+\r
+  //prevent mem leak for TPConly track\r
+  if(fTrackType==2&&tpcP){\r
+    delete tpcP;\r
+    tpcP = 0;\r
+  }\r
+\r
+  //Printf("%s:%d",(char*)__FILE__,__LINE__);\r
+  // CKB\r
+  // Vertex eff\r
+  if(tpcVertex){\r
+    fh1VtxEff->Fill("TPC VTX",1);\r
+  }\r
+  else{\r
+    fh1VtxEff->Fill("NO TPC VTX",1);\r
+  }\r
+\r
+  if(spdVertex){\r
+    fh1VtxEff->Fill("SPD VTX",1);\r
+  }\r
+  else{\r
+    fh1VtxEff->Fill("NO SPD VTX",1);\r
+  }\r
+\r
+\r
+  // Vertex correlation SPD vs. TPC\r
+  if(tpcVertex&&spdVertex){\r
+    fh2VtxTpcSpd->Fill(vtxTPC[2],vtxSPD[2]);\r
+  }\r
+  //  Printf("%s:%d",(char*)__FILE__,__LINE__);\r
+  // Multiplicity checks in the SPD\r
+  const AliMultiplicity *spdMult = esdE->GetMultiplicity();\r
+  // Multiplicity Analysis\r
+  Int_t nt = spdMult->GetNumberOfTracklets();\r
+  nfc++;\r
+  if (nfc == 520) {\r
+    nfc = 0;\r
+    for (Int_t ib = 1; ib <= 10; ib++) {\r
+      Int_t ic = Int_t(fEtaMulti->GetBinContent(ib));\r
+      Float_t eta = -1.8 + Float_t(ib-1) * 0.4;\r
+      fEtaMultiFluc->Fill(eta, Float_t(ic));\r
+      if (ib == 5) fEtaMultiH->Fill(Float_t(ic));\r
+    }\r
+\r
+    for (Int_t ib = 1; ib <= 64; ib++) {\r
+      Int_t ic = Int_t(fPhiMulti->GetBinContent(ib));\r
+      Float_t phi = 0.05 + Float_t(ib-1) * 0.1;\r
+      fPhiMultiFluc->Fill(phi, Float_t(ic));\r
+      if (ib  == 2) fPhiMultiH->Fill(Float_t(ic));\r
+    }\r
+\r
+    fEtaMulti->Reset();\r
+    fPhiMulti->Reset();    \r
+  }\r
+       \r
+  for (Int_t j = 0; j < nt; j++) {\r
+    fEtaMulti->Fill(spdMult->GetEta(j));\r
+    fPhiMulti->Fill(spdMult->GetPhi(j));\r
+  }\r
+\r
+  Short_t chips0 = spdMult->GetNumberOfFiredChips(0);\r
+  Short_t chips1 = spdMult->GetNumberOfFiredChips(1);\r
+  //Printf("%s:%d",(char*)__FILE__,__LINE__);\r
+  Int_t inputCount = 0;\r
+  \r
+  if(spdVertex){\r
+    // get multiplicity from ITS tracklets\r
+    for (Int_t i=0; i<spdMult->GetNumberOfTracklets(); ++i){\r
+      Float_t deltaPhi = spdMult->GetDeltaPhi(i);\r
+      // prevent values to be shifted by 2 Pi()\r
+      if (deltaPhi < -TMath::Pi())\r
+       deltaPhi += TMath::Pi() * 2;\r
+      if (deltaPhi > TMath::Pi())\r
+       deltaPhi -= TMath::Pi() * 2;      \r
+      if (TMath::Abs(deltaPhi) > 1)\r
+       printf("WARNING: Very high Delta Phi: %d %f %f %f\n", i, spdMult->GetTheta(i), spdMult->GetPhi(i), deltaPhi);\r
+      //trackled Deltaphi\r
+      Int_t label1=spdMult->GetLabel(i,0);\r
+      Int_t label2=spdMult->GetLabel(i,1);\r
+      if ( label1==label2){\r
+       if(stack->IsPhysicalPrimary(label1) == kTRUE)\r
+         fHistDeltaphiprimaries->Fill(deltaPhi);\r
+       else\r
+         fHistDeltaphisecondaries->Fill(deltaPhi);\r
+      }\r
+      else{\r
+       fHistDeltaphireject->Fill(deltaPhi);\r
+      }\r
+      ++inputCount;\r
+    }\r
+    // Printf("%s:%d",(char*)__FILE__,__LINE__);\r
+    fh2MultSpdChips[0]->Fill(chips0,inputCount);\r
+    fh2MultSpdChips[1]->Fill(chips1,inputCount);\r
+    //Printf("%s:%d",(char*)__FILE__,__LINE__);\r
+  }\r
+  // BKC\r
+\r
+  // MC analysis\r
+  Int_t nmc = MCEvent()->GetNumberOfTracks();\r
+  AliGenEventHeader*       header = MCEvent()->GenEventHeader();\r
+  // some test\r
+  AliGenCocktailEventHeader* cH = dynamic_cast<AliGenCocktailEventHeader*> (header);\r
+  AliGenPythiaEventHeader* pH;\r
+  if (cH == 0) \r
+    {\r
+      header->Dump();\r
+      pH = dynamic_cast<AliGenPythiaEventHeader*>(header);\r
+    } else {\r
+    TObject* entry = cH->GetHeaders()->FindObject("Pythia");\r
+    pH = dynamic_cast<AliGenPythiaEventHeader*>(entry);\r
+  }\r
+  Int_t iproc = -1;\r
+\r
+  if (pH) iproc = pH->ProcessType();\r
+\r
+  TArrayF mcV;\r
+  header->PrimaryVertex(mcV);\r
+  Float_t vzMC = mcV[2];\r
+\r
+  // CKB\r
+  // Fill the correlation with MC\r
+  if(spdVertex&&MCEvent()){\r
+    Float_t diff = vzMC-vtxSPD[2];\r
+    fh2VertexCorrelation[kSPD]->Fill(vzMC,vtxSPD[2]);\r
+    fh2VertexCorrelationShift[kSPD]->Fill(vzMC,diff);\r
+    fh1VertexShift[kSPD]->Fill(diff);\r
+    if(spdVertex->GetZRes()>0)fh1VertexShiftNorm[kSPD]->Fill(diff/spdVertex->GetZRes());\r
+  }\r
+  if(tpcVertex&&MCEvent()){\r
+    Float_t diff = vzMC-vtxTPC[2];\r
+    fh2VertexCorrelation[kTPC]->Fill(vzMC,vtxTPC[2]);\r
+    fh2VertexCorrelationShift[kTPC]->Fill(vzMC,diff);\r
+    fh1VertexShift[kTPC]->Fill(diff);\r
+    if(tpcVertex->GetZRes()>0)fh1VertexShiftNorm[kTPC]->Fill(diff/tpcVertex->GetZRes());\r
+  }\r
+  // BKC\r
+\r
+\r
+  // MC event loop\r
+  //printf("MC particle loop \n");\r
+  for (Int_t i = 0; i < nmc; i++) {\r
+    AliMCParticle* mcP = (AliMCParticle*) MCEvent()->GetTrack(i);\r
+    //printf("MC particle loop %5d \n", i);\r
+    // Primaries only\r
+    if (!(stack->IsPhysicalPrimary(mcP->Label()))) continue;\r
+    if (mcP->Particle()->GetPDG()->Charge() == 0) continue;\r
+    // Int_t charge=Int_t(mcP->Particle()->GetPDG()->Charge());\r
+    Float_t phi = mcP->Phi();\r
+    Float_t eta = mcP->Eta();\r
+    Float_t pt  = mcP->Pt();\r
+    if (TMath::Abs(eta) > 0.9) continue;\r
+    Int_t ntr = mcP->GetNumberOfTrackReferences();\r
+    Int_t nITS = 0;\r
+    Int_t nTPC = 0;\r
+    Int_t nFRA = 0;\r
+    Float_t  x = 0.;\r
+    Float_t  y = 0.;\r
+      \r
+    for (Int_t j = 0; j < ntr; j++) {\r
+         \r
+      AliTrackReference* ref = mcP->GetTrackReference(j);\r
+       \r
+      if (ref->DetectorId() == 0) nITS++;\r
+      if (ref->DetectorId() == 1) nTPC++;\r
+      if (ref->DetectorId() == 2) nFRA++;\r
+      if (nTPC == 1) {\r
+       x = ref->X();\r
+       y = ref->Y();\r
+       break;\r
+      }\r
+    }\r
+\r
+    fHistMCpt ->Fill(pt);\r
+    fHistMCeta->Fill(eta);\r
+    fHistMCphi->Fill(phi);\r
+    Int_t index=ConvertePDG(mcP->PdgCode());\r
+    if(index<8){\r
+      fHistMCptCharge [index] ->Fill(pt);\r
+      fHistMCetaCharge[index]->Fill(eta);\r
+      fHistMCphiCharge[index]->Fill(phi);\r
+    }\r
+\r
+\r
+    if (pt > 2.) {\r
+      fHistMCHPTeta->Fill(eta);\r
+      fHistMCHPTphi->Fill(phi);\r
+    }\r
+\r
+    Bool_t reco = kFALSE;\r
+    Int_t multR =  0;\r
+    Int_t jold  = -1;\r
+      \r
+    for (Int_t j = 0; j < esdE->GetNumberOfTracks(); j++) {\r
+      if (i == labels.At(j)) {\r
+       reco = kTRUE;\r
+       multR++;\r
+       //AliESDtrack* test = esdE->GetTrack(j);\r
+       if (multR > 1) {\r
+         Int_t nclus = 0;\r
+         AliESDtrack* track = esdE->GetTrack(jold);\r
+         nclus = track->GetTPCclusters(0);              \r
+         fHistNCluster->Fill(nclus);\r
+         \r
+                 \r
+         track = esdE->GetTrack(j);\r
+         nclus = track->GetTPCclusters(0);\r
+         fHistNCluster->Fill(nclus);\r
+         \r
+\r
+       }\r
+       jold = j;\r
+      }\r
+    }\r
+\r
+    fHistRecMult->Fill(Float_t(multR), 1.);\r
+    if (multR == 0) {\r
+      fHistMCNRpt ->Fill(pt);\r
+      fHistMCNReta->Fill(eta);\r
+      fHistMCNRphi->Fill(phi);\r
+      if(index<8){\r
+       fHistMCNRptCharge [index] ->Fill(pt);\r
+       fHistMCNRetaCharge[index]->Fill(eta);\r
+       fHistMCNRphiCharge[index]->Fill(phi);\r
+      }\r
+     \r
+\r
+      if (pt > 2.) {\r
+       fHistMCNRHPTeta->Fill(eta);\r
+       fHistMCNRHPTphi->Fill(phi);\r
+      }\r
+    }\r
+    else if (multR > 1)\r
+      {\r
+       fHistMULTpt ->Fill(pt);\r
+       fHistMULTeta->Fill(eta);\r
+       fHistMULTphi->Fill(phi);\r
+      }        \r
+  } // MC particle loop\r
+  //printf("End of MC particle loop \n");\r
+  \r
+  \r
+}      \r
+\r
+//________________________________________________________________________\r
+void AliAnalysisTaskEfficiency::Terminate(Option_t *) \r
+{\r
+\r
+}  \r
+\r
+\r
+Bool_t AliAnalysisTaskEfficiency::SelectJouri(AliESDtrack* track) \r
+{\r
+\r
+  Bool_t selected = kTRUE;\r
+  AliESDEvent* esdE = (AliESDEvent*) fInputEvent;    \r
+  // > 50 Clusters\r
+  if (track->GetTPCclusters(0) < 50) selected = kFALSE;\r
+\r
+  const AliESDVertex *vtx=esdE->GetPrimaryVertexSPD();\r
+  if (!vtx->GetStatus()) selected = kFALSE;\r
+    \r
+  Double_t zv=vtx->GetZv();\r
+\r
+    \r
+  const AliExternalTrackParam *ct=track->GetTPCInnerParam();\r
+  if (!ct)  selected = kFALSE;\r
+    \r
+  Double_t xyz[3];\r
+  ct->GetXYZ(xyz);\r
+  Float_t rv = TMath::Sqrt(xyz[0] * xyz[0] + xyz[1] * xyz[1]);\r
+  if (rv > 3.0)                     selected = kFALSE;\r
+  if (TMath::Abs(xyz[2] - zv)>2.5) selected = kFALSE;\r
+  return (selected);\r
+    \r
+}\r
+\r
+Int_t AliAnalysisTaskEfficiency::GetIndexAC(AliESDtrack *track){\r
+  if(!track)return -1;\r
+\r
+  // Crude selection for A and C side\r
+  // just with eta\r
+  if (track->Eta() > 0) { \r
+    if (track->Charge() > 0) \r
+      return kPosA;\r
+    else\r
+      return kNegA;    \r
+  }\r
+  else { // C side\r
+    if (track->Charge() > 0) \r
+      return kPosC;\r
+    else\r
+      return kNegC;    \r
+  }\r
+  \r
+  return -1;\r
+\r
+}\r
+\r
+\r
+Int_t AliAnalysisTaskEfficiency::ConvertePDG(Int_t pdg)\r
+{\r
+\r
+  // function converts the pdg values for nuclei (d, t, 3He, \r
+  // 4He and their antiparticles) into indices for histo-array\r
+  // filled in UserExec\r
+  \r
+  if      ( pdg ==  1000010020 )return 0;\r
+  else if ( pdg == -1000010020 )return 1;\r
+  else if ( pdg ==  1000010030 )return 2;\r
+  else if ( pdg == -1000010030 )return 3;\r
+  else if ( pdg ==  1000020030 )return 4;\r
+  else if ( pdg == -1000020030 )return 5;\r
+  else if ( pdg ==  1000020040 )return 6;\r
+  else if ( pdg == -1000020040 )return 7;\r
+  else return 9;\r
+  \r
+  \r
+}\r
diff --git a/test/vmctest/scripts/efficiency/AliAnalysisTaskEfficiency.h b/test/vmctest/scripts/efficiency/AliAnalysisTaskEfficiency.h
new file mode 100644 (file)
index 0000000..d5179e8
--- /dev/null
@@ -0,0 +1,149 @@
+#ifndef AliAnalysisTaskEfficiency_cxx\r
+#define AliAnalysisTaskEfficiency_cxx\r
\r
+class TH1F;\r
+class TH2F;\r
+class TList;\r
+class TProfile;\r
+\r
+class AliESDEvent;\r
+class AliESDtrack;\r
+class AliESDtrackCuts;\r
+\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class AliAnalysisTaskEfficiency : public AliAnalysisTaskSE {\r
+ public:\r
+  AliAnalysisTaskEfficiency(const char *name = "AliAnalysisTaskEfficiency");\r
+  virtual ~AliAnalysisTaskEfficiency() {}\r
+  \r
+  virtual void   UserCreateOutputObjects();\r
+  virtual void   UserExec(Option_t *option);\r
+  virtual void   Terminate(Option_t *);\r
+  virtual Bool_t SelectJouri(AliESDtrack* track);\r
+  virtual void   SetTrackType(Int_t type)      {fTrackType = type;} \r
+  virtual void   SetCuts(AliESDtrackCuts* cuts){fCuts = cuts;}\r
+  virtual void   SetFieldOn(Bool_t b = kTRUE)  {fFieldOn = b;} \r
+  \r
+ private:\r
+\r
+  // For crude division into A/C side and positive/negative\r
+  enum { kNegA = 0, kPosA, kNegC, kPosC, kTotalAC};\r
+  enum { kTPC = 0, kSPD, kTotalVtx};\r
+  Int_t GetIndexAC(AliESDtrack *track);\r
+  Int_t ConvertePDG(Int_t pdg);\r
+\r
+  Bool_t      fFieldOn;\r
+\r
+  TList       *fHists;          // List of histos\r
+  TH1F        *fHistRECpt;      // pt of reconstructed tracks\r
+  TH1F        *fHistMCpt;       // pt of MC primaries\r
+  TH1F        *fHistMCNRpt;     // pt of not reconstructed primaries\r
+  TH1F        *fHistFAKEpt;     // pt of fake track\r
+  TH1F        *fHistMULTpt;     // pt of multiply reconstructed tracks\r
+\r
+  TH1F        *fHistRECeta;     // eta of reconstructed tracks\r
+  TH1F        *fHistMCeta;      // eta of MC primaries\r
+  TH1F        *fHistMCNReta;    // eta of not reconstructed primaries\r
+  TH1F        *fHistFAKEeta;    // eta of fake track\r
+  TH1F        *fHistMULTeta;    // eta of multiply reconstructed tracks\r
+\r
+  TH1F        *fHistRECphi;     // phi of reconstructed tracks\r
+  TH1F        *fHistMCphi;      // phi of MC primaries\r
+  TH1F        *fHistMCNRphi;    // phi of not reconstructed primaries\r
+  TH1F        *fHistFAKEphi;    // phi of fake track\r
+  TH1F        *fHistMULTphi;    // phi of multiply reconstructed tracks\r
+  \r
+  TH1F        *fHistRECHPTeta;  // eta of reconstructed high-pT tracks\r
+  TH1F        *fHistMCHPTeta;   // eta of MC primary high-pT tracks \r
+  TH1F        *fHistMCNRHPTeta; // eta of not reconstructed primaries\r
+  TH1F        *fHistFAKEHPTeta; // eta of fake high-pT tracks\r
+  \r
+  TH1F        *fHistRECHPTphi;  // phi of reconstructed high-pT tracks\r
+  TH1F        *fHistMCHPTphi;   // phi of MC primary high-pT tracks \r
+  TH1F        *fHistMCNRHPTphi; // phi of not reconstructed primaries\r
+  TH1F        *fHistFAKEHPTphi; // phi of fake high-pT tracks\r
+\r
+  TH1F        *fHistRecMult;    // Reconstruction multiplicity\r
+  TH1F        *fHistNCluster;   // Number of clusters of suspicious tracks\r
+\r
+  // Checks on ESD data only\r
+  TH1F        *fh1VtxEff;              // Vtx Reconstruction eff for TPC and SPD\r
+  TH2F        *fh2PhiPadRow[kTotalAC]; // TPC pad row vs phi\r
+  TH2F        *fh2PhiLayer[kTotalAC];  // ITS layer vs phi\r
+  TH2F        *fh2MultSpdChips[2];     // Multiplicity vs fired chips\r
+  TH2F        *fh2VtxTpcSpd;           // Vtx correlation TPC <-> SPD\r
+\r
+  // Correlation histos \r
+  TH2F* fh2VertexCorrelation[kTotalVtx];                    // ESD z-vtx vs MC z-vtx\r
+  TH2F* fh2VertexCorrelationShift[kTotalVtx];               // (MC z-vtx - ESD z-vtx) vs MC z-vtx\r
+  TH1F* fh1VertexShift[kTotalVtx];                          // (MC z-vtx - ESD z-vtx) \r
+  TH1F* fh1VertexShiftNorm[kTotalVtx];                      // (MC z-vtx - ESD z-vtx) / (sigma_ESD-z-vtx) histogrammed\r
+\r
+  TH2F* fh2EtaCorrelation;                       // ESD eta vs MC eta\r
+  TH2F* fh2EtaCorrelationShift;                  // (MC eta - ESD eta) vs MC eta\r
+  TH2F* fh2PhiCorrelation;                       // ESD phi vs MC phi\r
+  TH2F* fh2PhiCorrelationShift;                  // (MC phi - ESD phi) vs MC phi\r
+  TH2F* fh2PtCorrelation;                        // ESD pt vs MC pt\r
+  TH2F* fh2PtCorrelationShift;                   // (MC pt - ESD pt) vs MC pt\r
+\r
+ //Tracked Deltaphi histos\r
+  TH1F* fHistDeltaphiprimaries;                // Tracklet dPhi for primaries\r
+  TH1F* fHistDeltaphisecondaries;              // Tracklet dPhi for secondaries\r
+  TH1F* fHistDeltaphireject;                   // Tracklet dPhi for rejected tracks\r
+\r
+  //ITS-TPC matching\r
+  TH1F* fHistTPCITSdeltax;                     // TPC-ITS matching Delta_x\r
+  TH1F* fHistTPCITSdeltay;                     // TPC-ITS matching Delta_y\r
+  TH1F* fHistTPCITSdeltaz;                     // TPC-ITS matching Delta_z\r
+  TH1F* fHistTPCITSdeltar;                     // TPC-ITS matching Delta_r\r
+  //TRD-TPC matching\r
+  TH1F* fHistTPCTRDdeltax;                     // TPC-TRD matching Delta_x\r
+  TH1F* fHistTPCTRDdeltay;                     // TPC-TRD matching Delta_y\r
+  TH1F* fHistTPCTRDdeltaz;                     // TPC-TRD matching Delta_z\r
+  TH1F* fHistTPCTRDdeltar;                     // TPC-TRD matching Delta_r\r
+  TH1F* fHistTPCTRDdeltaphi;                   // TPC-TRD matching Delta_phi\r
+  // Pulls\r
+  TProfile* fPtPullsPos;                          // Pt pulls\r
+  TProfile* fPtPullsNeg;                          // Pt pulls\r
+  TProfile* fPtShiftPos;                          // Pt shift\r
+  TProfile* fPtShiftNeg;                          // Pt shift\r
+  // Others\r
+  Int_t fTrackType;\r
+  AliESDtrackCuts* fCuts;                         // List of cuts\r
+  // SPD Vertex\r
+  TH2F* fVertexRvsZ;\r
+  TH2F* fVertexRvsZC;\r
+  // Multi Fluctuations\r
+  TProfile* fEtaMultiFluc;                    // multiplicity fluctuations\r
+  TH1F*     fEtaMulti;                        // multiplicity fluctuations\r
+  TH1F*     fEtaMultiH;                       // multiplicity fluctuations\r
+  TProfile* fPhiMultiFluc;                    // multiplicity fluctuations\r
+  TH1F*     fPhiMulti;                        // multiplicity fluctuations\r
+  TH1F*     fPhiMultiH;                       // multiplicity fluctuations\r
+\r
+  //for each nuclei type one histo (d, t, 3He, 4He)\r
+  TH1F        *fHistRECptCharge[8];       // pt of reconstructed tracks\r
+  TH1F        *fHistMCptCharge[8];        // pt of MC primaries\r
+  TH1F        *fHistMCNRptCharge[8];      // pt of not reconstructed primaries\r
+  TH1F        *fHistFAKEptCharge[8];      // pt of fake track\r
+\r
+  TH1F        *fHistRECetaCharge[8];      // eta of reconstructed tracks\r
+  TH1F        *fHistMCetaCharge[8];       // eta of MC primaries\r
+  TH1F        *fHistMCNRetaCharge[8];     // eta of not reconstructed primaries\r
+  TH1F        *fHistFAKEetaCharge[8];     // eta of fake track\r
+\r
+  TH1F        *fHistRECphiCharge[8];      // phi of reconstructed tracks\r
+  TH1F        *fHistMCphiCharge[8];       // phi of MC primaries\r
+  TH1F        *fHistMCNRphiCharge[8];     // phi of not reconstructed primaries\r
+  TH1F        *fHistFAKEphiCharge[8];     // phi of fake track\r
+\r
+\r
+  AliAnalysisTaskEfficiency(const AliAnalysisTaskEfficiency&); // not implemented\r
+  AliAnalysisTaskEfficiency& operator=(const AliAnalysisTaskEfficiency&); // not implemented\r
+  \r
+  ClassDef(AliAnalysisTaskEfficiency, 1); // example of analysis\r
+};\r
+\r
+#endif\r
diff --git a/test/vmctest/scripts/efficiency/CreateAnalysisPlugin.C b/test/vmctest/scripts/efficiency/CreateAnalysisPlugin.C
new file mode 100644 (file)
index 0000000..2164e73
--- /dev/null
@@ -0,0 +1,145 @@
+AliAnalysisGrid* CreateAnalysisPlugin(TString analysisMode="full") 
+{ 
+  AliAnalysisAlien *plugin = new AliAnalysisAlien();
+
+  // Overwrite all generated files, datasets and output results from a previous session
+  plugin->SetOverwriteMode();
+  // Set the run mode (can be "full", "test", "offline", "submit" or "terminate") 
+
+  plugin->SetRunMode(TString (analysisMode.Data()));
+
+  // Set versions of used packages 
+  plugin->SetAPIVersion("V1.1x"); 
+  plugin->SetROOTVersion("v5-28-00f");
+  plugin->SetAliROOTVersion("v4-21-33-AN"); 
+  // Declare input data to be processed. 
+
+  // Method 1: Create automatically XML collections using alien 'find' command. 
+  // Define production directory LFN 
+  //    plugin->SetGridDataDir("/alice/sim/LHC10a6"); 
+  //*   plugin->SetGridDataDir("/alice/data/2010/LHC10b"); 
+  // Set data search pattern 
+  //   plugin->SetDataPattern("*ESDs.root");  // simulated, tags not used
+  //*  plugin->SetDataPattern("*ESDs/pass2/*ESDs.root"); // real data check reco pass and data base directory
+  //*  plugin->SetRunPrefix("000");   // real data
+  //    plugin->SetDataPattern("*tag.root");  // Use ESD tags (same applies for AOD's)
+  // ...then add run numbers to be considered 
+  //    plugin->AddRunNumber(125020);    // simulated
+                                                       
+  //for ESDs
+  plugin->SetGridDataDir("/alice/data/2011/LHC11a"); 
+  plugin->SetDataPattern("*ESDs/pass1/*ESDs.root");
+  plugin->SetRunPrefix("000");
+
+  //ESDs sim
+  // plugin->SetGridDataDir("/alice/sim/LHC10f9b");
+  //plugin->SetDataPattern("*ESDs.root"); 
+
+  //for AODs
+  // plugin->SetGridDataDir("/alice/data/2010/LHC10c"); 
+  // plugin->SetRunPrefix("000");   // real data
+  //plugin->SetDataPattern("*ESDs/pass2_recovery_900GeV/AOD017/*AOD.root"); 
+  
+  //sim AODs
+  // plugin->SetGridDataDir("/alice/sim/LHC10d4a"); 
+  //plugin->SetDataPattern("*AOD012/*AOD.root");
+  
+  //  TString runs ="120824:120823:120822:120821:120820:120758:120750:120741:120671:120617:120616:120505:120504:120503:120244:120079:120076:120073:120072:120069:120067:119862:119859:119856:119853:119849:119846:119845:119844:119842:119841:119163:119161:119159:119086:119085:119084:119079:119077:119067:119061:119047:119041:119037"; // dont forget last two runs
+  
+  //TString runs ="120829:120825";
+  // TString runs="118506:118507:118512:118518:118556:118558:118560:118561";
+  TString runs ="146801";
+
+  TObjArray* array = runs.Tokenize ( ":" );
+  TObjString *str;
+  TString strr,strr2_1,strr2_2;
+  for ( Int_t i = 0;i < array->GetEntriesFast();i++ ) {
+    str = ( TObjString * ) array->At ( i );
+    strr = str->GetString();
+    if ( !strr.IsNull() ) {
+      plugin->AddRunNumber(strr.Atoi());
+    }
+  }  
+        
+  // Method 2: Declare existing data files (raw collections, xml collections, root file) 
+  // If no path mentioned data is supposed to be in the work directory (see SetGridWorkingDir()) 
+  // XML collections added via this method can be combined with the first method if 
+  // the content is compatible (using or not tags) 
+  //   plugin->AddDataFile("tag.xml"); 
+  //   plugin->AddDataFile("/alice/data/2008/LHC08c/000057657/raw/Run57657.Merged.RAW.tag.root"); 
+
+  // Define alien work directory where all files will be copied. Relative to alien $HOME. 
+  plugin->SetGridWorkingDir("146801"); 
+  // Declare alien output directory. Relative to working directory. 
+  plugin->SetGridOutputDir("output"); // In this case will be $HOME/work/output 
+  // Declare the analysis source files names separated by blancs. To be compiled runtime 
+  // using ACLiC on the worker nodes. 
+  plugin->SetAnalysisSource("AliAnalysisTaskEfficiency.cxx"); 
+  //    plugin->SetAdditionalRootLibs("CORRFW PWG2resonances");
+  //    plugin->SetAdditionalRootLibs("PWG2resonances");
+  //    plugin->SetAdditionalRootLibs("PWG2resonances");
+  // 
+  plugin->SetAdditionalLibs("AliAnalysisTaskEfficiency.h AliAnalysisTaskEfficiency.cxx");
+  //    plugin->EnablePackage("PWG2resonances");
+  //    plugin->EnablePackage("");
+  //    plugin->EnablePackage("");
+  // Declare all libraries (other than the default ones for the framework. These will be 
+  // loaded by the generated analysis macro. Add all extra files (task .cxx/.h) here. 
+
+  // No need for output file names. Procedure is automatic. 
+  //   plugin->SetOutputFiles("Pt.ESD.1.root"); 
+  //   plugin->SetDefaultOutputs(); 
+  // No need define the files to be archived. Note that this is handled automatically by the plugin.
+  //   plugin->SetOutputArchive("log_archive.zip:stdout,stderr"); 
+  // Set a name for the generated analysis macro (default MyAnalysis.C) Make this unique !
+  plugin->SetAnalysisMacro("AnalysisTest.C"); 
+  // Optionally set maximum number of input files/subjob (default 100, put 0 to ignore). The optimum for an analysis
+  // is correlated with the run time - count few hours TTL per job, not minutes !
+  plugin->SetSplitMaxInputFileNumber(100); 
+  // Optionally set number of failed jobs that will trigger killing waiting sub-jobs. 
+  plugin->SetMaxInitFailed(5); 
+  // Optionally resubmit threshold. 
+  plugin->SetMasterResubmitThreshold(90); 
+  // Optionally set time to live (default 30000 sec) 
+  plugin->SetTTL(20000); 
+  // Optionally set input format (default xml-single) 
+  plugin->SetInputFormat("xml-single"); 
+  // Optionally modify the name of the generated JDL (default analysis.jdl) 
+  plugin->SetJDLName("TaskRsn.jdl"); 
+  // Optionally modify job price (default 1) 
+  plugin->SetPrice(1);  
+  // Optionally modify split mode (default 'se')     
+  plugin->SetSplitMode("se"); 
+   
+   
+  //++++++++++++++++ PROOF ++++++++++++++++
+  //    Proof cluster
+   
+  plugin->SetProofCluster("alice-caf");
+  //plugin->SetProofCluster("skaf.saske.sk");
+  //   plugin->SetProofCluster("skaf-test.saske.sk");
+  // Dataset to be used
+  //     plugin->SetProofDataSet("/alice/sim/LHC10a12_104316#esdTree");
+  //       plugin->SetProofDataSet("/alice/sim/LHC10a12_104157#esdTree");
+  //     plugin->SetProofDataSet("ds.txt");
+  plugin->SetProofDataSet("g4g.txt");
+  // May need to reset proof. Supported modes: 0-no reset, 1-soft, 2-hard
+  plugin->SetProofReset(0);
+  // May limit the number of workers per slave. If used with SetNproofWorkers, SetParallel(nproofworkers) will be called after connection
+  //plugin->SetNproofWorkersPerSlave(1);
+  // May request connection to alien upon connection to grid
+  //     plugin->SetProofConnectGrid(kTRUE);
+    
+  // plugin->SetNproofWorkers(51);
+  // May use a specific version of root installed in proof
+  //     plugin->SetRootVersionForProof("current");
+  // May set the aliroot mode. Check http://aaf.cern.ch/node/83
+  plugin->SetAliRootMode("default"); // Loads AF libs by default
+  // May request ClearPackages (individual ClearPackage not supported)
+  plugin->SetClearPackages(kFALSE);
+  // Plugin test mode works only providing a file containing test file locations
+  //    plugin->SetFileForTestMode("AOD.txt");
+  plugin->SetFileForTestMode("test.txt");
+  //++++++++++++++ end PROOF ++++++++++++++++
+  return plugin; 
+} 
diff --git a/test/vmctest/scripts/efficiency/SetupAnalysis.C b/test/vmctest/scripts/efficiency/SetupAnalysis.C
new file mode 100644 (file)
index 0000000..94648ff
--- /dev/null
@@ -0,0 +1,148 @@
+void SetupAnalysis(TString mode,
+                  TString analysisMode="full",
+                  Bool_t useMC = kFALSE,
+                  Int_t nEvents=1.0*1e9, 
+                  Int_t nEventsSkip=0,
+                  TString format="esd")
+{
+  
+  // ALICE stuff
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) mgr = new AliAnalysisManager("CAF train");
+  
+  // Create and configure the alien handler plugin 
+  gROOT->LoadMacro("CreateAnalysisPlugin.C"); 
+  AliAnalysisGrid *alienHandler = CreateAnalysisPlugin(analysisMode);   
+  if (!alienHandler) return;
+  mgr->SetGridHandler(alienHandler);
+  
+  // input handler for esd or AOD, real or MC data
+  InputHandlerSetup(format,useMC);
+
+  // physics selection
+  if(!format.CompareTo("esd")){
+    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+    AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(kFALSE);
+    if(useMC) physSelTask->GetPhysicsSelection()->SetAnalyzeMC();   
+    AliPhysicsSelection* physSel = physSelTask->GetPhysicsSelection();
+    physSel->AddBackgroundIdentification(new AliBackgroundSelection());
+  }
+  
+  gROOT->ProcessLine(Form(".include %s/include",gSystem->ExpandPathName("$ALICE_ROOT")));
+  
+  gROOT->LoadMacro("AliAnalysisTaskEfficiency.cxx+g");
+  
+  // load and run AddTask macro
+  gROOT->LoadMacro("AddTaskEfficiency.C");
+  
+
+  AliAnalysisTaskSE* task1 = AddTaskEfficiency(-1);
+  if(!task1){
+    Printf("AddTask could not be run.");
+  }
+
+  // Run analysis
+  mgr->InitAnalysis();
+  
+  if ((!mode.CompareTo("proof")) ||(!mode.CompareTo("local"))) {
+    mgr->StartAnalysis(mode.Data(),nEvents,nEventsSkip);
+  }
+  else {
+    mgr->StartAnalysis(mode.Data());
+    
+  }
+  
+}
+
+
+TString GetFormatFromDataSet(TString dataset) {
+  
+  TString dsTreeName;
+  if (dataset.Contains("#")) {
+    Info("runSKAF.C",Form("Detecting format from dataset name '%s' ...",dataset.Data()));
+    dsTreeName=dataset(dataset.Last('#'),dataset.Length());
+  } else {
+    Info("runSKAF.C",Form("Detecting format from dataset '%s' (may take while, depends on network connection) ...",
+                         dataset.Data()));
+    TFileCollection *ds = gProof->GetDataSet(dataset.Data());
+    if (!ds) {
+      Error(Form("Dataset %s doesn't exist on proof cluster!!!!",dataset.Data()));
+      return "";
+    }
+    dsTreeName = ds->GetDefaultTreeName();
+  }
+  
+  if (dsTreeName.Contains("esdTree")) {
+    Info("runSKAF.C","ESD input format detected ...");
+    return "ESD";
+  } else if (dsTreeName.Contains("aodTree"))  {
+    Info("runSKAF.C","AOD input format detected ...");
+    return "AOD";
+  } else {
+    Error("runSKAF.C",Form("Tree %s is not supported !!!",dsTreeName.Data()));
+    Error("runSKAF.C",Form("Maybe set your DS to %s#esdTree or %s#aodTree",
+                          dataset.Data(),dataset.Data()));
+  }
+  
+  return "";
+}
+
+Bool_t InputHandlerSetup(TString format = "esd", Bool_t useKine = kTRUE)
+{
+  format.ToLower();
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+
+  AliAnalysisDataContainer *cin = mgr->GetCommonInputContainer();
+
+  if (cin) return;
+
+  if (!format.CompareTo("esd"))
+    {
+      AliESDInputHandler *esdInputHandler = 
+       dynamic_cast<AliESDInputHandler*>(AliAnalysisManager::GetAnalysisManager()->
+                                         GetInputEventHandler());
+
+      if (!esdInputHandler)
+       {
+         Info("CustomAnalysisTaskInputSetup", "Creating esdInputHandler ...");
+         esdInputHandler = new AliESDInputHandler();
+         mgr->SetInputEventHandler(esdInputHandler);
+       }
+
+      if (useKine)
+       {
+         AliMCEventHandler* mcInputHandler = 
+           dynamic_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->
+                                            GetMCtruthEventHandler());
+
+         if (!mcInputHandler)
+           {
+             Info("CustomAnalysisTaskInputSetup", "Creating mcInputHandler ...");
+             AliMCEventHandler* mcInputHandler = new AliMCEventHandler();
+             mgr->SetMCtruthEventHandler(mcInputHandler);
+           }
+       }
+
+    }
+  else if (!format.CompareTo("aod"))
+    {
+      AliAODInputHandler *aodInputHandler = 
+       dynamic_cast<AliAODInputHandler*>(AliAnalysisManager::GetAnalysisManager()->
+                                         GetInputEventHandler());
+
+      if (!aodInputHandler)
+       {
+         Info("CustomAnalysisTaskInputSetup", "Creating aodInputHandler ...");
+         aodInputHandler = new AliAODInputHandler();
+         mgr->SetInputEventHandler(aodInputHandler);
+       }
+    }
+  else
+    {
+      AliWarning("Wrong input format!!! Only ESD and AOD are supported. Skipping Task ...");
+      return kFALSE;
+    }
+
+  return kTRUE;
+}
diff --git a/test/vmctest/scripts/efficiency/doubleRatio.C b/test/vmctest/scripts/efficiency/doubleRatio.C
new file mode 100644 (file)
index 0000000..a0b15e4
--- /dev/null
@@ -0,0 +1,310 @@
+// plot macro for double ratio of reconstruction efficiencies of different particles
+// Author: Eva Sicking
+
+void doubleRatio(Int_t trackType=0, Int_t particle=3)
+{
+
+  //particle names
+  TString partName[4]={"Deu", "Tri", "He3", "He4"};
+
+  // track type names
+  TString trackName[4]={"Global", "TPC", "ITS_SA", "ITS"};
+
+  //open file and lists
+  //output of AliAnalysisTaskEfficiency
+  TFile *f = new TFile("AnalysisResults.root");
+  TList *list = (TList *)f->Get(Form("QAHists/QAHists_%s",trackName[trackType].Data()));
+
+  //style
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptStat(0);
+  const Int_t NRGBs = 5;
+  const Int_t NCont = 500;
+  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+  Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
+  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
+  Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
+  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
+  gStyle->SetNumberContours(NCont);
+    
+  // Pt     
+  TH1F* fHistRECpt  =  (TH1F*) list->
+    FindObject(Form("fHistRECptChargeAnti%s",partName[particle].Data()));
+  TH1F* fHistMCpt   =  (TH1F*) list->
+    FindObject(Form("fHistMCptChargeAnti%s",partName[particle].Data()));
+  TH1F* fHistFAKEpt =  (TH1F*) list->
+    FindObject(Form("fHistFAKEptChargeAnti%s",partName[particle].Data()));
+  TH1F* fHistMCNRpt =  (TH1F*) list->
+    FindObject(Form("fHistMCNRptChargeAnti%s",partName[particle].Data()));
+  
+    
+  c1 = new TCanvas(Form("pT_antinuclei_%s",trackName[trackType].Data()),
+                  Form("pT_antinuclei_%s",trackName[trackType].Data()),
+                  100, 100, 1000,600 );
+  c1->Divide(2,2);
+  c1->cd(1);
+  c1->GetPad(1)->SetLogy();
+    
+  fHistMCpt->SetXTitle("p_{T} [GeV]");
+  fHistMCpt->SetMinimum(1);
+  fHistMCpt->Draw();
+  fHistMCpt->SetTitle("");
+  fHistRECpt->SetLineColor(2);
+  fHistRECpt->Draw("same");
+  fHistFAKEpt->SetLineColor(4);
+  fHistFAKEpt->Draw("same");
+  fHistMCNRpt->SetLineColor(3);
+  fHistMCNRpt->Draw("same");
+
+  c1->cd(4);
+  fHistPtInEff = (TH1F*) fHistMCNRpt->Clone();
+  fHistPtInEff->Divide(fHistMCpt);
+  fHistPtInEff->SetXTitle("p_{T} [GeV]");
+  fHistPtInEff->SetTitle("Inefficiency from non-rec. MC tracks");
+  fHistPtInEff->Draw();
+  fHistPtInEff->SetMinimum(0);
+  fHistPtInEff->SetMaximum(1);
+
+  c1->cd(3);
+  fHistPtEff = (TH1F*) fHistRECpt->Clone();
+  fHistPtEff->Add(fHistFAKEpt, -1.);
+  fHistPtEff->Divide(fHistMCpt);
+  //fHistPtEff->Add(fHistPtInEff);
+  fHistPtEff->SetXTitle("p_{T} [GeV]");
+  fHistPtEff->SetTitle("Efficiency");
+  fHistPtEff->Draw();
+  fHistPtEff->SetMinimum(0);
+  fHistPtEff->SetMaximum(1);
+
+  c1->cd(2);
+  TLine *line = new TLine(0.07, 0.75, 0.2, 0.75);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.73, "MC");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+  line = new TLine(0.07, 0.55, 0.2, 0.55);
+  line->SetLineColor(2);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.53, "Reconstructed");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+  line = new TLine(0.07, 0.35, 0.2, 0.35);
+  line->SetLineColor(4);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.33, "Fake Tracks");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+
+  line = new TLine(0.07, 0.15, 0.2, 0.15);
+  line->SetLineColor(3);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.13, "MC not reconstructed");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+
+
+  // Pt     
+  TH1F* fHistRECpt2  =  (TH1F*) list->
+    FindObject(Form("fHistRECptCharge%s",partName[particle].Data()));
+  TH1F* fHistMCpt2   =  (TH1F*) list->
+    FindObject(Form("fHistMCptCharge%s",partName[particle].Data()));
+  TH1F* fHistFAKEpt2 =  (TH1F*) list->
+    FindObject(Form("fHistFAKEptCharge%s",partName[particle].Data()));
+  TH1F* fHistMCNRpt2 =  (TH1F*) list->
+    FindObject(Form("fHistMCNRptCharge%s",partName[particle].Data()));
+    
+  c2 = new TCanvas(Form("pT_nuclei_%s",trackName[trackType].Data()),
+                  Form("pT_nuclei_%s",trackName[trackType].Data()),
+                  100, 100, 1000,600 );
+  c2->Divide(2,2);
+  c2->cd(1);
+  c2->GetPad(1)->SetLogy();
+    
+  fHistMCpt2->SetXTitle("p_{T} [GeV]");
+  fHistMCpt2->SetMinimum(1);
+  fHistMCpt2->Draw();
+  fHistMCpt2->SetTitle("");
+  fHistRECpt2->SetLineColor(2);
+  fHistRECpt2->Draw("same");
+  fHistFAKEpt2->SetLineColor(4);
+  fHistFAKEpt2->Draw("same");
+  fHistMCNRpt2->SetLineColor(3);
+  fHistMCNRpt2->Draw("same");
+
+  c2->cd(4);
+  fHistPtInEff2 = (TH1F*) fHistMCNRpt2->Clone();
+  fHistPtInEff2->Divide(fHistMCpt2);
+  fHistPtInEff2->SetXTitle("p_{T} [GeV]");
+  fHistPtInEff2->SetTitle("Inefficiency from non-rec. MC tracks");
+  fHistPtInEff2->Draw();
+  fHistPtInEff2->SetMinimum(0);
+  fHistPtInEff2->SetMaximum(1);
+
+  c2->cd(3);
+  fHistPtEff2 = (TH1F*) fHistRECpt2->Clone();
+  fHistPtEff2->Add(fHistFAKEpt2, -1.);
+  fHistPtEff2->Divide(fHistMCpt2);
+  //    fHistPtEff->Add(fHistPtInEff);
+  fHistPtEff2->SetXTitle("p_{T} [GeV]");
+  fHistPtEff2->SetTitle("Efficiency");
+  fHistPtEff2->Draw();
+  fHistPtEff2->SetMinimum(0);
+  fHistPtEff2->SetMaximum(1);
+
+  c2->cd(2);
+  TLine *line = new TLine(0.07, 0.75, 0.2, 0.75);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.73, "MC");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+  line = new TLine(0.07, 0.55, 0.2, 0.55);
+  line->SetLineColor(2);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.53, "Reconstructed");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+  line = new TLine(0.07, 0.35, 0.2, 0.35);
+  line->SetLineColor(4);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.33, "Fake Tracks");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+
+  line = new TLine(0.07, 0.15, 0.2, 0.15);
+  line->SetLineColor(3);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.13, "MC not reconstructed");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+
+
+
+  new TCanvas;
+  
+  fHistPtEffPosNeg = (TH1F*) fHistPtEff->Clone();
+  fHistPtEffPosNeg->Divide(fHistPtEff2);
+  fHistPtEffPosNeg->SetLineWidth(2);
+  fHistPtEffPosNeg->Draw();
+  fHistPtEffPosNeg->SetMinimum(0);
+  fHistPtEffPosNeg->SetMaximum(1.2);
+  fHistPtEff->Draw("same");//anti-particle
+  fHistPtEff->SetLineColor(kRed);
+  fHistPtEff->SetLineWidth(2);
+  fHistPtEff2->Draw("same");
+  fHistPtEff2->SetLineColor(kBlack);
+  fHistPtEff2->SetLineWidth(2);
+  TF1 *func = new TF1("func", "1",-5,20);
+  func->SetLineWidth(1);
+  func->Draw("same");
+
+
+
+
+
+
+  TLegend *legp= new TLegend(0.65,0.78,0.9,0.98);
+  legp->SetFillColor(kWhite);
+  legp->SetBorderSize(0);
+  Int_t number =4;
+  TF1 *fun[2];
+  for(Int_t i=0;i<2;i++){
+    fun[i]= new TF1(Form("fun%d",i),"gaus",-5.0,5.0);
+    fun[i]->SetLineColor(2**i);
+    fun[i]->SetLineStyle(1);
+  }
+  legp->AddEntry(fun[0],Form("%s",partName[particle].Data() ),"l");   
+  legp->AddEntry(fun[1],Form("Anti-%s",partName[particle].Data() ),"l");   
+
+
+
+  
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptStat(0);
+  gStyle->SetPalette(1);
+
+
+  TCanvas *c = new TCanvas("c","PtAll" , 100, 100, 880, 680);
+  c->Draw();
+  c ->cd();
+  TPad *pad_spectrum = new TPad( "pad_spectrum","pad_spectrum", 0.0,0.39,1.0,0.95);
+  c ->cd();
+  TPad *pad_ratio = new TPad( "pad_ratio","pad_ratio", 0.0,0.03,1.0,0.39);  
+  
+  
+  c ->cd();
+  pad_spectrum->SetTopMargin(0.);
+  pad_spectrum->SetBottomMargin(0.);
+  pad_spectrum->SetLeftMargin(0.12);
+  pad_spectrum->SetRightMargin(0.055);
+  pad_spectrum->Draw();
+  pad_spectrum->cd(); 
+  TH1F * histo = new TH1F("histo", "",150, -0.5,149.5 );
+  histo->GetXaxis()->SetRangeUser(0,11);
+  //  histo->GetXaxis()->SetRangeUser(25,35);
+  histo->GetYaxis()->SetRangeUser(0.01,1.3);
+  histo->SetXTitle("N_{charged}");
+  histo->SetTitle("");
+  histo->SetYTitle("efficiency  ");
+  histo->Draw();
+  histo->GetYaxis()->SetNdivisions(10);
+  histo->GetYaxis()->SetLabelSize(0.07);
+  histo->GetYaxis()->SetTitleSize(0.08);
+  histo->GetYaxis()->SetTitleOffset(0.6);
+  fHistPtEff->Draw("same");
+  fHistPtEff2->Draw("same");
+
+  legp->Draw();
+  
+  
+
+  c ->cd();
+  pad_ratio->SetTopMargin(0.);
+  pad_ratio->SetBottomMargin(0.36);
+  pad_ratio->SetLeftMargin(0.12);
+  pad_ratio->SetRightMargin(0.055);
+  pad_ratio->Draw();
+  pad_ratio->cd(); 
+  TH1F * histo2 = new TH1F("histo2", "",150, -0.5,149.5 );
+  histo2->GetXaxis()->SetRangeUser(0,90);
+  //  histo->GetXaxis()->SetRangeUser(25,35);
+  histo2->SetXTitle("p_{T} (GeV/c)");
+  histo2->SetTitle("");
+  histo2->SetYTitle("Ratio Neg / Pos");
+  histo2->Draw();
+  histo2->GetXaxis()->SetRangeUser(0.0,11);
+  histo2->SetLabelSize(0.1,"X");
+  histo2->SetLabelSize(0.1,"Y");
+  histo2->SetTitleSize(0.09,"X");
+  histo2->SetTitleSize(0.09,"Y");
+  histo2->SetTitleOffset(0.3,"Y");
+  histo2->GetYaxis()->SetRangeUser(0.59,1.16);
+  histo2->GetYaxis()->SetNdivisions(4);
+  histo2->GetYaxis()->SetLabelSize(0.12);
+  histo2->GetXaxis()->SetLabelSize(0.12);
+  histo2->GetXaxis()->SetTitleSize(0.15);
+  histo2->GetYaxis()->SetTitleSize(0.11);
+
+  fHistPtEffPosNeg->Draw("same");
+  func->Draw("same");
+  func->SetLineStyle(2);
+
+  c->Print(Form("%s.png",partName[particle].Data()));
+
+}
diff --git a/test/vmctest/scripts/efficiency/macrosalida.C b/test/vmctest/scripts/efficiency/macrosalida.C
new file mode 100644 (file)
index 0000000..2e3bf1a
--- /dev/null
@@ -0,0 +1,290 @@
+// monitor reconstruction efficiencies for different track types
+// use output of AliAnalysisTaskEfficiency
+// Authors: Veronica Canoa Roman, Eva Sicking
+
+void macrosalida(Int_t prod=0,Int_t trackType=0 )
+{
+
+  //postfix of analysis output files ("X")
+  TString prodName[6]={"A","B","C","D","E","F"};
+  //names of physics lists using in analysed productions (LHC11d6x)
+  TString prodNameDetail[6]={"Geant4_QGSP_BERT_EMV_p02",
+                            "Geant4_QGSP_BERT_CHIPS_p02",
+                            "Geant4_QGSP_BERT_EMV_b01",
+                            "Geant4_QGSP_BERT_CHIPS_b01",
+                            "Geant4_QGSP_FTFP_BERT_b01",
+                            "Geant3"};
+  // name of track types
+  // ITS = left over tracks from global tracking, efficiency does not make sense
+  TString trackName[4]={"Global", "TPC", "ITS_SA", "ITS"};
+  
+  //open file
+  // TFile *f = new TFile(Form("AnalysisResults%s.root",prodName[prod].Data()));
+  TFile *f = new TFile(Form("AnalysisResults.root",prodName[prod].Data()));
+  TList *list = (TList *)f->Get(Form("QAHists/QAHists_%s",trackName[trackType].Data()));
+
+  //style
+  gROOT->SetStyle("Plain");
+  gStyle->SetOptStat(0);
+  const Int_t NRGBs = 5;
+  const Int_t NCont = 500;
+  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
+  Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
+  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
+  Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
+  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
+  gStyle->SetNumberContours(NCont);
+    
+
+  
+  // Pt efficiency histos     
+  TH1F* fHistRECpt  =  (TH1F*) list->FindObject("fHistRECpt");
+  TH1F* fHistMCpt   =  (TH1F*) list->FindObject("fHistMCpt");
+  TH1F* fHistFAKEpt =  (TH1F*) list->FindObject("fHistFAKEpt");
+  TH1F* fHistMCNRpt =  (TH1F*) list->FindObject("fHistMCNRpt");
+  
+    
+  c1 = new TCanvas(Form("pT_%s_%s",prodNameDetail[prod].Data(),
+                       trackName[trackType].Data()),
+                  Form("pT_%s_%s",prodNameDetail[prod].Data(),
+                       trackName[trackType].Data()),
+                  100, 100, 1000,800 );
+  c1->Divide(2,2);
+  c1->cd(1);
+  c1->GetPad(1)->SetLogy();
+    
+  fHistMCpt->SetXTitle("p_{T} [GeV]");
+  fHistMCpt->SetMinimum(1);
+  fHistMCpt->Draw();
+  fHistMCpt->SetTitle("");
+  fHistMCpt->SetLineWidth(2);
+  fHistRECpt->SetLineColor(2);
+  fHistRECpt->SetLineWidth(2);
+  fHistRECpt->Draw("same");
+  fHistFAKEpt->SetLineColor(4);
+  fHistFAKEpt->SetLineWidth(2);
+  fHistFAKEpt->Draw("same");
+  fHistMCNRpt->SetLineColor(3);
+  fHistMCNRpt->SetLineWidth(2);
+  fHistMCNRpt->Draw("same");
+
+  c1->cd(4);
+  fHistPtInEff = (TH1F*) fHistMCNRpt->Clone();
+  fHistPtInEff->Divide(fHistMCpt);
+  fHistPtInEff->SetXTitle("p_{T} [GeV]");
+  fHistPtInEff->SetTitle("Inefficiency from non-rec. MC tracks");
+  fHistPtInEff->Draw();
+  fHistPtInEff->SetMinimum(0);
+  fHistPtInEff->SetMaximum(1);
+
+  c1->cd(3);
+  fHistPtEff = (TH1F*) fHistRECpt->Clone();
+  fHistPtEff->Add(fHistFAKEpt, -1.);
+  fHistPtEff->Divide(fHistMCpt);
+  //    fHistPtEff->Add(fHistPtInEff);
+  fHistPtEff->SetXTitle("p_{T} [GeV]");
+  fHistPtEff->SetTitle("Efficiency");
+  fHistPtEff->Draw();
+  fHistPtEff->SetMinimum(0);
+  fHistPtEff->SetMaximum(1);
+
+  c1->cd(2);
+  TLine *line = new TLine(0.07, 0.75, 0.2, 0.75);
+  line->SetLineWidth(3);
+  line->Draw();
+  tex = new TLatex(0.25, 0.73, "MC");
+  tex->SetLineWidth(2);
+  tex->SetTextSize(0.08);
+  tex->Draw();
+
+  line = new TLine(0.07, 0.55, 0.2, 0.55);
+  line->SetLineColor(2);
+  line->SetLineWidth(3);
+  line->Draw();
+  tex = new TLatex(0.25, 0.53, "Reconstructed");
+  tex->SetLineWidth(2);
+  tex->SetTextSize(0.08);
+  tex->Draw();
+
+  line = new TLine(0.07, 0.35, 0.2, 0.35);
+  line->SetLineColor(4);
+  line->SetLineWidth(3);
+  line->Draw();
+  tex = new TLatex(0.25, 0.33, "Fake Tracks");
+  tex->SetLineWidth(3);
+  tex->SetTextSize(0.08);
+  tex->Draw();
+
+
+  line = new TLine(0.07, 0.15, 0.2, 0.15);
+  line->SetLineColor(3);
+  line->SetLineWidth(3);
+  line->Draw();
+  tex = new TLatex(0.25, 0.13, "MC not reconstructed");
+  tex->SetLineWidth(2);
+  tex->SetTextSize(0.08);
+  tex->Draw();
+
+  // c1->Print(Form("pT_%s_%s.png",prodNameDetail[prod].Data(),trackName[trackType].Data()));
+
+  // Eta efficiency histos     
+  c2 = new TCanvas(Form("eta_%s_%s",prodNameDetail[prod].Data(),
+                       trackName[trackType].Data()),
+                  Form("eta_%s_%s",prodNameDetail[prod].Data(),
+                       trackName[trackType].Data()),
+                  100, 100, 1000,600 );  
+  c2->Divide(2,2);
+  c2->cd(1);
+
+  TH1F* fHistRECeta  =  (TH1F*) list->FindObject("fHistRECeta");
+  TH1F* fHistMCeta   =  (TH1F*) list->FindObject("fHistMCeta");
+  TH1F* fHistFAKEeta =  (TH1F*) list->FindObject("fHistFAKEeta");
+  TH1F* fHistMCNReta =  (TH1F*) list->FindObject("fHistMCNReta");
+  fHistMCeta->SetMinimum(0.);
+  fHistMCeta->SetXTitle("#eta");
+  fHistMCeta->Draw();
+  fHistMCeta->SetTitle("");
+  fHistRECeta->SetLineColor(2);
+  fHistRECeta->Draw("same");
+  fHistFAKEeta->SetLineColor(4);
+  fHistFAKEeta->Draw("same");
+  fHistMCNReta->SetLineColor(3);
+  fHistMCNReta->Draw("same");
+  c2->cd(4);
+  fHistEtaInEff = (TH1F*) fHistMCNReta->Clone();
+  fHistEtaInEff->Divide(fHistMCeta);
+  fHistEtaInEff->SetXTitle("#eta");
+  fHistEtaInEff->SetTitle("Inefficiency from non-rec. MC tracks");
+  fHistEtaInEff->Draw();
+  fHistEtaInEff->SetMinimum(0);
+  fHistEtaInEff->SetMaximum(1);
+
+  c2->cd(3);
+  fHistEtaEff = (TH1F*) fHistRECeta->Clone();
+  fHistEtaEff->Add(fHistFAKEeta, -1.);
+  fHistEtaEff->Divide(fHistMCeta);
+  fHistEtaEff->SetXTitle("#eta");
+  fHistEtaEff->SetTitle("Efficiency");
+  fHistEtaEff->Draw();
+  fHistEtaEff->SetMinimum(0);
+  fHistEtaEff->SetMaximum(1);
+
+  c2->cd(2);
+  TLine *line = new TLine(0.07, 0.75, 0.2, 0.75);
+  line->Draw();
+  line->SetLineWidth(2);
+  tex = new TLatex(0.25, 0.73, "MC");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+  line = new TLine(0.07, 0.55, 0.2, 0.55);
+  line->SetLineColor(2);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.53, "Reconstructed");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+  line = new TLine(0.07, 0.35, 0.2, 0.35);
+  line->SetLineColor(4);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.33, "Fake Tracks");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+
+  line = new TLine(0.07, 0.15, 0.2, 0.15);
+  line->SetLineColor(3);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.13, "MC not reconstructed");
+  tex->SetLineWidth(2);
+  tex->Draw();
+  
+  //c2->Print(Form("eta_%s_%s.png",prodNameDetail[prod].Data(),trackName[trackType].Data()));
+
+
+  // Phi efficiency histos     
+  c3 = new TCanvas(Form("phi_%s_%s",prodNameDetail[prod].Data(),
+                       trackName[trackType].Data()),
+                  Form("phi_%s_%s",prodNameDetail[prod].Data(),
+                       trackName[trackType].Data()),
+                  100, 100, 1000,600 );  
+  c3->Divide(2,2);
+  c3->cd(1);
+
+  TH1F* fHistRECphi  =  (TH1F*) list->FindObject("fHistRECphi");
+  TH1F* fHistMCphi   =  (TH1F*) list->FindObject("fHistMCphi");
+  TH1F* fHistFAKEphi =  (TH1F*) list->FindObject("fHistFAKEphi");
+  TH1F* fHistMCNRphi =  (TH1F*) list->FindObject("fHistMCNRphi");
+    
+  fHistMCphi->SetMinimum(0.);
+  fHistMCphi->SetXTitle("#phi");
+  fHistMCphi->Draw();
+  fHistMCphi->SetTitle("");
+  fHistRECphi->SetLineColor(2);
+  fHistRECphi->Draw("same");
+  fHistFAKEphi->SetLineColor(4);
+  fHistFAKEphi->Draw("same");
+  fHistMCNRphi->SetLineColor(3);
+  fHistMCNRphi->Draw("same");
+
+  c3->cd(4);
+    
+  fHistPhiInEff = (TH1F*) fHistMCNRphi->Clone();
+  fHistPhiInEff->Divide(fHistMCphi);
+  fHistPhiInEff->SetXTitle("#phi");
+  fHistPhiInEff->SetTitle("Inefficiency from non-rec. MC tracks");
+  fHistPhiInEff->Draw();
+  fHistPhiInEff->SetMinimum(0);
+  fHistPhiInEff->SetMaximum(1);
+
+  c3->cd(3);
+  fHistPhiEff = (TH1F*) fHistRECphi->Clone();
+  fHistPhiEff->Add(fHistFAKEphi, -1.);
+  fHistPhiEff->Divide(fHistMCphi);
+  fHistPhiEff->SetXTitle("#phi");
+  fHistPhiEff->SetTitle("Efficiency");
+  fHistPhiEff->Draw();
+  fHistPhiEff->SetMinimum(0);
+  fHistPhiEff->SetMaximum(1);
+
+
+  c3->cd(2);
+  TLine *line = new TLine(0.07, 0.75, 0.2, 0.75);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.73, "MC");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+  line = new TLine(0.07, 0.55, 0.2, 0.55);
+  line->SetLineColor(2);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.53, "Reconstructed");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+  line = new TLine(0.07, 0.35, 0.2, 0.35);
+  line->SetLineColor(4);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.33, "Fake Tracks");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+
+  line = new TLine(0.07, 0.15, 0.2, 0.15);
+  line->SetLineColor(3);
+  line->SetLineWidth(2);
+  line->Draw();
+  tex = new TLatex(0.25, 0.13, "MC not reconstructed");
+  tex->SetLineWidth(2);
+  tex->Draw();
+
+  // c3->Print(Form("phi_%s_%s.png",prodNameDetail[prod].Data(),trackName[trackType].Data()));
+
+
+}
diff --git a/test/vmctest/scripts/efficiency/runALICE.C b/test/vmctest/scripts/efficiency/runALICE.C
new file mode 100644 (file)
index 0000000..e918087
--- /dev/null
@@ -0,0 +1,65 @@
+void runALICE(TString mode="local",
+             TString analysisMode="full",
+             Bool_t useMC = kTRUE,
+             Int_t nEvents=1.0*1e6, 
+             Int_t nEventsSkip=0,
+             TString format="esd") {
+  
+  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("libCORRFW.so");
+
+  
+  // Create and configure the alien handler plugin 
+  gROOT->LoadMacro("SetupAnalysis.C"); 
+  SetupAnalysis(mode.Data(),
+               analysisMode.Data(),
+               useMC,
+               nEvents,
+               nEventsSkip,
+               format.Data());
+  
+}
+
+Int_t setupPar(const char* pararchivename) {
+  ///////////////////
+  // Setup PAR File//
+  ///////////////////
+  if (pararchivename) {
+    char processline[1024];
+    sprintf(processline,".! tar xzf %s.par",pararchivename);
+    gROOT->ProcessLine(processline);
+    const char* ocwd = gSystem->WorkingDirectory();
+    gSystem->ChangeDirectory(pararchivename);
+                                                                                                                                               
+    // check for BUILD.sh and execute
+    if (!gSystem->AccessPathName("PROOF-INF/BUILD.sh")) {
+      //       printf("*******************************\n");
+      //       printf("*** Building PAR archive    ***\n");
+      //       printf("*******************************\n");
+      
+      if (gSystem->Exec("PROOF-INF/BUILD.sh")) {
+        Error("runAnalysis","Cannot Build the PAR Archive! - Abort!");
+        return -1;
+      }
+    }
+    // check for SETUP.C and execute
+    if (!gSystem->AccessPathName("PROOF-INF/SETUP.C")) {
+      //       printf("*******************************\n");
+      //       printf("*** Setup PAR archive       ***\n");
+      //       printf("*******************************\n");
+      gROOT->Macro("PROOF-INF/SETUP.C");
+    }
+    
+    gSystem->ChangeDirectory("../");
+  }
+  return 1;
+}