Add Kinematics reader
authorakisiel <akisiel@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 14 Oct 2011 15:54:40 +0000 (15:54 +0000)
committerakisiel <akisiel@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 14 Oct 2011 15:54:40 +0000 (15:54 +0000)
16 files changed:
PWG2/CMakelibPWG2femtoscopy.pkg
PWG2/CMakelibPWG2femtoscopyUser.pkg
PWG2/FEMTOSCOPY/AliFemto/AliAnalysisTaskFemto.cxx
PWG2/FEMTOSCOPY/AliFemto/AliAnalysisTaskFemto.h
PWG2/FEMTOSCOPY/AliFemto/AliFemtoCutMonitorEventMult.cxx
PWG2/FEMTOSCOPY/AliFemto/AliFemtoCutMonitorEventMult.h
PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.cxx
PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderESDChainKine.h
PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderKinematicsChain.cxx [new file with mode: 0644]
PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderKinematicsChain.h [new file with mode: 0644]
PWG2/FEMTOSCOPY/AliFemto/AliFemtoMCTrackCut.cxx [new file with mode: 0644]
PWG2/FEMTOSCOPY/AliFemto/AliFemtoMCTrackCut.h [new file with mode: 0644]
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoPairCutAntiGamma.cxx
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoPairCutAntiGamma.h
PWG2/PWG2femtoscopyLinkDef.h
PWG2/PWG2femtoscopyUserLinkDef.h

index 54864aa..a3ca1b4 100644 (file)
@@ -83,7 +83,9 @@ set ( SRCS
     FEMTOSCOPY/AliFemto/AliFemtoCutMonitor.cxx 
     FEMTOSCOPY/AliFemto/AliFemtoCorrFctn.cxx 
     FEMTOSCOPY/AliFemto/AliAnalysisTaskFemto.cxx 
-    FEMTOSCOPY/AliFemto/AliTwoTrackRes.cxx 
+    FEMTOSCOPY/AliFemto/AliTwoTrackRes.cxx
+    FEMTOSCOPY/AliFemto/AliFemtoMCTrackCut.cxx
+    FEMTOSCOPY/AliFemto/AliFemtoEventReaderKinematicsChain.cxx
     )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
@@ -129,6 +131,8 @@ set ( HDRS ${HDRS}
     FEMTOSCOPY/AliFemto/AliFemtoPicoEventCollectionVector.h 
     FEMTOSCOPY/AliFemto/AliFemtoEventWriter.h 
     FEMTOSCOPY/AliFemto/AliFemtoEventWriterCollection.h 
+    FEMTOSCOPY/AliFemto/AliFemtoMCTrackCut.h
+    FEMTOSCOPY/AliFemto/AliFemtoEventReaderKinematicsChain.h
     )
 
 set ( FSRCS  FEMTOSCOPY/AliFemto/AliFemtoFsiTools.F FEMTOSCOPY/AliFemto/AliFemtoFsiWeightLednicky.F)
index 426a042..39b5d62 100644 (file)
@@ -25,7 +25,7 @@
 # SHLIBS - Shared Libraries and objects for linking (Executables only)           #
 #--------------------------------------------------------------------------------#
 
-set ( SRCS  FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityKTPairCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityTPCEntranceSepPairCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoTPCInnerCorrFctn.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoChi2CorrFctn.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnSource.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnDEtaDPhi.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelGausRinvFreezeOutGenerator.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnTrueQ.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnNonIdDR.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctn3DSpherical.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctn3DLCMSSpherical.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCutMonitorParticleMomRes.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelAllHiddenInfo.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCutMonitorTrackTPCchiNdof.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoQATrackCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoQAEventCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCutMonitorTrackTPCncls.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityTPCEntranceSepQAPairCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityQAPairCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityCorrFctn.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoYlm.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnDirectYlm.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnDirectYlm.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCutMonitorParticlePtPDG.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnTPCNcls.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnDEtaDPhi.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnGammaMonitor.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoQinvCorrFctnEMCIC.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctn3DSphericalEMCIC.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoBPLCMS3DCorrFctnEMCIC.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoPairCutAntiGamma.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCutMonitorParticleEtCorr.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoPairCutRadialDistance.cxx)
+set ( SRCS  FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityPairCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityKTPairCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityTPCEntranceSepPairCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoESDTrackCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoTPCInnerCorrFctn.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoChi2CorrFctn.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnSource.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnDEtaDPhi.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelGausRinvFreezeOutGenerator.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnTrueQ.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnNonIdDR.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctn3DSpherical.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctn3DLCMSSpherical.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCutMonitorParticleMomRes.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelAllHiddenInfo.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCutMonitorTrackTPCchiNdof.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoQATrackCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoQAEventCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCutMonitorTrackTPCncls.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityTPCEntranceSepQAPairCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityQAPairCut.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoShareQualityCorrFctn.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoYlm.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnDirectYlm.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnDirectYlm.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCutMonitorParticlePtPDG.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnTPCNcls.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnDEtaDPhi.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnGammaMonitor.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoQinvCorrFctnEMCIC.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctn3DSphericalEMCIC.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoBPLCMS3DCorrFctnEMCIC.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoPairCutAntiGamma.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoCutMonitorParticleEtCorr.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoPairCutRadialDistance.cxx FEMTOSCOPY/AliFemtoUser/AliFemtoPairCutPt.cxx)
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
 
index c9de3f5..b54e5b4 100644 (file)
-//------------------------------------------------------\r
-// AliAnalysisTaskFemto - A task for the analysis framework\r
-// from the FEMTOSCOPY analysis of PWG2. Creates the necessary\r
-// connection between the ESD or AOD input and the femtoscopic\r
-// code.\r
-// Author: Adam Kisiel, OSU; Adam.Kisiel@cern.ch\r
-//------------------------------------------------------\r
-#include "TROOT.h"\r
-#include "TChain.h"\r
-#include "TH1.h"\r
-#include "TCanvas.h"\r
-#include "TSystem.h"\r
-#include "TFile.h"\r
-#include "TInterpreter.h"\r
-\r
-#include "AliAnalysisTask.h"\r
-\r
-#include "AliESDEvent.h"\r
-\r
-#include "AliFemtoAnalysis.h"\r
-#include "AliAnalysisTaskFemto.h"\r
-#include "AliVHeader.h"\r
-#include "AliGenEventHeader.h"\r
-#include "AliGenHijingEventHeader.h"\r
-#include "AliGenCocktailEventHeader.h"\r
-\r
-ClassImp(AliAnalysisTaskFemto)\r
-\r
-// Default name for the setup macro of femto analysis  \r
-// This function MUST be defined in the separate file !!!\r
-// extern AliFemtoManager *ConfigFemtoAnalysis();\r
-\r
-//________________________________________________________________________\r
-AliAnalysisTaskFemto::AliAnalysisTaskFemto(const char *name, const char *aConfigMacro, const char *aConfigParams):\r
-    AliAnalysisTask(name,""), \r
-    fESD(0), \r
-    fESDpid(0),\r
-    fAOD(0),\r
-    fAODpidUtil(0),\r
-    fStack(0),\r
-    fOutputList(0), \r
-    fReader(0x0),\r
-    fManager(0x0),\r
-    fAnalysisType(0),\r
-    fConfigMacro(0),\r
-    fConfigParams(0)\r
-{\r
-  // Constructor.\r
-  // Input slot #0 works with an Ntuple\r
-  DefineInput(0, TChain::Class());\r
-  // Output slot #0 writes into a TH1 container\r
-  DefineOutput(0, TList::Class());\r
-  fConfigMacro = (char *) malloc(sizeof(char) * strlen(aConfigMacro));\r
-  strcpy(fConfigMacro, aConfigMacro);\r
-  fConfigParams = (char *) malloc(sizeof(char) * strlen(aConfigParams));\r
-  strcpy(fConfigParams, aConfigParams);\r
-}\r
-//________________________________________________________________________\r
-AliAnalysisTaskFemto::AliAnalysisTaskFemto(const char *name, const char *aConfigMacro="ConfigFemtoAnalysis.C"): \r
-    AliAnalysisTask(name,""), \r
-    fESD(0), \r
-    fESDpid(0),\r
-    fAOD(0),\r
-    fAODpidUtil(0),\r
-    fStack(0),\r
-    fOutputList(0), \r
-    fReader(0x0),\r
-    fManager(0x0),\r
-    fAnalysisType(0),\r
-    fConfigMacro(0),\r
-    fConfigParams(0)\r
-{\r
-  // Constructor.\r
-  // Input slot #0 works with an Ntuple\r
-  DefineInput(0, TChain::Class());\r
-  // Output slot #0 writes into a TH1 container\r
-  DefineOutput(0, TList::Class());\r
-  fConfigMacro = (char *) malloc(sizeof(char) * strlen(aConfigMacro));\r
-  strcpy(fConfigMacro, aConfigMacro);\r
-  fConfigParams = (char *) malloc(sizeof(char) * 2);\r
-  strcpy(fConfigParams, "");\r
-}\r
-\r
-AliAnalysisTaskFemto::AliAnalysisTaskFemto(const AliAnalysisTaskFemto& aFemtoTask):\r
-    AliAnalysisTask(aFemtoTask), \r
-    fESD(0), \r
-    fESDpid(0),\r
-    fAOD(0),\r
-    fAODpidUtil(0),\r
-    fStack(0),\r
-    fOutputList(0), \r
-    fReader(0x0),\r
-    fManager(0x0),\r
-    fAnalysisType(0),\r
-    fConfigMacro(0),\r
-    fConfigParams(0)\r
-{\r
-  // copy constructor\r
-  fESD = aFemtoTask.fESD; \r
-  fESDpid = aFemtoTask.fESDpid; \r
-  fAOD = aFemtoTask.fAOD; \r
-  fAODpidUtil = aFemtoTask.fAODpidUtil;\r
-  fStack = aFemtoTask.fStack;\r
-  fOutputList = aFemtoTask.fOutputList;   \r
-  fReader = aFemtoTask.fReader;       \r
-  fManager = aFemtoTask.fManager;      \r
-  fAnalysisType = aFemtoTask.fAnalysisType; \r
-  fConfigMacro = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigMacro));\r
-  strcpy(fConfigMacro, aFemtoTask.fConfigMacro);\r
-  fConfigParams = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigParams));\r
-  strcpy(fConfigParams, aFemtoTask.fConfigParams);\r
-}\r
-\r
-\r
-AliAnalysisTaskFemto& AliAnalysisTaskFemto::operator=(const AliAnalysisTaskFemto& aFemtoTask){\r
-  // assignment operator\r
-  if (this == &aFemtoTask)\r
-    return *this;\r
-\r
-  fESD = aFemtoTask.fESD; \r
-  fESDpid = aFemtoTask.fESDpid;\r
-  fAOD = aFemtoTask.fAOD; \r
-  fAODpidUtil = aFemtoTask.fAODpidUtil;\r
-  fStack = aFemtoTask.fStack;\r
-  fOutputList = aFemtoTask.fOutputList;   \r
-  fReader = aFemtoTask.fReader;       \r
-  fManager = aFemtoTask.fManager;      \r
-  fAnalysisType = aFemtoTask.fAnalysisType; \r
-  if (fConfigMacro) free(fConfigMacro);\r
-  fConfigMacro = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigMacro));\r
-  strcpy(fConfigMacro, aFemtoTask.fConfigMacro);\r
-  if (fConfigParams) free(fConfigParams);\r
-  fConfigParams = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigParams));\r
-  strcpy(fConfigParams, aFemtoTask.fConfigParams);\r
-\r
-  return *this;\r
-}\r
-\r
-AliAnalysisTaskFemto::~AliAnalysisTaskFemto() \r
-{\r
-  if (fConfigMacro) free(fConfigMacro);\r
-  if (fConfigParams) free(fConfigParams);\r
-}\r
-\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskFemto::ConnectInputData(Option_t *) {\r
-  AliInfo(Form("   ConnectInputData %s\n", GetName()));\r
-\r
-  fESD = 0;\r
-  fESDpid = 0;\r
-  fAOD = 0;\r
-  fAODpidUtil = 0;\r
-  fAnalysisType = 0;\r
-\r
-  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));\r
-  if (!tree) {\r
-    AliWarning("Could not read chain from input slot 0");\r
-    return;\r
-  } \r
-\r
-  AliFemtoEventReaderESDChain *femtoReader = dynamic_cast<AliFemtoEventReaderESDChain *> (fReader);\r
-  if ((dynamic_cast<AliFemtoEventReaderESDChain *> (fReader)) ||\r
-      (dynamic_cast<AliFemtoEventReaderESDChainKine *> (fReader))) {\r
-    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
-    \r
-    if(esdH) {\r
-      AliInfo("Selected ESD analysis");\r
-      fAnalysisType = 1;\r
-      \r
-//       if (!esdH) {\r
-//     AliWarning("Could not get ESDInputHandler");\r
-//       } \r
-//       else {\r
-       fESD = esdH->GetEvent();\r
-        fESDpid = esdH->GetESDpid();\r
-        femtoReader->SetESDPid(fESDpid);\r
-//       }\r
-    }\r
-  }\r
-  \r
-    AliFemtoEventReaderAODChain *femtoReaderAOD = dynamic_cast<AliFemtoEventReaderAODChain *> (fReader);\r
-  if (dynamic_cast<AliFemtoEventReaderAODChain *> (fReader)) {\r
-    AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
-       \r
-    if (!aodH) {\r
-      TObject *handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();\r
-      AliInfo("Has output handler ");\r
-      if( handler && handler->InheritsFrom("AliAODHandler") ) {\r
-       AliInfo("Selected AOD analysis");\r
-\r
-       fAOD = ((AliAODHandler*)handler)->GetAOD();\r
-       fAnalysisType = 2;\r
-      }\r
-      else {\r
-       AliWarning("Selected AOD reader but no AOD handler found");\r
-      }\r
-    } \r
-    else {\r
-      AliInfo("Selected AOD analysis");\r
-      fAnalysisType = 2;\r
-      \r
-      fAOD = aodH->GetEvent();\r
-\r
-      fAODpidUtil = aodH->GetAODpidUtil();\r
-      //      printf("aodH->GetAODpidUtil(): %x",aodH->GetAODpidUtil());\r
-      femtoReaderAOD->SetAODpidUtil(fAODpidUtil);\r
-    }\r
-  }\r
-\r
-  if ((!fAOD) && (!fESD)) {\r
-    AliWarning("Wrong analysis type: Only ESD and AOD types are allowed!");\r
-  }\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskFemto::CreateOutputObjects() {\r
-  AliInfo("Creating Femto Analysis objects\n");\r
-\r
-  gSystem->SetIncludePath("-I$ROOTSYS/include -I./STEERBase/ -I./ESD/ -I./AOD/ -I./ANALYSIS/ -I./ANALYSISalice/ -I./PWG2AOD/AOD -I./PWG2femtoscopy/FEMTOSCOPY/AliFemto -I./PWG2femtoscopyUser/FEMTOSCOPY/AliFemtoUser");\r
-  //  char fcm[2000];\r
-//   sprintf(fcm, "%s++", fConfigMacro);\r
-//   gROOT->LoadMacro(fcm);\r
-  gROOT->LoadMacro(fConfigMacro);\r
-  //  fJetFinder = (AliJetFinder*) gInterpreter->ProcessLine("ConfigJetAnalysis()");\r
-  if (!fConfigParams)\r
-    SetFemtoManager((AliFemtoManager *) gInterpreter->ProcessLine("ConfigFemtoAnalysis()"));\r
-  else\r
-    SetFemtoManager((AliFemtoManager *) gInterpreter->ProcessLine(Form("ConfigFemtoAnalysis(%s)", fConfigParams)));\r
-\r
-  TList *tOL;\r
-  fOutputList = fManager->Analysis(0)->GetOutputList();\r
-\r
-  for (unsigned int ian = 1; ian<fManager->AnalysisCollection()->size(); ian++) {\r
-    tOL = fManager->Analysis(ian)->GetOutputList();\r
-\r
-    TIter nextListCf(tOL);\r
-    while (TObject *obj = nextListCf()) {\r
-      fOutputList->Add(obj);\r
-    }\r
-\r
-    delete tOL;\r
-  }\r
-\r
-  PostData(0, fOutputList);\r
-}\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskFemto::Exec(Option_t *) {\r
-  // Task making a femtoscopic analysis.\r
-\r
-  if (fAnalysisType==1) {\r
-    if (!fESD) {\r
-      AliWarning("fESD not available");\r
-      return;\r
-    }\r
-\r
-    //Get MC data\r
-    AliMCEventHandler*    mctruth = (AliMCEventHandler*) \r
-      ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());\r
-    \r
-    AliGenHijingEventHeader *hdh = 0;\r
-    if(mctruth) {\r
-      fStack = mctruth->MCEvent()->Stack();\r
-\r
-      AliGenCocktailEventHeader *hd = dynamic_cast<AliGenCocktailEventHeader *> (mctruth->MCEvent()->GenEventHeader());\r
-      \r
-      if (hd) {\r
-       \r
-       //      AliInfo ("Got MC cocktail event header %p\n", (void *) hd);\r
-       TList *lhd = hd->GetHeaders();\r
-       //      AliInfo ("Got list of headers %d\n", lhd->GetEntries());\r
-       \r
-       for (int iterh=0; iterh<lhd->GetEntries(); iterh++) \r
-         {\r
-           hdh = dynamic_cast<AliGenHijingEventHeader *> (lhd->At(iterh));\r
-           //      AliInfo ("HIJING header at %i is %p\n", iterh, (void *) hdh);\r
-         }\r
-      }    \r
-    }\r
-\r
-    // Get ESD\r
-    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
-    \r
-    if (!esdH) {\r
-      AliWarning("Could not get ESDInputHandler");\r
-      return;\r
-    } \r
-    else {\r
-      fESD = esdH->GetEvent();\r
-      fESDpid = esdH->GetESDpid();   \r
-    }\r
-\r
-    AliInfo(Form("Tracks in ESD: %d \n",fESD->GetNumberOfTracks()));\r
-\r
-    if (fESD->GetNumberOfTracks() >= 0) {\r
-    \r
-      if (!fReader) {\r
-       AliWarning("No ESD reader for ESD analysis !\n");\r
-      }\r
-      \r
-      AliFemtoEventReaderESDChain* fesdc = dynamic_cast<AliFemtoEventReaderESDChain *> (fReader);\r
-      if (fesdc)\r
-       {\r
-         // Process the event with no Kine information\r
-         fesdc->SetESDSource(fESD);\r
-         fManager->ProcessEvent();\r
-       }\r
-      AliFemtoEventReaderESDChainKine* fesdck = dynamic_cast<AliFemtoEventReaderESDChainKine *> (fReader);\r
-      if (fesdck) \r
-       {\r
-         // Process the event with Kine information\r
-         fesdck->SetESDSource(fESD);\r
-         fesdck->SetStackSource(fStack);\r
-         \r
-         fesdck->SetGenEventHeader(hdh);\r
-         fManager->ProcessEvent();\r
-       }\r
-      AliFemtoEventReaderStandard* fstd = dynamic_cast<AliFemtoEventReaderStandard *> (fReader);\r
-      if (fstd) \r
-       {\r
-         // Process the event with Kine information\r
-         fstd->SetESDSource(fESD);\r
-         if (mctruth) {\r
-           fstd->SetStackSource(fStack);\r
-           fstd->SetGenEventHeader(hdh);\r
-           fstd->SetInputType(AliFemtoEventReaderStandard::kESDKine);\r
-         }\r
-         else\r
-           fstd->SetInputType(AliFemtoEventReaderStandard::kESD);\r
-         fManager->ProcessEvent();\r
-       }\r
-    } \r
-\r
-    // Post the output histogram list\r
-    PostData(0, fOutputList);\r
-  }\r
-\r
-  if (fAnalysisType==2) {    \r
-    if (!fAOD) {\r
-      AliWarning("fAOD not available");\r
-      return;\r
-    }\r
-\r
-    // Get AOD\r
-//     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());\r
-      \r
-//     if (!aodH) {\r
-//       AliWarning("Could not get AODInputHandler");\r
-//       return;\r
-//     } \r
-//     else {\r
-\r
-//       fAOD = aodH->GetEvent();\r
-//     }\r
-\r
-    AliInfo(Form("Tracks in AOD: %d \n",fAOD->GetNumberOfTracks()));\r
-    \r
-    if (fAOD->GetNumberOfTracks() > 0) {\r
-      if (!fReader) {\r
-       AliWarning("No AOD reader for AOD analysis! \n");\r
-      }\r
-      else {\r
-       AliFemtoEventReaderAODChain* faodc = dynamic_cast<AliFemtoEventReaderAODChain *> (fReader);\r
-\r
-       if (faodc) {\r
-         // Process the event\r
-         faodc->SetAODSource(fAOD);\r
-         fManager->ProcessEvent();\r
-       }\r
-       AliFemtoEventReaderStandard* fstd = dynamic_cast<AliFemtoEventReaderStandard *> (fReader);\r
-\r
-       if (fstd) {\r
-         // Process the event\r
-         fstd->SetAODSource(fAOD);\r
-         fstd->SetInputType(AliFemtoEventReaderStandard::kAOD);\r
-         fManager->ProcessEvent();\r
-       }\r
-      }\r
-    } \r
-\r
-    // Post the output histogram list\r
-    PostData(0, fOutputList);\r
-  }\r
-}      \r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskFemto::Terminate(Option_t *) {\r
-  // Do the final processing\r
-  if (fManager) {\r
-    fManager->Finish();\r
-  }\r
-}\r
-//________________________________________________________________________\r
-void AliAnalysisTaskFemto:: FinishTaskOutput() {\r
-  // Do the final processing\r
-  if (fManager) {\r
-    fManager->Finish();\r
-  }\r
-}\r
-//________________________________________________________________________\r
-void AliAnalysisTaskFemto::SetFemtoReaderESD(AliFemtoEventReaderESDChain *aReader)\r
-{\r
-  AliInfo("Selecting Femto reader for ESD\n");\r
-  fReader = aReader;\r
-}\r
-//________________________________________________________________________\r
-void AliAnalysisTaskFemto::SetFemtoReaderESDKine(AliFemtoEventReaderESDChainKine *aReader)\r
-{\r
-  AliInfo("Selecting Femto reader for ESD with Kinematics information\n");\r
-  fReader = aReader;\r
-}\r
-//________________________________________________________________________\r
-void AliAnalysisTaskFemto::SetFemtoReaderAOD(AliFemtoEventReaderAODChain *aReader)\r
-{\r
-  AliInfo("Selecting Femto reader for AOD\n");\r
-  fReader = aReader;\r
-}\r
-void AliAnalysisTaskFemto::SetFemtoReaderStandard(AliFemtoEventReaderStandard *aReader)\r
-{\r
-  AliInfo("Selecting Standard all-purpose Femto reader\n");\r
-  fReader = aReader;\r
-}\r
-//________________________________________________________________________\r
-void AliAnalysisTaskFemto::SetFemtoManager(AliFemtoManager *aManager)\r
-{\r
-  fManager = aManager;\r
-  AliInfo(Form("Got reader %p\n", (void *) aManager->EventReader()));\r
-  AliFemtoEventReaderESDChain     *tReaderESDChain     = dynamic_cast<AliFemtoEventReaderESDChain *> (aManager->EventReader());\r
-  AliFemtoEventReaderESDChainKine *tReaderESDChainKine = dynamic_cast<AliFemtoEventReaderESDChainKine *> (aManager->EventReader());\r
-  AliFemtoEventReaderAODChain     *tReaderAODChain     = dynamic_cast<AliFemtoEventReaderAODChain *> (aManager->EventReader());\r
-  AliFemtoEventReaderStandard     *tReaderStandard     = dynamic_cast<AliFemtoEventReaderStandard *> (aManager->EventReader());\r
-\r
-  if ((!tReaderESDChain) && (!tReaderESDChainKine) && (!tReaderAODChain) && (!tReaderStandard)) {\r
-    AliWarning("No AliFemto event reader created. Will not run femto analysis.\n");\r
-    return;\r
-  }\r
-  if (tReaderESDChain) SetFemtoReaderESD(tReaderESDChain);\r
-  if (tReaderESDChainKine) SetFemtoReaderESDKine(tReaderESDChainKine);\r
-  if (tReaderAODChain) SetFemtoReaderAOD(tReaderAODChain);\r
-  if (tReaderStandard) SetFemtoReaderStandard(tReaderStandard);\r
-}\r
-\r
+//------------------------------------------------------
+// AliAnalysisTaskFemto - A task for the analysis framework
+// from the FEMTOSCOPY analysis of PWG2. Creates the necessary
+// connection between the ESD or AOD input and the femtoscopic
+// code.
+// Author: Adam Kisiel, OSU; Adam.Kisiel@cern.ch
+//------------------------------------------------------
+#include "TROOT.h"
+#include "TChain.h"
+#include "TH1.h"
+#include "TCanvas.h"
+#include "TSystem.h"
+#include "TFile.h"
+#include "TInterpreter.h"
+
+#include "AliAnalysisTask.h"
+
+#include "AliESDEvent.h"
+
+#include "AliFemtoAnalysis.h"
+#include "AliAnalysisTaskFemto.h"
+#include "AliVHeader.h"
+#include "AliGenEventHeader.h"
+#include "AliGenHijingEventHeader.h"
+#include "AliGenCocktailEventHeader.h"
+
+ClassImp(AliAnalysisTaskFemto)
+
+// Default name for the setup macro of femto analysis  
+// This function MUST be defined in the separate file !!!
+// extern AliFemtoManager *ConfigFemtoAnalysis();
+
+//________________________________________________________________________
+AliAnalysisTaskFemto::AliAnalysisTaskFemto(const char *name, const char *aConfigMacro, const char *aConfigParams):
+    AliAnalysisTask(name,""), 
+    fESD(0), 
+    fESDpid(0),
+    fAOD(0),
+    fAODpidUtil(0),
+    fStack(0),
+    fOutputList(0), 
+    fReader(0x0),
+    fManager(0x0),
+    fAnalysisType(0),
+    fConfigMacro(0),
+    fConfigParams(0)
+{
+  // Constructor.
+  // Input slot #0 works with an Ntuple
+  DefineInput(0, TChain::Class());
+  // Output slot #0 writes into a TH1 container
+  DefineOutput(0, TList::Class());
+  fConfigMacro = (char *) malloc(sizeof(char) * strlen(aConfigMacro));
+  strcpy(fConfigMacro, aConfigMacro);
+  fConfigParams = (char *) malloc(sizeof(char) * strlen(aConfigParams));
+  strcpy(fConfigParams, aConfigParams);
+}
+//________________________________________________________________________
+AliAnalysisTaskFemto::AliAnalysisTaskFemto(const char *name, const char *aConfigMacro="ConfigFemtoAnalysis.C"): 
+    AliAnalysisTask(name,""), 
+    fESD(0), 
+    fESDpid(0),
+    fAOD(0),
+    fAODpidUtil(0),
+    fStack(0),
+    fOutputList(0), 
+    fReader(0x0),
+    fManager(0x0),
+    fAnalysisType(0),
+    fConfigMacro(0),
+    fConfigParams(0)
+{
+  // Constructor.
+  // Input slot #0 works with an Ntuple
+  DefineInput(0, TChain::Class());
+  // Output slot #0 writes into a TH1 container
+  DefineOutput(0, TList::Class());
+  fConfigMacro = (char *) malloc(sizeof(char) * strlen(aConfigMacro));
+  strcpy(fConfigMacro, aConfigMacro);
+  fConfigParams = (char *) malloc(sizeof(char) * 2);
+  strcpy(fConfigParams, "");
+}
+
+AliAnalysisTaskFemto::AliAnalysisTaskFemto(const AliAnalysisTaskFemto& aFemtoTask):
+    AliAnalysisTask(aFemtoTask), 
+    fESD(0), 
+    fESDpid(0),
+    fAOD(0),
+    fAODpidUtil(0),
+    fStack(0),
+    fOutputList(0), 
+    fReader(0x0),
+    fManager(0x0),
+    fAnalysisType(0),
+    fConfigMacro(0),
+    fConfigParams(0)
+{
+  // copy constructor
+  fESD = aFemtoTask.fESD; 
+  fESDpid = aFemtoTask.fESDpid; 
+  fAOD = aFemtoTask.fAOD; 
+  fAODpidUtil = aFemtoTask.fAODpidUtil;
+  fStack = aFemtoTask.fStack;
+  fOutputList = aFemtoTask.fOutputList;   
+  fReader = aFemtoTask.fReader;       
+  fManager = aFemtoTask.fManager;      
+  fAnalysisType = aFemtoTask.fAnalysisType; 
+  fConfigMacro = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigMacro));
+  strcpy(fConfigMacro, aFemtoTask.fConfigMacro);
+  fConfigParams = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigParams));
+  strcpy(fConfigParams, aFemtoTask.fConfigParams);
+}
+
+
+AliAnalysisTaskFemto& AliAnalysisTaskFemto::operator=(const AliAnalysisTaskFemto& aFemtoTask){
+  // assignment operator
+  if (this == &aFemtoTask)
+    return *this;
+
+  fESD = aFemtoTask.fESD; 
+  fESDpid = aFemtoTask.fESDpid;
+  fAOD = aFemtoTask.fAOD; 
+  fAODpidUtil = aFemtoTask.fAODpidUtil;
+  fStack = aFemtoTask.fStack;
+  fOutputList = aFemtoTask.fOutputList;   
+  fReader = aFemtoTask.fReader;       
+  fManager = aFemtoTask.fManager;      
+  fAnalysisType = aFemtoTask.fAnalysisType; 
+  if (fConfigMacro) free(fConfigMacro);
+  fConfigMacro = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigMacro));
+  strcpy(fConfigMacro, aFemtoTask.fConfigMacro);
+  if (fConfigParams) free(fConfigParams);
+  fConfigParams = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigParams));
+  strcpy(fConfigParams, aFemtoTask.fConfigParams);
+
+  return *this;
+}
+
+AliAnalysisTaskFemto::~AliAnalysisTaskFemto() 
+{
+  if (fConfigMacro) free(fConfigMacro);
+  if (fConfigParams) free(fConfigParams);
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskFemto::ConnectInputData(Option_t *) {
+  AliInfo(Form("   ConnectInputData %s\n", GetName()));
+
+  fESD = 0;
+  fESDpid = 0;
+  fAOD = 0;
+  fAODpidUtil = 0;
+  fAnalysisType = 0;
+
+  TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
+  if (!tree) {
+    AliWarning("Could not read chain from input slot 0");
+    return;
+  } 
+
+  AliFemtoEventReaderESDChain *femtoReader = dynamic_cast<AliFemtoEventReaderESDChain *> (fReader);
+  if ((dynamic_cast<AliFemtoEventReaderESDChain *> (fReader))) {
+    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    
+    if(esdH) {
+      AliInfo("Selected ESD analysis");
+      fAnalysisType = 1;
+      
+//       if (!esdH) {
+//     AliWarning("Could not get ESDInputHandler");
+//       } 
+//       else {
+       fESD = esdH->GetEvent();
+        fESDpid = esdH->GetESDpid();
+        femtoReader->SetESDPid(fESDpid);
+//       }
+    }
+}
+  else if ((dynamic_cast<AliFemtoEventReaderKinematicsChain *> (fReader))) {
+    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    
+    if(esdH) {
+      AliInfo("Selected ESD analysis");
+      fAnalysisType = 1;
+      
+//       if (!esdH) {
+//     AliWarning("Could not get ESDInputHandler");
+//       } 
+//       else {
+       fESD = esdH->GetEvent();
+        //fESDpid = esdH->GetESDpid();
+        //femtoReader->SetESDPid(fESDpid);
+//       }
+    }
+  }
+
+
+  AliFemtoEventReaderESDChainKine *femtoReaderESDKine = dynamic_cast<AliFemtoEventReaderESDChainKine *> (fReader);
+  if ((dynamic_cast<AliFemtoEventReaderESDChainKine *> (fReader))) {
+    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    
+    if(esdH) {
+      AliInfo("Selected ESD analysis");
+      fAnalysisType = 1;
+      
+//       if (!esdH) {
+//     AliWarning("Could not get ESDInputHandler");
+//       } 
+//       else {
+       fESD = esdH->GetEvent();
+        fESDpid = esdH->GetESDpid();
+        femtoReaderESDKine->SetESDPid(fESDpid);
+//       }
+    }
+}
+
+  AliFemtoEventReaderKinematicsChain *femtoReaderKine = dynamic_cast<AliFemtoEventReaderKinematicsChain *> (fReader);
+if ((dynamic_cast<AliFemtoEventReaderKinematicsChain *> (fReader))) {
+    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    
+    if(esdH) {
+      AliInfo("Selected ESD analysis");
+      fAnalysisType = 1;
+      
+//       if (!esdH) {
+//     AliWarning("Could not get ESDInputHandler");
+//       } 
+//       else {
+       fESD = esdH->GetEvent();
+        //fESDpid = esdH->GetESDpid();
+        //femtoReader->SetESDPid(fESDpid);
+//       }
+    }
+ }
+
+  
+    AliFemtoEventReaderAODChain *femtoReaderAOD = dynamic_cast<AliFemtoEventReaderAODChain *> (fReader);
+  if (dynamic_cast<AliFemtoEventReaderAODChain *> (fReader)) {
+    AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+       
+    if (!aodH) {
+      TObject *handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
+      AliInfo("Has output handler ");
+      if( handler && handler->InheritsFrom("AliAODHandler") ) {
+       AliInfo("Selected AOD analysis");
+
+       fAOD = ((AliAODHandler*)handler)->GetAOD();
+       fAnalysisType = 2;
+      }
+      else {
+       AliWarning("Selected AOD reader but no AOD handler found");
+      }
+    } 
+    else {
+      AliInfo("Selected AOD analysis");
+      fAnalysisType = 2;
+      
+      fAOD = aodH->GetEvent();
+
+      fAODpidUtil = aodH->GetAODpidUtil();
+      //      printf("aodH->GetAODpidUtil(): %x",aodH->GetAODpidUtil());
+      femtoReaderAOD->SetAODpidUtil(fAODpidUtil);
+    }
+  }
+
+  if ((!fAOD) && (!fESD)) {
+    AliWarning("Wrong analysis type: Only ESD and AOD types are allowed!");
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskFemto::CreateOutputObjects() {
+  AliInfo("Creating Femto Analysis objects\n");
+
+  gSystem->SetIncludePath("-I$ROOTSYS/include -I./STEERBase/ -I./ESD/ -I./AOD/ -I./ANALYSIS/ -I./ANALYSISalice/ -I./PWG2AOD/AOD -I./PWG2femtoscopy/FEMTOSCOPY/AliFemto -I./PWG2femtoscopyUser/FEMTOSCOPY/AliFemtoUser");
+  //  char fcm[2000];
+//   sprintf(fcm, "%s++", fConfigMacro);
+//   gROOT->LoadMacro(fcm);
+  gROOT->LoadMacro(fConfigMacro);
+  //  fJetFinder = (AliJetFinder*) gInterpreter->ProcessLine("ConfigJetAnalysis()");
+  if (!fConfigParams)
+    SetFemtoManager((AliFemtoManager *) gInterpreter->ProcessLine("ConfigFemtoAnalysis()"));
+  else
+    SetFemtoManager((AliFemtoManager *) gInterpreter->ProcessLine(Form("ConfigFemtoAnalysis(%s)", fConfigParams)));
+
+  TList *tOL;
+  fOutputList = fManager->Analysis(0)->GetOutputList();
+
+  for (unsigned int ian = 1; ian<fManager->AnalysisCollection()->size(); ian++) {
+    tOL = fManager->Analysis(ian)->GetOutputList();
+
+    TIter nextListCf(tOL);
+    while (TObject *obj = nextListCf()) {
+      fOutputList->Add(obj);
+    }
+
+    delete tOL;
+  }
+
+  PostData(0, fOutputList);
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskFemto::Exec(Option_t *) {
+  // Task making a femtoscopic analysis.
+
+  if (fAnalysisType==1) {
+    if (!fESD) {
+      AliWarning("fESD not available");
+      return;
+    }
+
+    //Get MC data
+    AliMCEventHandler*    mctruth = (AliMCEventHandler*) 
+      ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
+    
+    AliGenHijingEventHeader *hdh = 0;
+    if(mctruth) {
+      fStack = mctruth->MCEvent()->Stack();
+
+      AliGenCocktailEventHeader *hd = dynamic_cast<AliGenCocktailEventHeader *> (mctruth->MCEvent()->GenEventHeader());
+      
+      if (hd) {
+       
+       //      AliInfo ("Got MC cocktail event header %p\n", (void *) hd);
+       TList *lhd = hd->GetHeaders();
+       //      AliInfo ("Got list of headers %d\n", lhd->GetEntries());
+       
+       for (int iterh=0; iterh<lhd->GetEntries(); iterh++) 
+         {
+           hdh = dynamic_cast<AliGenHijingEventHeader *> (lhd->At(iterh));
+           //      AliInfo ("HIJING header at %i is %p\n", iterh, (void *) hdh);
+         }
+      }    
+    }
+
+    // Get ESD
+    AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    
+    if (!esdH) {
+      AliWarning("Could not get ESDInputHandler");
+      return;
+    } 
+    else {
+      fESD = esdH->GetEvent();
+      fESDpid = esdH->GetESDpid();   
+    }
+
+    AliInfo(Form("Tracks in ESD: %d \n",fESD->GetNumberOfTracks()));
+
+    if (fESD->GetNumberOfTracks() >= 0) {
+    
+      if (!fReader) {
+       AliWarning("No ESD reader for ESD analysis !\n");
+      }
+      
+      AliFemtoEventReaderESDChain* fesdc = dynamic_cast<AliFemtoEventReaderESDChain *> (fReader);
+      if (fesdc)
+       {
+         // Process the event with no Kine information
+         fesdc->SetESDSource(fESD);
+         fManager->ProcessEvent();
+       }
+
+       AliFemtoEventReaderKinematicsChain* fkinec = dynamic_cast<AliFemtoEventReaderKinematicsChain *> (fReader);
+      if (fkinec)
+       {
+         // Process the event with Kine information only
+         fkinec->SetStackSource(fStack);
+         fManager->ProcessEvent();
+       }
+
+
+      AliFemtoEventReaderESDChainKine* fesdck = dynamic_cast<AliFemtoEventReaderESDChainKine *> (fReader);
+      if (fesdck) 
+       {
+         // Process the event with Kine information
+         fesdck->SetESDSource(fESD);
+         fesdck->SetStackSource(fStack);
+         
+         fesdck->SetGenEventHeader(hdh);
+         fManager->ProcessEvent();
+       }
+      AliFemtoEventReaderStandard* fstd = dynamic_cast<AliFemtoEventReaderStandard *> (fReader);
+      if (fstd) 
+       {
+         // Process the event with Kine information
+         fstd->SetESDSource(fESD);
+         if (mctruth) {
+           fstd->SetStackSource(fStack);
+           fstd->SetGenEventHeader(hdh);
+           fstd->SetInputType(AliFemtoEventReaderStandard::kESDKine);
+         }
+         else
+           fstd->SetInputType(AliFemtoEventReaderStandard::kESD);
+         fManager->ProcessEvent();
+       }
+    } 
+
+    // Post the output histogram list
+    PostData(0, fOutputList);
+  }
+
+  if (fAnalysisType==2) {    
+    if (!fAOD) {
+      AliWarning("fAOD not available");
+      return;
+    }
+
+    // Get AOD
+//     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+      
+//     if (!aodH) {
+//       AliWarning("Could not get AODInputHandler");
+//       return;
+//     } 
+//     else {
+
+//       fAOD = aodH->GetEvent();
+//     }
+
+    AliInfo(Form("Tracks in AOD: %d \n",fAOD->GetNumberOfTracks()));
+    
+    if (fAOD->GetNumberOfTracks() > 0) {
+      if (!fReader) {
+       AliWarning("No AOD reader for AOD analysis! \n");
+      }
+      else {
+       AliFemtoEventReaderAODChain* faodc = dynamic_cast<AliFemtoEventReaderAODChain *> (fReader);
+
+       if (faodc) {
+         // Process the event
+         faodc->SetAODSource(fAOD);
+         fManager->ProcessEvent();
+       }
+       AliFemtoEventReaderStandard* fstd = dynamic_cast<AliFemtoEventReaderStandard *> (fReader);
+
+       if (fstd) {
+         // Process the event
+         fstd->SetAODSource(fAOD);
+         fstd->SetInputType(AliFemtoEventReaderStandard::kAOD);
+         fManager->ProcessEvent();
+       }
+      }
+    } 
+
+    // Post the output histogram list
+    PostData(0, fOutputList);
+  }
+}      
+
+//________________________________________________________________________
+void AliAnalysisTaskFemto::Terminate(Option_t *) {
+  // Do the final processing
+  if (fManager) {
+    fManager->Finish();
+  }
+}
+//________________________________________________________________________
+void AliAnalysisTaskFemto:: FinishTaskOutput() {
+  // Do the final processing
+  if (fManager) {
+    fManager->Finish();
+  }
+}
+//________________________________________________________________________
+void AliAnalysisTaskFemto::SetFemtoReaderESD(AliFemtoEventReaderESDChain *aReader)
+{
+  AliInfo("Selecting Femto reader for ESD\n");
+  fReader = aReader;
+}
+//________________________________________________________________________
+void AliAnalysisTaskFemto::SetFemtoReaderESDKine(AliFemtoEventReaderESDChainKine *aReader)
+{
+  AliInfo("Selecting Femto reader for ESD with Kinematics information\n");
+  fReader = aReader;
+}
+//________________________________________________________________________
+void AliAnalysisTaskFemto::SetFemtoReaderAOD(AliFemtoEventReaderAODChain *aReader)
+{
+  AliInfo("Selecting Femto reader for AOD\n");
+  fReader = aReader;
+}
+void AliAnalysisTaskFemto::SetFemtoReaderStandard(AliFemtoEventReaderStandard *aReader)
+{
+  AliInfo("Selecting Standard all-purpose Femto reader\n");
+  fReader = aReader;
+}
+void AliAnalysisTaskFemto::SetFemtoReaderKinematics(AliFemtoEventReaderKinematicsChain *aReader)
+{
+  printf("Selecting Femto reader for Kinematics (Monte Carlo) information\n");
+  fReader = aReader;
+}
+//________________________________________________________________________
+void AliAnalysisTaskFemto::SetFemtoManager(AliFemtoManager *aManager)
+{
+  fManager = aManager;
+  AliInfo(Form("Got reader %p\n", (void *) aManager->EventReader()));
+  AliFemtoEventReaderESDChain     *tReaderESDChain     = dynamic_cast<AliFemtoEventReaderESDChain *> (aManager->EventReader());
+  AliFemtoEventReaderESDChainKine *tReaderESDChainKine = dynamic_cast<AliFemtoEventReaderESDChainKine *> (aManager->EventReader());
+  AliFemtoEventReaderAODChain     *tReaderAODChain     = dynamic_cast<AliFemtoEventReaderAODChain *> (aManager->EventReader());
+  AliFemtoEventReaderStandard     *tReaderStandard     = dynamic_cast<AliFemtoEventReaderStandard *> (aManager->EventReader());
+  AliFemtoEventReaderKinematicsChain *tReaderKineChain = dynamic_cast<AliFemtoEventReaderKinematicsChain *> (aManager->EventReader());
+
+ if ((!tReaderESDChain) && (!tReaderESDChainKine) && (!tReaderAODChain) && (!tReaderStandard) && (!tReaderKineChain)) {
+    AliWarning("No AliFemto event reader created. Will not run femto analysis.\n");
+    return;
+  }
+  if (tReaderESDChain) SetFemtoReaderESD(tReaderESDChain);
+  if (tReaderESDChainKine) SetFemtoReaderESDKine(tReaderESDChainKine);
+  if (tReaderAODChain) SetFemtoReaderAOD(tReaderAODChain);
+  if (tReaderStandard) SetFemtoReaderStandard(tReaderStandard);
+  if (tReaderKineChain) SetFemtoReaderKinematics(tReaderKineChain);
+}
+
index a96f2de..711f450 100644 (file)
@@ -26,6 +26,7 @@
 #include "AliFemtoEventReaderESDChainKine.h"
 #include "AliFemtoEventReaderAODChain.h"
 #include "AliFemtoEventReaderStandard.h"
+#include "AliFemtoEventReaderKinematicsChain.h"
 #include "AliFemtoManager.h"
 
 #include "AliESDpid.h"
@@ -53,6 +54,7 @@ class AliAnalysisTaskFemto : public AliAnalysisTask {
   void SetFemtoReaderAOD(AliFemtoEventReaderAODChain *aReader);
   void SetFemtoReaderStandard(AliFemtoEventReaderStandard *aReader);
   void SetFemtoManager(AliFemtoManager *aManager);
+  void SetFemtoReaderKinematics(AliFemtoEventReaderKinematicsChain *aReader);
 
  private:
   AliESDEvent                 *fESD;          //! ESD object
index c5e5497..4af476c 100644 (file)
@@ -24,7 +24,9 @@ AliFemtoCutMonitorEventMult::AliFemtoCutMonitorEventMult():
   fEst2Est3(0),
   fEst1Norm(0),
   fEst2Norm(0),
-  fEst3Norm(0)
+  fEst3Norm(0),
+  freadMC(kFALSE),
+  faddhists(kFALSE)
 {
   // Default constructor
   fEvMult = new TH1D("EvMult", "Event Multiplicity", 5001, -0.5, 5000.5);
@@ -45,7 +47,9 @@ AliFemtoCutMonitorEventMult::AliFemtoCutMonitorEventMult(const char *aName):
   fEst2Est3(0),
   fEst1Norm(0),
   fEst2Norm(0),
-  fEst3Norm(0)
+  fEst3Norm(0),
+  freadMC(kFALSE),
+  faddhists(kFALSE)
 {
   // Normal constructor
   char name[200];
@@ -55,38 +59,43 @@ AliFemtoCutMonitorEventMult::AliFemtoCutMonitorEventMult(const char *aName):
   snprintf(name, 200, "NormEvMult%s", aName);
   fNormEvMult = new TH1D(name, "Normalized Event Multiplicity", 5001, -0.5, 5000.5);
 
-  snprintf(name, 200, "SPDEvMult%s", aName);
-  fSPDMult = new TH1D(name, "SPD Tracklet Multiplicity", 5001, -0.5, 5000.5);
+  if(!freadMC) {
+    snprintf(name, 200, "SPDEvMult%s", aName);
+    fSPDMult = new TH1D(name, "SPD Tracklet Multiplicity", 5001, -0.5, 5000.5);
+  }
 
   snprintf(name, 200, "EvMultTotPt%s", aName);
   fMultSumPt = new TH2D(name,"Event Multiplicity vs Total pT",501,-0.5,500.5,1000,0.0,100.0);
 
-  snprintf(name, 200, "EvMultEstITSTPC%s", aName);
-  fEstimateITSTPC = new TH1D(name, "ITS+TPC Multiplicity Estimate", 5001, -0.5, 5000.5);
+  if(faddhists)
+    {
+      snprintf(name, 200, "EvMultEstITSTPC%s", aName);
+      fEstimateITSTPC = new TH1D(name, "ITS+TPC Multiplicity Estimate", 5001, -0.5, 5000.5);
 
-  snprintf(name, 200, "EvMultEstTracklets%s", aName);
-  fEstimateTracklets = new TH1D(name, "Tracklets Multiplicity Estimate", 5001, -0.5, 5000.5);
+      snprintf(name, 200, "EvMultEstTracklets%s", aName);
+      fEstimateTracklets = new TH1D(name, "Tracklets Multiplicity Estimate", 5001, -0.5, 5000.5);
 
-  snprintf(name, 200, "EvMultEstITSPure%s", aName);
-  fEstimateITSPure = new TH1D(name, "ITS Pure Multiplicity Estimate", 8001, -0.5, 8000.5);
+      snprintf(name, 200, "EvMultEstITSPure%s", aName);
+      fEstimateITSPure = new TH1D(name, "ITS Pure Multiplicity Estimate", 8001, -0.5, 8000.5);
 
-  snprintf(name, 200, "EstITSTPCEstTracklet%s", aName);
-  fEst1Est2 = new TH2D(name,"ITS+TPC vs Tracklets",501,-0.5,5000.5,501,-0.5,500.5);
+      snprintf(name, 200, "EstITSTPCEstTracklet%s", aName);
+      fEst1Est2 = new TH2D(name,"ITS+TPC vs Tracklets",501,-0.5,5000.5,501,-0.5,500.5);
 
-  snprintf(name, 200, "EstITSTPCEstITSPure%s", aName);
-  fEst1Est3 = new TH2D(name,"ITS+TPC vs ITS Pure",501,-0.5,5000.5,801,-0.5,8000.5);
+      snprintf(name, 200, "EstITSTPCEstITSPure%s", aName);
+      fEst1Est3 = new TH2D(name,"ITS+TPC vs ITS Pure",501,-0.5,5000.5,801,-0.5,8000.5);
 
-  snprintf(name, 200, "EstTrackletEstITSPure%s", aName);
-  fEst2Est3 = new TH2D(name,"Tracklets vs ITS Pure",501,-0.5,5000.5,801,-0.5,8000.5);
+      snprintf(name, 200, "EstTrackletEstITSPure%s", aName);
+      fEst2Est3 = new TH2D(name,"Tracklets vs ITS Pure",501,-0.5,5000.5,801,-0.5,8000.5);
 
-  snprintf(name, 200, "EstITSTPCNormMult%s", aName);
-  fEst1Norm = new TH2D(name,"ITS+TPC vs Normalized Mult",501,-0.5,5000.5,501,-0.5,5000.5);
+      snprintf(name, 200, "EstITSTPCNormMult%s", aName);
+      fEst1Norm = new TH2D(name,"ITS+TPC vs Normalized Mult",501,-0.5,5000.5,501,-0.5,5000.5);
 
-  snprintf(name, 200, "EstTrackletsNormMult%s", aName);
-  fEst2Norm = new TH2D(name,"Tracklets vs Normalized Mult",501,-0.5,5000.5,501,-0.5,5000.5);
+      snprintf(name, 200, "EstTrackletsNormMult%s", aName);
+      fEst2Norm = new TH2D(name,"Tracklets vs Normalized Mult",501,-0.5,5000.5,501,-0.5,5000.5);
 
-  snprintf(name, 200, "EstITSPureNormMult%s", aName);
-  fEst3Norm = new TH2D(name,"ITS Pure vs Normalized Mult",501,-0.5,5000.5,501,-0.5,5000.5);
+      snprintf(name, 200, "EstITSPureNormMult%s", aName);
+      fEst3Norm = new TH2D(name,"ITS Pure vs Normalized Mult",501,-0.5,5000.5,501,-0.5,5000.5);
+    }
 
 }
 
@@ -104,7 +113,9 @@ AliFemtoCutMonitorEventMult::AliFemtoCutMonitorEventMult(const AliFemtoCutMonito
   fEst2Est3(0),
   fEst1Norm(0),
   fEst2Norm(0),
-  fEst3Norm(0)
+  fEst3Norm(0),
+  freadMC(kFALSE),
+  faddhists(kFALSE)
 {
   // copy constructor
   if (fEvMult) delete fEvMult;
@@ -113,38 +124,45 @@ AliFemtoCutMonitorEventMult::AliFemtoCutMonitorEventMult(const AliFemtoCutMonito
   if (fNormEvMult) delete fNormEvMult;
   fNormEvMult = new TH1D(*aCut.fNormEvMult);
 
-  if (fSPDMult) delete fSPDMult;
-  fSPDMult = new TH1D(*aCut.fSPDMult);
+  
+  if(!freadMC){
+    if (fSPDMult) delete fSPDMult;
+    fSPDMult = new TH1D(*aCut.fSPDMult);
+  }
 
   if (fMultSumPt) delete fMultSumPt;
   fMultSumPt = new TH2D(*aCut.fMultSumPt);
 
-  if (fEstimateITSTPC) delete fEstimateITSTPC;
-  fEstimateITSTPC = new TH1D(*aCut.fEstimateITSTPC);
 
-  if (fEstimateTracklets) delete fEstimateTracklets;
-  fEstimateTracklets = new TH1D(*aCut.fEstimateTracklets);
+  if(faddhists)
+    {
+      if (fEstimateITSTPC) delete fEstimateITSTPC;
+      fEstimateITSTPC = new TH1D(*aCut.fEstimateITSTPC);
+
+      if (fEstimateTracklets) delete fEstimateTracklets;
+      fEstimateTracklets = new TH1D(*aCut.fEstimateTracklets);
 
-  if (fEstimateITSPure) delete fEstimateITSPure;
-  fEstimateITSPure = new TH1D(*aCut.fEstimateITSPure);
+      if (fEstimateITSPure) delete fEstimateITSPure;
+      fEstimateITSPure = new TH1D(*aCut.fEstimateITSPure);
 
-  if (fEst1Est2) delete fEst1Est2;
-  fEst1Est2 = new TH2D(*aCut.fEst1Est2);
+      if (fEst1Est2) delete fEst1Est2;
+      fEst1Est2 = new TH2D(*aCut.fEst1Est2);
 
-  if (fEst1Est3) delete fEst1Est3;
-  fEst1Est3 = new TH2D(*aCut.fEst1Est3);
+      if (fEst1Est3) delete fEst1Est3;
+      fEst1Est3 = new TH2D(*aCut.fEst1Est3);
 
-  if (fEst2Est3) delete fEst2Est3;
-  fEst2Est3 = new TH2D(*aCut.fEst2Est3);
+      if (fEst2Est3) delete fEst2Est3;
+      fEst2Est3 = new TH2D(*aCut.fEst2Est3);
 
-  if (fEst1Norm) delete fEst1Norm;
-  fEst1Norm = new TH2D(*aCut.fEst1Norm);
+      if (fEst1Norm) delete fEst1Norm;
+      fEst1Norm = new TH2D(*aCut.fEst1Norm);
 
-  if (fEst2Norm) delete fEst2Norm;
-  fEst2Norm = new TH2D(*aCut.fEst2Norm);
+      if (fEst2Norm) delete fEst2Norm;
+      fEst2Norm = new TH2D(*aCut.fEst2Norm);
 
-  if (fEst3Norm) delete fEst3Norm;
-  fEst3Norm = new TH2D(*aCut.fEst3Norm);
+      if (fEst3Norm) delete fEst3Norm;
+      fEst3Norm = new TH2D(*aCut.fEst3Norm);
+    }
 }
 
 AliFemtoCutMonitorEventMult::~AliFemtoCutMonitorEventMult()
@@ -152,17 +170,23 @@ AliFemtoCutMonitorEventMult::~AliFemtoCutMonitorEventMult()
   // Destructor
   delete fEvMult;
   delete fNormEvMult;
-  delete fSPDMult;
+  if(!freadMC){
+    delete fSPDMult;
+  }
   delete fMultSumPt;
-  delete fEstimateITSTPC;
-  delete fEstimateTracklets;
-  delete fEstimateITSPure;
-  delete fEst1Est2;
-  delete fEst1Est3;
-  delete fEst2Est3;
-  delete fEst1Norm;
-  delete fEst2Norm;
-  delete fEst3Norm;      
+
+  if(faddhists)
+    {
+      delete fEstimateITSTPC;
+      delete fEstimateTracklets;
+      delete fEstimateITSPure;
+      delete fEst1Est2;
+      delete fEst1Est3;
+      delete fEst2Est3;
+      delete fEst1Norm;
+      delete fEst2Norm;
+      delete fEst3Norm;
+    }      
 }
 
 AliFemtoCutMonitorEventMult& AliFemtoCutMonitorEventMult::operator=(const AliFemtoCutMonitorEventMult& aCut)
@@ -177,38 +201,44 @@ AliFemtoCutMonitorEventMult& AliFemtoCutMonitorEventMult::operator=(const AliFem
   if (fNormEvMult) delete fNormEvMult;
   fNormEvMult = new TH1D(*aCut.fNormEvMult);
   
-  if (fSPDMult) delete fSPDMult;
-  fSPDMult = new TH1D(*aCut.fSPDMult);
+  if(!freadMC){
+    if (fSPDMult) delete fSPDMult;
+    fSPDMult = new TH1D(*aCut.fSPDMult);
+  }
   
   if (fMultSumPt) delete fMultSumPt;
   fMultSumPt = new TH2D(*aCut.fMultSumPt);
 
-  if (fEstimateITSTPC) delete fEstimateITSTPC;
-  fEstimateITSTPC = new TH1D(*aCut.fEstimateITSTPC);
 
-  if (fEstimateTracklets) delete fEstimateTracklets;
-  fEstimateTracklets = new TH1D(*aCut.fEstimateTracklets);
+  if(faddhists)
+    {
+      if (fEstimateITSTPC) delete fEstimateITSTPC;
+      fEstimateITSTPC = new TH1D(*aCut.fEstimateITSTPC);
 
-  if (fEstimateITSPure) delete fEstimateITSPure;
-  fEstimateITSPure = new TH1D(*aCut.fEstimateITSPure);
+      if (fEstimateTracklets) delete fEstimateTracklets;
+      fEstimateTracklets = new TH1D(*aCut.fEstimateTracklets);
 
-  if (fEst1Est2) delete fEst1Est2;
-  fEst1Est2 = new TH2D(*aCut.fEst1Est2);
+      if (fEstimateITSPure) delete fEstimateITSPure;
+      fEstimateITSPure = new TH1D(*aCut.fEstimateITSPure);
 
-  if (fEst1Est3) delete fEst1Est3;
-  fEst1Est3 = new TH2D(*aCut.fEst1Est3);
+      if (fEst1Est2) delete fEst1Est2;
+      fEst1Est2 = new TH2D(*aCut.fEst1Est2);
 
-  if (fEst2Est3) delete fEst2Est3;
-  fEst2Est3 = new TH2D(*aCut.fEst2Est3);
+      if (fEst1Est3) delete fEst1Est3;
+      fEst1Est3 = new TH2D(*aCut.fEst1Est3);
 
-  if (fEst1Norm) delete fEst1Norm;
-  fEst1Norm = new TH2D(*aCut.fEst1Norm);
+      if (fEst2Est3) delete fEst2Est3;
+      fEst2Est3 = new TH2D(*aCut.fEst2Est3);
 
-  if (fEst2Norm) delete fEst2Norm;
-  fEst2Norm = new TH2D(*aCut.fEst2Norm);
+      if (fEst1Norm) delete fEst1Norm;
+      fEst1Norm = new TH2D(*aCut.fEst1Norm);
 
-  if (fEst3Norm) delete fEst3Norm;
-  fEst3Norm = new TH2D(*aCut.fEst3Norm);
+      if (fEst2Norm) delete fEst2Norm;
+      fEst2Norm = new TH2D(*aCut.fEst2Norm);
+
+      if (fEst3Norm) delete fEst3Norm;
+      fEst3Norm = new TH2D(*aCut.fEst3Norm);
+    }
 
   return *this;
 }
@@ -225,18 +255,23 @@ void AliFemtoCutMonitorEventMult::Fill(const AliFemtoEvent* aEvent)
   // Fill in the monitor histograms with the values from the current track
   fEvMult->Fill(aEvent->NumberOfTracks());
   fNormEvMult->Fill(aEvent->UncorrectedNumberOfPrimaries());
-  fSPDMult->Fill(aEvent->SPDMultiplicity());
+  if(!freadMC){
+    fSPDMult->Fill(aEvent->SPDMultiplicity());
+  }
   fMultSumPt->Fill(aEvent->UncorrectedNumberOfPrimaries(), aEvent->ZDCEMEnergy());
 
-  fEstimateITSTPC->Fill(aEvent->MultiplicityEstimateITSTPC());
-  fEstimateTracklets->Fill(aEvent->MultiplicityEstimateTracklets());
-  fEstimateITSPure->Fill(aEvent->MultiplicityEstimateITSPure());
-  fEst1Est2->Fill(aEvent->MultiplicityEstimateITSTPC(),aEvent->MultiplicityEstimateTracklets());
-  fEst1Est3->Fill(aEvent->MultiplicityEstimateITSTPC(),aEvent->MultiplicityEstimateITSPure());
-  fEst2Est3->Fill(aEvent->MultiplicityEstimateTracklets(),aEvent->MultiplicityEstimateITSPure());
-  fEst1Norm->Fill(aEvent->MultiplicityEstimateITSTPC(),aEvent->UncorrectedNumberOfPrimaries());
-  fEst2Norm->Fill(aEvent->MultiplicityEstimateTracklets(),aEvent->UncorrectedNumberOfPrimaries());
-  fEst3Norm->Fill(aEvent->MultiplicityEstimateITSPure(),aEvent->UncorrectedNumberOfPrimaries());
+  if(faddhists)
+    {
+      fEstimateITSTPC->Fill(aEvent->MultiplicityEstimateITSTPC());
+      fEstimateTracklets->Fill(aEvent->MultiplicityEstimateTracklets());
+      fEstimateITSPure->Fill(aEvent->MultiplicityEstimateITSPure());
+      fEst1Est2->Fill(aEvent->MultiplicityEstimateITSTPC(),aEvent->MultiplicityEstimateTracklets());
+      fEst1Est3->Fill(aEvent->MultiplicityEstimateITSTPC(),aEvent->MultiplicityEstimateITSPure());
+      fEst2Est3->Fill(aEvent->MultiplicityEstimateTracklets(),aEvent->MultiplicityEstimateITSPure());
+      fEst1Norm->Fill(aEvent->MultiplicityEstimateITSTPC(),aEvent->UncorrectedNumberOfPrimaries());
+      fEst2Norm->Fill(aEvent->MultiplicityEstimateTracklets(),aEvent->UncorrectedNumberOfPrimaries());
+      fEst3Norm->Fill(aEvent->MultiplicityEstimateITSPure(),aEvent->UncorrectedNumberOfPrimaries());
+    }
 
 }
 
@@ -245,18 +280,23 @@ void AliFemtoCutMonitorEventMult::Write()
   // Write out the relevant histograms
   fEvMult->Write();
   fNormEvMult->Write();
-  fSPDMult->Write();
+  if(!freadMC){
+    fSPDMult->Write();
+  }
   fMultSumPt->Write();
 
-  fEstimateITSTPC->Write();
-  fEstimateTracklets->Write();
-  fEstimateITSPure->Write();
-  fEst1Est2->Write();
-  fEst1Est3->Write();
-  fEst2Est3->Write();
-  fEst1Norm->Write();
-  fEst2Norm->Write();
-  fEst3Norm->Write();
+  if(faddhists)
+    {
+      fEstimateITSTPC->Write();
+      fEstimateTracklets->Write();
+      fEstimateITSPure->Write();
+      fEst1Est2->Write();
+      fEst1Est3->Write();
+      fEst2Est3->Write();
+      fEst1Norm->Write();
+      fEst2Norm->Write();
+      fEst3Norm->Write();
+    }
 
 }
 
@@ -268,15 +308,29 @@ TList *AliFemtoCutMonitorEventMult::GetOutputList()
   tOutputList->Add(fSPDMult);
   tOutputList->Add(fMultSumPt);
 
-  tOutputList->Add(fEstimateITSTPC);
-  tOutputList->Add(fEstimateTracklets);
-  tOutputList->Add(fEstimateITSPure);
-  tOutputList->Add(fEst1Est2);
-  tOutputList->Add(fEst1Est3);
-  tOutputList->Add(fEst2Est3);
-  tOutputList->Add(fEst1Norm);
-  tOutputList->Add(fEst2Norm);
-  tOutputList->Add(fEst3Norm);
+  if(faddhists)
+    {
+      tOutputList->Add(fEstimateITSTPC);
+      tOutputList->Add(fEstimateTracklets);
+      tOutputList->Add(fEstimateITSPure);
+      tOutputList->Add(fEst1Est2);
+      tOutputList->Add(fEst1Est3);
+      tOutputList->Add(fEst2Est3);
+      tOutputList->Add(fEst1Norm);
+      tOutputList->Add(fEst2Norm);
+      tOutputList->Add(fEst3Norm);
+    }
   
   return tOutputList;
 }
+
+void AliFemtoCutMonitorEventMult::SetReadMC(Bool_t mc)
+{
+  freadMC=mc;
+}
+
+void AliFemtoCutMonitorEventMult::AdditionalMultHistsOn(Bool_t addhists)
+{
+  faddhists=addhists;
+}
+
index f5baf87..143109d 100644 (file)
@@ -39,6 +39,8 @@ class AliFemtoCutMonitorEventMult : public AliFemtoCutMonitor{
   virtual void Fill(const AliFemtoEvent* aEvent,const AliFemtoParticleCollection* aCollection)
   {AliFemtoCutMonitor::Fill(aEvent, aCollection);}
 
+  void SetReadMC(Bool_t mc);
+  void AdditionalMultHistsOn(Bool_t addhists);
   void Write();
 
   virtual TList *GetOutputList();
@@ -48,6 +50,9 @@ class AliFemtoCutMonitorEventMult : public AliFemtoCutMonitor{
   TH1D *fNormEvMult; // Normalized event multiplicity distribution
   TH1D *fSPDMult;    // SPD tracklet multiplicity
   TH2D *fMultSumPt;  // Event total pT vs. multiplicity
+
+  Bool_t freadMC;     // If true - add only one histogram to the output
+  Bool_t faddhists;   // If true - add only additional multiplicity histograms
   
   TH1D *fEstimateITSTPC;     // Multiplicity estimate ITS+TPC
   TH1D *fEstimateTracklets;  // Multiplicity estimate Tracklets
index 78f43b5..48fb6cd 100644 (file)
 #include "AliESDtrack.h"
 #include "AliESDVertex.h"
 
+#include "AliMultiplicity.h"
+#include "AliCentrality.h"
+#include "AliEventplane.h"
+
 #include "AliFmPhysicalHelixD.h"
 #include "AliFmThreeVectorF.h"
 
@@ -30,6 +34,8 @@
 
 #include "AliVertexerTracks.h"
 
+#include "AliPID.h"
+
 ClassImp(AliFemtoEventReaderESDChainKine)
 
 #if !(ST_NO_NAMESPACES)
@@ -48,7 +54,10 @@ AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine():
   fEvent(0x0),
   fStack(0x0),
   fGenHeader(0x0),
-  fRotateToEventPlane(0)
+  fRotateToEventPlane(0),
+  fESDpid(0),
+  fIsPidOwner(0)
+
 {
   //constructor with 0 parameters , look at default settings 
 }
@@ -65,7 +74,9 @@ AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoE
   fEvent(0x0),
   fStack(0x0),
   fGenHeader(0x0),
-  fRotateToEventPlane(0)
+  fRotateToEventPlane(0),
+  fESDpid(0),
+  fIsPidOwner(0)
 {
   // Copy constructor
   fConstrained = aReader.fConstrained;
@@ -76,6 +87,8 @@ AliFemtoEventReaderESDChainKine::AliFemtoEventReaderESDChainKine(const AliFemtoE
   fEvent = new AliESDEvent();
   fStack = aReader.fStack;
   fRotateToEventPlane = aReader.fRotateToEventPlane;
+  fESDpid = aReader.fESDpid;
+  fIsPidOwner = aReader.fIsPidOwner;
 }
 //__________________
 AliFemtoEventReaderESDChainKine::~AliFemtoEventReaderESDChainKine()
@@ -101,6 +114,9 @@ AliFemtoEventReaderESDChainKine& AliFemtoEventReaderESDChainKine::operator=(cons
   fStack = aReader.fStack;
   fGenHeader = aReader.fGenHeader;
   fRotateToEventPlane = aReader.fRotateToEventPlane;
+  fESDpid = aReader.fESDpid;
+  fIsPidOwner = aReader.fIsPidOwner;
+
   return *this;
 }
 //__________________
@@ -124,6 +140,31 @@ bool AliFemtoEventReaderESDChainKine::GetConstrained() const
   return fConstrained;
 }
 //__________________
+void AliFemtoEventReaderESDChainKine::SetReadTPCInner(const bool readinner)
+{
+  fReadInner=readinner;
+}
+
+bool AliFemtoEventReaderESDChainKine::GetReadTPCInner() const
+{
+  return fReadInner;
+}
+
+//__________________
+void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
+{
+  fUseTPCOnly=usetpconly;
+}
+
+bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
+{
+  return fUseTPCOnly;
+}
+void AliFemtoEventReaderESDChainKine::SetUseMultiplicity(EstEventMult aType)
+{
+  fEstEventMult = aType;
+}
+//__________________
 AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
 {
   // Get the event, read all the relevant information
@@ -149,8 +190,19 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
   hbtEvent->SetZDCEMEnergy(fEvent->GetZDCEMEnergy());
   hbtEvent->SetZDCParticipants(fEvent->GetZDCParticipants());
   hbtEvent->SetTriggerMask(fEvent->GetTriggerMask());
-  hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
-       
+  //hbtEvent->SetTriggerCluster(fEvent->GetTriggerCluster());
+
+  if ((fEvent->IsTriggerClassFired("CINT1WU-B-NOPF-ALL")) ||
+      (fEvent->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) ||
+      (fEvent->IsTriggerClassFired("CINT1-B-NOPF-ALLNOTRD")) ||
+      (fEvent->IsTriggerClassFired("CINT1-B-NOPF-FASTNOTRD")))
+    hbtEvent->SetTriggerCluster(1);
+  else if ((fEvent->IsTriggerClassFired("CSH1WU-B-NOPF-ALL")) ||
+          (fEvent->IsTriggerClassFired("CSH1-B-NOPF-ALLNOTRD")))
+    hbtEvent->SetTriggerCluster(2);
+  else 
+    hbtEvent->SetTriggerCluster(0);
+
   //Vertex
   double fV1[3];
   double fVCov[6];
@@ -161,42 +213,13 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
       fVCov[4] = -1001.0;
   }
   else {
-    if (fEvent->GetPrimaryVertex()) {
-      fEvent->GetPrimaryVertex()->GetXYZ(fV1);
-      fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
-
-      if (!fEvent->GetPrimaryVertex()->GetStatus()) {
-       // Get the vertex from SPD
-       fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
-       fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
-       
-       
-       if (!fEvent->GetPrimaryVertexSPD()->GetStatus())
-         fVCov[4] = -1001.0;
-       else {
-         fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
-         fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
-       }
-      }
-    }
-    else {
-      if (fEvent->GetPrimaryVertexSPD()) {
-       fEvent->GetPrimaryVertexSPD()->GetXYZ(fV1);
-       fEvent->GetPrimaryVertexSPD()->GetCovMatrix(fVCov);
-      }
-    }
-    if ((!fEvent->GetPrimaryVertex()) && (!fEvent->GetPrimaryVertexSPD()))
-      {
-       cout << "No vertex found !!!" << endl;
-       fV1[0] = 10000.0;
-       fV1[1] = 10000.0;
-       fV1[2] = 10000.0;
-       fVCov[4] = -1001.0;
-      }
+    fEvent->GetPrimaryVertex()->GetXYZ(fV1);
+    fEvent->GetPrimaryVertex()->GetCovMatrix(fVCov);
+    if (!fEvent->GetPrimaryVertex()->GetStatus())
+      fVCov[4] = -1001.0;
   }
 
   AliFmThreeVectorF vertex(fV1[0],fV1[1],fV1[2]);
-  
   hbtEvent->SetPrimVertPos(vertex);
   hbtEvent->SetPrimVertCov(fVCov);
 
@@ -222,11 +245,38 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
 
   hbtEvent->SetReactionPlaneAngle(tReactionPlane);
 
+  Int_t spdetaonecount = 0;
+  
+  for (int iter=0; iter<fEvent->GetMultiplicity()->GetNumberOfTracklets(); iter++) 
+    if (fabs(fEvent->GetMultiplicity()->GetEta(iter)) < 1.0)
+      spdetaonecount++;
+
+  //  hbtEvent->SetSPDMult(fEvent->GetMultiplicity()->GetNumberOfTracklets());
+  hbtEvent->SetSPDMult(spdetaonecount);
+
   //starting to reading tracks
   int nofTracks=0;  //number of reconstructed tracks in event
   nofTracks=fEvent->GetNumberOfTracks();
   int realnofTracks=0;//number of track which we use ina analysis
 
+  int tNormMult = 0;
+  int tNormMultPos = 0;
+  int tNormMultNeg = 0;
+
+  Float_t tTotalPt = 0.0;
+
+  Float_t b[2];
+  Float_t bCov[3];
+
+  Int_t tTracklet=0, tITSTPC=0, tITSPure=0;
+  
+  fEvent->EstimateMultiplicity(tTracklet, tITSTPC, tITSPure, 1.2);
+  
+  hbtEvent->SetMultiplicityEstimateITSTPC(tITSTPC);
+  hbtEvent->SetMultiplicityEstimateTracklets(tTracklet);
+  //  hbtEvent->SetMultiplicityEstimateITSPure(tITSPure);
+  hbtEvent->SetMultiplicityEstimateITSPure(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
+
   Int_t *motherids;
   if (fStack) {
     motherids = new Int_t[fStack->GetNtrack()];
@@ -264,10 +314,67 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
       //      cout << "Reading track " << i << endl;
       bool  tGoodMomentum=true; //flaga to chcek if we can read momentum of this track
                
-      AliFemtoTrack* trackCopy = new AliFemtoTrack();  
+       
       const AliESDtrack *esdtrack=fEvent->GetTrack(i);//getting next track
       //      const AliESDfriendTrack *tESDfriendTrack = esdtrack->GetFriendTrack();
 
+
+      if ((esdtrack->GetStatus() & AliESDtrack::kTPCrefit) &&
+         (esdtrack->GetStatus() & AliESDtrack::kITSrefit)) {
+       if (esdtrack->GetTPCNcls() > 70) 
+         if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
+           if (TMath::Abs(esdtrack->Eta()) < 1.2) {
+             esdtrack->GetImpactParameters(b,bCov);
+             if ((b[0]<0.2) && (b[1] < 0.25)) {
+               tNormMult++;
+               tTotalPt += esdtrack->Pt();
+             }
+           }
+         }
+      }
+      else if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) {
+       if (esdtrack->GetTPCNcls() > 100) 
+         if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 4.0) {
+           if (TMath::Abs(esdtrack->Eta()) < 1.2) {
+             esdtrack->GetImpactParameters(b,bCov);
+             if ((b[0]<2.4) && (b[1] < 3.2)) {
+               tNormMult++;
+               tTotalPt += esdtrack->Pt();
+             }
+           }
+         }
+      }
+      
+      hbtEvent->SetZDCEMEnergy(tTotalPt);
+      //       if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit)
+      //       if (esdtrack->GetTPCNcls() > 80) 
+      //         if (esdtrack->GetTPCchi2()/esdtrack->GetTPCNcls() < 6.0) 
+      //           if (esdtrack->GetConstrainedParam())
+      //             if (fabs(esdtrack->GetConstrainedParam()->Eta()) < 0.5)
+      //               if (esdtrack->GetConstrainedParam()->Pt() < 1.0) {
+      //                 if (esdtrack->GetSign() > 0)
+      //                   tNormMultPos++;
+      //                 else if (esdtrack->GetSign() < 0)
+      //                   tNormMultNeg--;
+      //               }
+
+      // If reading ITS-only tracks, reject all with TPC
+      if (fTrackType == kITSOnly) {
+       if (esdtrack->GetStatus() & AliESDtrack::kTPCrefit) continue;
+       if (!(esdtrack->GetStatus() & AliESDtrack::kITSrefit)) continue;
+       if (esdtrack->GetStatus() & AliESDtrack::kTPCin) continue;
+       UChar_t iclm = esdtrack->GetITSClusterMap();
+       Int_t incls = 0;
+       for (int iter=0; iter<6; iter++) if (iclm&(1<<iter)) incls++;
+       if (incls<=3) {
+         if (Debug()>1) cout << "Rejecting track with " << incls << " clusters" << endl;
+         continue;
+       }
+      }
+
+
+
+      AliFemtoTrack* trackCopy = new AliFemtoTrack();
       trackCopy->SetCharge((short)esdtrack->GetSign());
 
       //in aliroot we have AliPID 
@@ -280,6 +387,76 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
       trackCopy->SetPidProbPion(esdpid[2]);
       trackCopy->SetPidProbKaon(esdpid[3]);
       trackCopy->SetPidProbProton(esdpid[4]);
+
+      esdpid[0] = -100000.0;
+      esdpid[1] = -100000.0;
+      esdpid[2] = -100000.0;
+      esdpid[3] = -100000.0;
+      esdpid[4] = -100000.0;
+
+      double tTOF = 0.0;
+
+      if (esdtrack->GetStatus()&AliESDtrack::kTOFpid) {
+       tTOF = esdtrack->GetTOFsignal();
+       esdtrack->GetIntegratedTimes(esdpid);
+      }
+
+      trackCopy->SetTofExpectedTimes(tTOF-esdpid[2], tTOF-esdpid[3], tTOF-esdpid[4]);
+
+
+      //////  TPC ////////////////////////////////////////////
+      float nsigmaTPCK=-1000.;                                                        
+      float nsigmaTPCPi=-1000.;                                                        
+      float nsigmaTPCP=-1000.;  
+
+     if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTPCpid))){
+        nsigmaTPCK = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kKaon);
+        nsigmaTPCPi = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kPion);
+        nsigmaTPCP = fESDpid->NumberOfSigmasTPC(esdtrack,AliPID::kProton);
+
+      }
+      trackCopy->SetNSigmaTPCPi(nsigmaTPCPi);
+      trackCopy->SetNSigmaTPCK(nsigmaTPCK);
+      trackCopy->SetNSigmaTPCP(nsigmaTPCP);
+
+      ///// TOF ///////////////////////////////////////////////
+
+      float vp=-1000.;
+      float nsigmaTOFPi=-1000.;
+      float nsigmaTOFK=-1000.;
+      float nsigmaTOFP=-1000.;
+
+       if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
+           (esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
+           (esdtrack->GetStatus()&AliESDtrack::kTIME) &&
+           !(esdtrack->GetStatus()&AliESDtrack::kTOFmismatch))
+         {
+
+           //if ((esdtrack->GetStatus()&AliESDtrack::kTOFpid) &&
+           //(esdtrack->GetStatus()&AliESDtrack::kTOFout) &&
+           //(esdtrack->GetStatus()&AliESDtrack::kTIME)){
+           // collect info from ESDpid class
+         
+             if ((fESDpid) && (esdtrack->IsOn(AliESDtrack::kTOFpid))) {
+             
+             
+             double tZero = fESDpid->GetTOFResponse().GetStartTime(esdtrack->P());
+             
+             nsigmaTOFPi = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kPion,tZero);
+             nsigmaTOFK = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kKaon,tZero);
+             nsigmaTOFP = fESDpid->NumberOfSigmasTOF(esdtrack,AliPID::kProton,tZero);
+             
+             Double_t len=esdtrack->GetIntegratedLength();
+             Double_t tof=esdtrack->GetTOFsignal();
+             if(tof > 0.) vp=len/tof/0.03;
+           }
+         }
+
+      trackCopy->SetVTOF(vp);
+      trackCopy->SetNSigmaTOFPi(nsigmaTOFPi);
+      trackCopy->SetNSigmaTOFK(nsigmaTOFK);
+      trackCopy->SetNSigmaTOFP(nsigmaTOFP);
+
                                                
       double pxyz[3];
       double rxyz[3];
@@ -298,6 +475,17 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
        param->PropagateToDCA(fEvent->GetPrimaryVertexTPC(), (fEvent->GetMagneticField()), 10000, impact, covimpact);
        param->GetPxPyPz(pxyz);//reading noconstarined momentum
 
+         if (fReadInner == true) {
+           AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
+           tInfo->SetPDGPid(211);
+           tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
+           tInfo->SetMass(0.13957);
+           //    tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
+           //    tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
+           tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
+           trackCopy->SetHiddenInfo(tInfo);
+         }
+
        if (fRotateToEventPlane) {
          double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
          double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
@@ -312,7 +500,6 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
          delete trackCopy;
          continue;
        }
-
        trackCopy->SetP(v);//setting momentum
        trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
 
@@ -333,11 +520,47 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
        delete param;
       }
       else {
+         if (fReadInner == true) {
+         
+           if (esdtrack->GetTPCInnerParam()) {
+             AliExternalTrackParam *param = new AliExternalTrackParam(*esdtrack->GetTPCInnerParam());
+             param->GetXYZ(rxyz);
+             //            param->PropagateToDCA(fEvent->GetPrimaryVertex(), (fEvent->GetMagneticField()), 10000);
+             param->GetPxPyPz(pxyz);//reading noconstarined momentum
+             delete param;
+           
+             AliFemtoModelHiddenInfo *tInfo = new AliFemtoModelHiddenInfo();
+             tInfo->SetPDGPid(211);
+             tInfo->SetTrueMomentum(pxyz[0], pxyz[1], pxyz[2]);
+             tInfo->SetMass(0.13957);
+             //            tInfo->SetEmissionPoint(rxyz[0], rxyz[1], rxyz[2], 0.0);
+             //tInfo->SetEmissionPoint(fV1[0], fV1[1], fV1[2], 0.0);
+             tInfo->SetEmissionPoint(rxyz[0]-fV1[0], rxyz[1]-fV1[1], rxyz[2]-fV1[2], 0.0);
+             trackCopy->SetHiddenInfo(tInfo);
+           }
+         }
+
+
+       if(fTrackType == kGlobal) {
        if (fConstrained==true)             
          tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
        else
          tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
-       
+       }
+         else if (fTrackType == kTPCOnly) {
+           if (esdtrack->GetTPCInnerParam())
+             esdtrack->GetTPCInnerParam()->GetPxPyPz(pxyz);
+           else {
+             delete trackCopy;
+             continue;
+           }
+         }
+         else if (fTrackType == kITSOnly) {
+           if (fConstrained==true)                 
+             tGoodMomentum=esdtrack->GetConstrainedPxPyPz(pxyz); //reading constrained momentum
+           else
+             tGoodMomentum=esdtrack->GetPxPyPz(pxyz);//reading noconstarined momentum
+         }
 
        if (fRotateToEventPlane) {
          double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
@@ -412,6 +635,7 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
       }
       trackCopy->SetKinkIndexes(indexes);
 
+
       // Fill the hidden information with the simulated data
       if (TMath::Abs(esdtrack->GetLabel()) < fStack->GetNtrack()) {
        TParticle *tPart = fStack->Particle(TMath::Abs(esdtrack->GetLabel()));
@@ -533,6 +757,52 @@ AliFemtoEvent* AliFemtoEventReaderESDChainKine::ReturnHbtEvent()
   delete [] motherids;
   
   hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event  
+
+  AliCentrality *cent = fEvent->GetCentrality();
+  if (cent) {
+    hbtEvent->SetCentralityV0(cent->GetCentralityPercentile("V0M"));
+    //    hbtEvent->SetCentralityFMD(cent->GetCentralityPercentile("FMD"));
+    hbtEvent->SetCentralitySPD1(cent->GetCentralityPercentile("CL1"));
+    //    hbtEvent->SetCentralityTrk(cent->GetCentralityPercentile("TRK"));
+
+    if (Debug()>1) printf("  FemtoReader Got Event with %f %f %f %f\n", cent->GetCentralityPercentile("V0M"), 0.0, cent->GetCentralityPercentile("CL1"), 0.0);
+  }
+
+  if (fEstEventMult == kGlobalCount) 
+    hbtEvent->SetNormalizedMult(tNormMult);
+  else if (fEstEventMult == kTracklet)
+    hbtEvent->SetNormalizedMult(tTracklet);
+  else if (fEstEventMult == kITSTPC)
+    hbtEvent->SetNormalizedMult(tITSTPC);
+  else if (fEstEventMult == kITSPure)
+    hbtEvent->SetNormalizedMult(tITSPure);
+  else if (fEstEventMult == kSPDLayer1)
+    hbtEvent->SetNormalizedMult(fEvent->GetMultiplicity()->GetNumberOfITSClusters(1));
+  else if (fEstEventMult == kV0Centrality) {
+    // centrality between 0 (central) and 1 (very peripheral)
+
+    if (cent) {
+      if (cent->GetCentralityPercentile("V0M") < 0.00001)
+       hbtEvent->SetNormalizedMult(-1);
+      else
+       hbtEvent->SetNormalizedMult(lrint(10.0*cent->GetCentralityPercentile("V0M")));
+      if (Debug()>1) printf ("Set Centrality %i %f %li\n", hbtEvent->UncorrectedNumberOfPrimaries(), 
+                            10.0*cent->GetCentralityPercentile("V0M"), lrint(10.0*cent->GetCentralityPercentile("V0M")));
+    }
+  }
+
+  if (tNormMultPos > tNormMultNeg)
+    hbtEvent->SetZDCParticipants(tNormMultPos);
+  else
+    hbtEvent->SetZDCParticipants(tNormMultNeg);
+  
+  AliEventplane* ep = fEvent->GetEventplane();
+  if (ep) {
+    hbtEvent->SetEP(ep);
+    hbtEvent->SetReactionPlaneAngle(ep->GetEventplane("Q"));
+  }
+
+
   fCurEvent++; 
   //  cout<<"end of reading nt "<<nofTracks<<" real number "<<realnofTracks<<endl;
   return hbtEvent; 
@@ -558,21 +828,10 @@ void AliFemtoEventReaderESDChainKine::SetGenEventHeader(AliGenEventHeader *aGenH
   // You must provide the address where it can be found
   fGenHeader = aGenHeader;
 }
-
-//__________________
-void AliFemtoEventReaderESDChainKine::SetUseTPCOnly(const bool usetpconly)
-{
-  fUseTPCOnly=usetpconly;
-}
 void AliFemtoEventReaderESDChainKine::SetRotateToEventPlane(short dorotate)
 {
   fRotateToEventPlane=dorotate;
 }
-
-bool AliFemtoEventReaderESDChainKine::GetUseTPCOnly() const
-{
-  return fUseTPCOnly;
-}
 Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double *covar)
 {
   // Calculates the number of sigma to the vertex.
@@ -614,3 +873,9 @@ Float_t AliFemtoEventReaderESDChainKine::GetSigmaToVertex(double *impact, double
   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
   return d;
 }
+
+void AliFemtoEventReaderESDChainKine::SetReadTrackType(ReadTrackType aType)
+{
+  fTrackType = aType;
+}
+
index 4af1998..09f07bc 100644 (file)
 #include <list>
 #include <AliGenEventHeader.h>
 
+#include "AliESDpid.h"
+
 class AliFemtoEvent;
 
 class AliFemtoEventReaderESDChainKine : public AliFemtoEventReader 
 {
  public:
+  enum TrackType {kGlobal=0, kTPCOnly=1, kITSOnly=2, kSPDTracklet=3};
+  typedef enum TrackType ReadTrackType;
+
+  enum EventMult {kTracklet=0, kITSTPC=1, kITSPure=2, kGlobalCount=3, kSPDLayer1=4, kV0Centrality=5 };
+  typedef enum EventMult EstEventMult;
+
   AliFemtoEventReaderESDChainKine();
   AliFemtoEventReaderESDChainKine(const AliFemtoEventReaderESDChainKine& aReader);
   ~AliFemtoEventReaderESDChainKine();
@@ -38,17 +46,26 @@ class AliFemtoEventReaderESDChainKine : public AliFemtoEventReader
   bool GetConstrained() const;
   void SetUseTPCOnly(const bool usetpconly);
   bool GetUseTPCOnly() const;
+  void SetReadTPCInner(const bool readinner);
+  bool GetReadTPCInner() const;
+
+  void SetUseMultiplicity(EstEventMult aType);
+
+  void SetReadTrackType(ReadTrackType aType);
 
   void SetESDSource(AliESDEvent *aESD);
   void SetStackSource(AliStack *aStack);
   void SetGenEventHeader(AliGenEventHeader *aGenHeader);
   void SetRotateToEventPlane(short dorotate);
 
+  void SetESDPid(AliESDpid *esdPid) { fESDpid = esdPid; }
+
  protected:
 
  private:
   string         fFileName;      // name of current ESD file
   bool           fConstrained;   // flag to set which momentum from ESD file will be use
+  bool           fReadInner;     // flag to set if one wants to read TPC-only momentum
   bool           fUseTPCOnly;    // flog to set to read TPC only momentum instead of the full
   int            fNumberofEvent; // number of Events in ESD file
   int            fCurEvent;      // number of current event
@@ -56,11 +73,16 @@ class AliFemtoEventReaderESDChainKine : public AliFemtoEventReader
   AliESDEvent*   fEvent;         // ESD event
   AliStack      *fStack;         // Kinematics stack pointer
   AliGenEventHeader *fGenHeader; // Link to the generator event header
+  ReadTrackType  fTrackType;     // Type of track read
+  EstEventMult   fEstEventMult;  // Type of the event multiplicity estimator
 
   short          fRotateToEventPlane; // Rotate the event so that event plane is at x=0
 
   Float_t GetSigmaToVertex(double *impact, double *covar);
 
+  AliESDpid *fESDpid;
+  Bool_t fIsPidOwner;
+
 #ifdef __ROOT__
   ClassDef(AliFemtoEventReaderESDChainKine, 1)
 #endif
diff --git a/PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderKinematicsChain.cxx b/PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderKinematicsChain.cxx
new file mode 100644 (file)
index 0000000..cc90a13
--- /dev/null
@@ -0,0 +1,347 @@
+/////////////////////////////////////////////////////////////////////////////////////
+//                                                                                 //
+// AliFemtoEventReaderKinematicsChain - the reader class for the Alice ESD and     //
+// the model Kinematics information tailored for the Task framework and the        //
+// Reads in AliESDfriend to create shared hit/quality information                  //
+// Authors: Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch     //
+//          Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch //
+//                                                                                //
+/////////////////////////////////////////////////////////////////////////////////////
+
+#include "AliFemtoEventReaderKinematicsChain.h"
+
+#include "TFile.h"
+#include "TTree.h"
+#include "TList.h"
+
+#include "AliFmPhysicalHelixD.h"
+#include "AliFmThreeVectorF.h"
+
+#include "SystemOfUnits.h"
+
+#include "AliFemtoEvent.h"
+
+#include "TParticle.h"
+#include "AliStack.h"
+#include "TParticlePDG.h"
+#include "AliFemtoModelHiddenInfo.h"
+#include "AliFemtoModelGlobalHiddenInfo.h"
+#include "AliGenHijingEventHeader.h"
+#include "AliGenCocktailEventHeader.h"
+
+#include "AliVertexerTracks.h"
+
+ClassImp(AliFemtoEventReaderKinematicsChain)
+
+#if !(ST_NO_NAMESPACES)
+  using namespace units;
+#endif
+
+using namespace std;
+//____________________________
+AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain():
+  fFileName(" "),
+  fConstrained(true),
+  fNumberofEvent(0),
+  fCurEvent(0),
+  fCurFile(0),
+  fStack(0x0),
+  fGenHeader(0x0),
+  fRotateToEventPlane(0)
+{
+  //constructor with 0 parameters , look at default settings 
+}
+
+//__________________
+AliFemtoEventReaderKinematicsChain::AliFemtoEventReaderKinematicsChain(const AliFemtoEventReaderKinematicsChain& aReader):
+  AliFemtoEventReader(aReader),
+  fFileName(" "),
+  fConstrained(true),
+  fNumberofEvent(0),
+  fCurEvent(0),
+  fCurFile(0),
+  fStack(0x0),
+  fGenHeader(0x0),
+  fRotateToEventPlane(0)
+{
+  // Copy constructor
+  fConstrained = aReader.fConstrained;
+  fNumberofEvent = aReader.fNumberofEvent;
+  fCurEvent = aReader.fCurEvent;
+  fCurFile = aReader.fCurFile;
+  fStack = aReader.fStack;
+  fRotateToEventPlane = aReader.fRotateToEventPlane;
+}
+//__________________
+AliFemtoEventReaderKinematicsChain::~AliFemtoEventReaderKinematicsChain()
+{
+  //Destructor
+  //delete fEvent;
+}
+
+//__________________
+AliFemtoEventReaderKinematicsChain& AliFemtoEventReaderKinematicsChain::operator=(const AliFemtoEventReaderKinematicsChain& aReader)
+{
+  // Assignment operator
+  if (this == &aReader)
+    return *this;
+
+  fConstrained = aReader.fConstrained;
+  fNumberofEvent = aReader.fNumberofEvent;
+  fCurEvent = aReader.fCurEvent;
+  fCurFile = aReader.fCurFile;
+  fStack = aReader.fStack;
+  fGenHeader = aReader.fGenHeader;
+  fRotateToEventPlane = aReader.fRotateToEventPlane;
+  return *this;
+}
+//__________________
+// Simple report
+AliFemtoString AliFemtoEventReaderKinematicsChain::Report()
+{
+  AliFemtoString temp = "\n This is the AliFemtoEventReaderKinematicsChain\n";
+  return temp;
+}
+
+//__________________
+void AliFemtoEventReaderKinematicsChain::SetConstrained(const bool constrained)
+{
+  // Select whether to read constrained or not constrained momentum
+  fConstrained=constrained;
+}
+//__________________
+bool AliFemtoEventReaderKinematicsChain::GetConstrained() const
+{
+  // Check whether we read constrained or not constrained momentum
+  return fConstrained;
+}
+//__________________
+AliFemtoEvent* AliFemtoEventReaderKinematicsChain::ReturnHbtEvent()
+{
+  // Get the event, read all the relevant information from the stack
+  // and fill the AliFemtoEvent class
+  // Returns a valid AliFemtoEvent
+  AliFemtoEvent *hbtEvent = 0;
+  string tFriendFileName;
+
+  cout << "AliFemtoEventReaderKinematlaicsChain::Starting to read event: "<<fCurEvent<<endl;
+       
+  hbtEvent = new AliFemtoEvent;
+  //setting basic things
+  //  hbtEvent->SetEventNumber(fEvent->GetEventNumber());
+  hbtEvent->SetRunNumber(0); //No Run number in Kinematics!
+  hbtEvent->SetMagneticField(0*kilogauss);//to check if here is ok
+  hbtEvent->SetZDCN1Energy(0);
+  hbtEvent->SetZDCP1Energy(0);
+  hbtEvent->SetZDCN2Energy(0);
+  hbtEvent->SetZDCP2Energy(0);
+  hbtEvent->SetZDCEMEnergy(0);
+  hbtEvent->SetZDCParticipants(0);
+  hbtEvent->SetTriggerMask(0);
+  hbtEvent->SetTriggerCluster(0);
+       
+  //Vertex
+  double fV1[3];
+  double fVCov[6];
+
+
+  AliFmThreeVectorF vertex(0,0,0);
+    
+  
+  hbtEvent->SetPrimVertPos(vertex);
+  hbtEvent->SetPrimVertCov(fVCov);
+
+  Double_t tReactionPlane = 0;
+
+  AliGenHijingEventHeader *hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader); 
+  if (!hdh) {
+    AliGenCocktailEventHeader *cdh = dynamic_cast<AliGenCocktailEventHeader *> (fGenHeader);
+    if (cdh) {
+      TList *tGenHeaders = cdh->GetHeaders();
+      for (int ihead = 0; ihead<tGenHeaders->GetEntries(); ihead++) {
+       hdh = dynamic_cast<AliGenHijingEventHeader *> (fGenHeader);     
+       if (hdh) break;
+      }
+    }
+  }
+
+  if (hdh)
+    {
+      tReactionPlane = hdh->ReactionPlaneAngle();
+      cout << "Got reaction plane " << tReactionPlane << endl;
+    }
+
+  hbtEvent->SetReactionPlaneAngle(tReactionPlane);
+
+  //starting to reading tracks
+  int nofTracks=0;  //number of all tracks in MC event
+  nofTracks=fStack->GetNtrack(); 
+  int realnofTracks=0;//number of track which we use in analysis
+
+
+  int tNormMult = 0;
+  for (int i=0;i<nofTracks;i++)
+    {
+
+      //take only primaries
+      if(!fStack->IsPhysicalPrimary(i)) {continue;}
+                 
+      AliFemtoTrack* trackCopy = new AliFemtoTrack();  
+       
+         //getting next track
+      TParticle *kinetrack= fStack->Particle(i);
+
+      //setting multiplicity
+        realnofTracks++;//real number of tracks (only primary particles)
+
+      //setting normalized multiplicity
+      if (kinetrack->Eta() < 0.9)
+       if(kinetrack->GetPDG()->Charge()/3!=0)
+         tNormMult++;
+         
+         
+         //charge
+      trackCopy->SetCharge((short)(fStack->Particle(i)->GetPDG()->Charge()/3));
+
+
+      //in aliroot we have AliPID 
+      //0-electron 1-muon 2-pion 3-kaon 4-proton 5-photon 6-pi0 7-neutron 8-kaon0 9-eleCon   
+      //we use only 5 first
+      double kinepid[5];
+      for(int pid_iter=0;pid_iter<5;pid_iter++)
+         kinepid[pid_iter]=0;
+
+      int pdgcode = kinetrack->GetPdgCode();
+      //proton
+      if(pdgcode==2212 || pdgcode==-2212)
+        kinepid[4]=1000;
+      //kaon
+      if(pdgcode==321 || pdgcode==-321 )
+        kinepid[3]=1000;
+      //pion
+      if( pdgcode==211 || pdgcode==-211)
+        kinepid[2]=1000;
+      //electron
+      if(pdgcode==11 || pdgcode==-11)
+        kinepid[0]=1000;
+      //muon
+      if(pdgcode==13 || pdgcode==-13)
+        kinepid[1]=1000;
+
+      trackCopy->SetPidProbElectron(kinepid[0]);
+      trackCopy->SetPidProbMuon(kinepid[1]);
+      trackCopy->SetPidProbPion(kinepid[2]);
+      trackCopy->SetPidProbKaon(kinepid[3]);
+      trackCopy->SetPidProbProton(kinepid[4]);
+                                       
+                                       
+       //Momentum
+      double pxyz[3];
+      double rxyz[3];
+     
+       pxyz[0]=kinetrack->Px();
+       pxyz[1]=kinetrack->Py();
+       pxyz[2]=kinetrack->Pz();
+
+       rxyz[0]=kinetrack->Vx();
+       rxyz[1]=kinetrack->Vy();
+       rxyz[2]=kinetrack->Vz();
+
+       if (fRotateToEventPlane) {
+         double tPhi = TMath::ATan2(pxyz[1], pxyz[0]);
+         double tRad = TMath::Hypot(pxyz[0], pxyz[1]);
+         
+         pxyz[0] = tRad*TMath::Cos(tPhi - tReactionPlane);
+         pxyz[1] = tRad*TMath::Sin(tPhi - tReactionPlane);
+       }
+
+       AliFemtoThreeVector v(pxyz[0],pxyz[1],pxyz[2]);
+       if (v.Mag() < 0.0001) {
+         //      cout << "Found 0 momentum ???? "  << pxyz[0] << " " << pxyz[1] << " " << pxyz[2] << endl;
+         delete trackCopy;
+         continue;
+       }
+
+       trackCopy->SetP(v);//setting momentum
+       trackCopy->SetPt(sqrt(pxyz[0]*pxyz[0]+pxyz[1]*pxyz[1]));
+       const AliFmThreeVectorD kP(pxyz[0],pxyz[1],pxyz[2]);
+       const AliFmThreeVectorD kOrigin(fV1[0],fV1[1],fV1[2]);
+
+       //label
+    trackCopy->SetLabel(i);
+
+
+       hbtEvent->TrackCollection()->push_back(trackCopy);//adding track to analysis
+       
+               
+    }
+  
+  hbtEvent->SetNumberOfTracks(realnofTracks);//setting number of track which we read in event
+  hbtEvent->SetNormalizedMult(tNormMult);
+  fCurEvent++; 
+
+  return hbtEvent; 
+}
+
+//___________________
+void AliFemtoEventReaderKinematicsChain::SetStackSource(AliStack *aStack)
+{
+  // The chain loads the stack for us
+  // You must provide the address where it can be found
+  fStack = aStack;
+}
+//___________________
+void AliFemtoEventReaderKinematicsChain::SetGenEventHeader(AliGenEventHeader *aGenHeader)
+{
+  // The chain loads the generator event header for us
+  // You must provide the address where it can be found
+  fGenHeader = aGenHeader;
+}
+
+//__________________
+void AliFemtoEventReaderKinematicsChain::SetRotateToEventPlane(short dorotate)
+{
+  fRotateToEventPlane=dorotate;
+}
+
+Float_t AliFemtoEventReaderKinematicsChain::GetSigmaToVertex(double *impact, double *covar)
+{
+  // Calculates the number of sigma to the vertex.
+
+  Float_t b[2];
+  Float_t bRes[2];
+  Float_t bCov[3];
+
+  b[0] = impact[0];
+  b[1] = impact[1];
+  bCov[0] = covar[0];
+  bCov[1] = covar[1];
+  bCov[2] = covar[2];
+
+  bRes[0] = TMath::Sqrt(bCov[0]);
+  bRes[1] = TMath::Sqrt(bCov[2]);
+
+  // -----------------------------------
+  // How to get to a n-sigma cut?
+  //
+  // The accumulated statistics from 0 to d is
+  //
+  // ->  Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
+  // ->  1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
+  //
+  // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
+  // Can this be expressed in a different way?
+
+  if (bRes[0] == 0 || bRes[1] ==0)
+    return -1;
+
+  Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
+
+  // stupid rounding problem screws up everything:
+  // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
+  if (TMath::Exp(-d * d / 2) < 1e-10)
+    return 1000;
+
+  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
+  return d;
+}
diff --git a/PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderKinematicsChain.h b/PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderKinematicsChain.h
new file mode 100644 (file)
index 0000000..2b61045
--- /dev/null
@@ -0,0 +1,69 @@
+/////////////////////////////////////////////////////////////////////////////////////
+//                                                                                 //
+// AliFemtoEventReaderKinematicsChain - the reader class for the Alice ESD and     //
+// the model Kinematics information tailored for the Task framework and the        //
+// Reads in AliESDfriend to create shared hit/quality information                  //
+// Authors: Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch     //
+//          Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch //
+//                                                                                //
+/////////////////////////////////////////////////////////////////////////////////////
+
+
+#ifndef ALIFEMTOEVENTREADERKINEMATICSCHAIN_H
+#define ALIFEMTOEVENTREADERKINEMATICSCHAIN_H
+
+#include "AliFemtoEventReader.h"
+#include "AliFemtoEnumeration.h"
+
+#include <string>
+#include <vector>
+#include <TTree.h>
+#include <AliStack.h>
+#include <list>
+#include <AliGenEventHeader.h>
+
+class AliFemtoEvent;
+
+class AliFemtoEventReaderKinematicsChain : public AliFemtoEventReader 
+{
+ public:
+  AliFemtoEventReaderKinematicsChain();
+  AliFemtoEventReaderKinematicsChain(const AliFemtoEventReaderKinematicsChain& aReader);
+  ~AliFemtoEventReaderKinematicsChain();
+
+  AliFemtoEventReaderKinematicsChain& operator=(const AliFemtoEventReaderKinematicsChain& aReader);
+
+  AliFemtoEvent* ReturnHbtEvent();
+  AliFemtoString Report();
+  void SetConstrained(const bool constrained);
+  bool GetConstrained() const;
+
+  //void SetESDSource(AliESDEvent *aESD);
+  void SetStackSource(AliStack *aStack);
+  void SetGenEventHeader(AliGenEventHeader *aGenHeader);
+  void SetRotateToEventPlane(short dorotate);
+
+ protected:
+
+ private:
+  string         fFileName;      // name of current ESD file
+  bool           fConstrained;   // flag to set which momentum from ESD file will be use
+  int            fNumberofEvent; // number of Events in ESD file
+  int            fCurEvent;      // number of current event
+  unsigned int   fCurFile;       // number of current file
+  AliStack       *fStack;         // Kinematics stack pointer
+  AliGenEventHeader *fGenHeader; // Link to the generator event header
+
+  short          fRotateToEventPlane; // Rotate the event so that event plane is at x=0
+
+  Float_t GetSigmaToVertex(double *impact, double *covar);
+
+#ifdef __ROOT__
+  ClassDef(AliFemtoEventReaderKinematicsChain, 1)
+#endif
+
+    };
+  
+#endif
+
+
diff --git a/PWG2/FEMTOSCOPY/AliFemto/AliFemtoMCTrackCut.cxx b/PWG2/FEMTOSCOPY/AliFemto/AliFemtoMCTrackCut.cxx
new file mode 100644 (file)
index 0000000..21c51b0
--- /dev/null
@@ -0,0 +1,185 @@
+/////////////////////////////////////////////////////////////////////////////////////
+//                                                                                 //
+// AliFemtoMCTrackCut: A basic track cut that used information from                //
+// ALICE MC to accept or reject the track.                                         //  
+// Enables the selection on charge, transverse momentum, rapidity,                 //
+// and PDG of the particle                                                        //
+// Authors: Malgorzata Janik, Warsaw University of Technology, majanik@cern.ch     //
+//          Lukasz Graczykowski, Warsaw University of Technology, lgraczyk@cern.ch //
+//                                                                                //
+/////////////////////////////////////////////////////////////////////////////////////
+
+
+
+#include "AliFemtoMCTrackCut.h"
+#include <cstdio>
+
+#ifdef __ROOT__ 
+ClassImp(AliFemtoMCTrackCut)
+#endif
+
+
+AliFemtoMCTrackCut::AliFemtoMCTrackCut() :
+    fCharge(0),
+    fLabel(0),
+    fNTracksPassed(0),
+    fNTracksFailed(0)
+{
+  // Default constructor
+  fNTracksPassed = fNTracksFailed = 0;
+  fCharge = 0;  // takes both charges 0
+  fPt[0]=0.0;              fPt[1] = 100.0;//100
+  fPDGcode = 0;
+  fRapidity[0]=-2;       fRapidity[1]=2;//-2 2
+  fEta[0]=-2;       fEta[1]=2;//-2 2
+  fLabel=false;
+}
+//------------------------------
+AliFemtoMCTrackCut::~AliFemtoMCTrackCut(){
+  /* noop */
+}
+//------------------------------
+bool AliFemtoMCTrackCut::Pass(const AliFemtoTrack* track)
+{
+
+  if (fLabel)
+    {
+      if(track->Label()<0)
+       {
+         fNTracksFailed++;
+         return false;
+       }    
+    }
+
+   if (fCharge!=0)
+    {               
+    if (fCharge==10)   
+       {
+       if(track->Charge()==0){
+         fNTracksFailed++;
+         return false;
+         }
+       } 
+      else if (track->Charge()!= fCharge)      
+       {
+         fNTracksFailed++;
+         return false;
+       }
+    }
+
+   if (fPDGcode!=0)
+     {
+     
+     if(fPDGcode==11 || fPDGcode==-11 )
+     { if(!fMass) fMass=0.000511;
+        if (track->PidProbElectron()!=1000)
+          {
+            fNTracksFailed++;
+            return false;
+          }
+          }
+       if(fPDGcode==13 || fPDGcode==-13)
+       {
+        if (track->PidProbMuon()!=1000)
+          {if(!fMass) fMass=0.105658;
+            fNTracksFailed++;
+            return false;
+          }
+          }
+       if(fPDGcode==211 || fPDGcode==-211 )
+       {
+        if (track->PidProbPion()!=1000)
+          {if(!fMass) fMass= 0.1395699;
+            fNTracksFailed++;
+            return false;
+          }
+          }
+       if(fPDGcode==2212 || fPDGcode==-2212 )
+       { if(!fMass) fMass=0.938272013;
+        if (track->PidProbProton()!=1000)
+          {
+            fNTracksFailed++;
+            return false;
+          }
+          }
+      if(fPDGcode==321 || fPDGcode==-321 )
+      { if(!fMass) fMass=0.493677;
+        if (track->PidProbKaon()!=1000)
+          {
+            fNTracksFailed++;
+            return false;
+          }
+          }
+     }
+
+  float tEnergy = ::sqrt(track->P().Mag2()+fMass*fMass);
+  //cout<<"MCTrackCut: tEnergy: "<<tEnergy<<endl;
+  //cout<<"MCTrackCut: track->P().z(): "<<track->P().z()<<endl;
+  //cout<<"MCTrackCut: tEnergy-track->P().z(): "<<tEnergy-track->P().z()<<endl;
+  float tRapidity;
+  if(tEnergy-track->P().z() == 0 || (tEnergy+track->P().z())/(tEnergy-track->P().z()) == 0)
+    {
+    fNTracksFailed++;
+    return false;
+    }
+  else
+    tRapidity = 0.5*::log((tEnergy+track->P().z())/(tEnergy-track->P().z()));
+  float tPt = ::sqrt((track->P().x())*(track->P().x())+(track->P().y())*(track->P().y()));
+  float tEta = track->P().PseudoRapidity();
+
+  if ((tRapidity<fRapidity[0])||(tRapidity>fRapidity[1]))
+    {
+      fNTracksFailed++;
+      return false;
+    }
+  if ((tEta<fEta[0])||(tEta>fEta[1]))
+    {
+      fNTracksFailed++;
+      return false;
+    }
+  if ((tPt<fPt[0])||(tPt>fPt[1]))
+    {
+      fNTracksFailed++;
+      return false;
+    }
+
+  fNTracksPassed++ ;
+  return true;
+    
+    
+}
+//------------------------------
+AliFemtoString AliFemtoMCTrackCut::Report()
+{
+  // Prepare report from the execution
+  string tStemp;
+  char tCtemp[100];
+  sprintf(tCtemp,"Particle mass:\t%E\n",this->Mass());
+  tStemp=tCtemp;
+  sprintf(tCtemp,"Particle charge:\t%d\n",fCharge);
+  tStemp+=tCtemp;
+  sprintf(tCtemp,"Particle pT:\t%E - %E\n",fPt[0],fPt[1]);
+  tStemp+=tCtemp;
+  sprintf(tCtemp,"Particle rapidity:\t%E - %E\n",fRapidity[0],fRapidity[1]);
+  tStemp+=tCtemp; 
+  sprintf(tCtemp,"Particle eta:\t%E - %E\n",fEta[0],fEta[1]);
+  tStemp+=tCtemp;
+ sprintf(tCtemp,"Number of tracks which passed:\t%ld  Number which failed:\t%ld\n",fNTracksPassed,fNTracksFailed);
+  tStemp += tCtemp;
+  AliFemtoString returnThis = tStemp;
+  return returnThis;
+}
+TList *AliFemtoMCTrackCut::ListSettings()
+{
+  // return a list of settings in a writable form
+  TList *tListSetttings = new TList();
+  char buf[200];
+  snprintf(buf, 200, "AliFemtoMCTrackCut.mass=%f", this->Mass());
+  tListSetttings->AddLast(new TObjString(buf));
+
+  snprintf(buf, 200, "AliFemtoMCTrackCut.charge=%i", fCharge);
+  tListSetttings->AddLast(new TObjString(buf));
+  return tListSetttings;
+}
+
+                           
diff --git a/PWG2/FEMTOSCOPY/AliFemto/AliFemtoMCTrackCut.h b/PWG2/FEMTOSCOPY/AliFemto/AliFemtoMCTrackCut.h
new file mode 100644 (file)
index 0000000..61c54d0
--- /dev/null
@@ -0,0 +1,68 @@
+///////////////////////////////////////////////////////////////////////////
+//                                                                       //
+// AliFemtoMCTrackCut: A basic track cut that used information from     //
+// ALICE MC to accept or reject the track.                              //  
+// Enables the selection on charge, transverse momentum, rapidity,       //
+// and PDG of the particle                                                                                      //
+// Authors: Malgorzata Janik (WUT)    majanik@cern.ch                                  //
+//                     Lukasz Graczykowski (WUT) lgraczyk@cern.ch                                       //
+//                                                                       //
+///////////////////////////////////////////////////////////////////////////
+
+#ifndef ALIFEMTOMCTRACKCUT_H
+#define ALIFEMTOMCTRACKCUT_H
+
+
+#include "AliFemtoTrackCut.h"
+
+class AliFemtoMCTrackCut : public AliFemtoTrackCut 
+{
+
+ public:
+  AliFemtoMCTrackCut();
+  virtual ~AliFemtoMCTrackCut();
+
+  virtual bool Pass(const AliFemtoTrack* aTrack);
+
+  virtual AliFemtoString Report();
+  virtual TList *ListSettings();
+  virtual AliFemtoParticleType Type(){return hbtTrack;}
+
+  void SetPt(const float& lo, const float& hi);
+  void SetRapidity(const float& lo, const float& hi);
+  void SetEta(const float& lo, const float& hi);
+  void SetCharge(const int& ch);
+  void SetPDG(const int& pdg);
+  void SetLabel(const bool& flag);
+
+
+ private:   // here are the quantities I want to cut on...
+
+  int               fCharge;             // particle charge
+  float             fPt[2];              // bounds for transverse momentum
+  float             fRapidity[2];        // bounds for rapidity
+  float             fEta[2];             // bounds for pseudorapidity
+  bool              fLabel;              // if true label<0 will not pass throught 
+  int               fPDGcode;            // PDG code of the particle
+  long              fNTracksPassed;      // passed tracks count
+  long              fNTracksFailed;      // failed tracks count
+
+
+
+#ifdef __ROOT__ 
+  ClassDef(AliFemtoMCTrackCut, 1)
+#endif
+    };
+
+
+inline void AliFemtoMCTrackCut::SetPt(const float& lo, const float& hi){fPt[0]=lo; fPt[1]=hi;}
+inline void AliFemtoMCTrackCut::SetRapidity(const float& lo,const float& hi){fRapidity[0]=lo; fRapidity[1]=hi;}
+inline void AliFemtoMCTrackCut::SetEta(const float& lo,const float& hi){fEta[0]=lo; fEta[1]=hi;}
+inline void AliFemtoMCTrackCut::SetCharge(const int& ch){fCharge = ch;}
+inline void AliFemtoMCTrackCut::SetPDG(const int& pdg){fPDGcode = pdg;}
+inline void AliFemtoMCTrackCut::SetLabel(const bool& flag){fLabel=flag;}
+
+
+#endif
+
index a18156c..2a9abb2 100644 (file)
@@ -34,17 +34,20 @@ ClassImp(AliFemtoPairCutAntiGamma)
 AliFemtoPairCutAntiGamma::AliFemtoPairCutAntiGamma():
   AliFemtoShareQualityPairCut(),
   fMaxEEMinv(0.0),
-  fMaxDTheta(0.0)
+  fMaxDTheta(0.0),
+  fDTPCMin(0)
 {
 }
 //__________________
 AliFemtoPairCutAntiGamma::AliFemtoPairCutAntiGamma(const AliFemtoPairCutAntiGamma& c) : 
   AliFemtoShareQualityPairCut(c),
   fMaxEEMinv(0.0),
-  fMaxDTheta(0.0)
+  fMaxDTheta(0.0),
+  fDTPCMin(0)
 { 
   fMaxEEMinv = c.fMaxEEMinv;
   fMaxDTheta = c.fMaxDTheta;
+  fDTPCMin = c.fDTPCMin;
 }
 
 //__________________
@@ -75,15 +78,28 @@ bool AliFemtoPairCutAntiGamma::Pass(const AliFemtoPair* pair){
     if ((minv < fMaxEEMinv) && (dtheta < fMaxDTheta)) temp = false;
   }
 
-  if (temp) {
+  bool tempTPCEntrance = true;
+  
+  double distx = pair->Track1()->Track()->NominalTpcEntrancePoint().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
+  double disty = pair->Track1()->Track()->NominalTpcEntrancePoint().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
+  double distz = pair->Track1()->Track()->NominalTpcEntrancePoint().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
+  double dist = sqrt(distx*distx + disty*disty + distz*distz);
+
+  tempTPCEntrance = dist > fDTPCMin;
+
+
+  if (temp && tempTPCEntrance) {
     temp = AliFemtoShareQualityPairCut::Pass(pair);
     if (temp) fNPairsPassed++;
     else fNPairsFailed++;
+    return temp;
   }
   else
+    {
     fNPairsFailed++;
+    return false;
+    }
 
-  return temp;
 }
 //__________________
 AliFemtoString AliFemtoPairCutAntiGamma::Report(){
@@ -119,3 +135,8 @@ void AliFemtoPairCutAntiGamma::SetMaxThetaDiff(Double_t maxdtheta)
 {
   fMaxDTheta = maxdtheta;
 }
+
+void AliFemtoPairCutAntiGamma::SetTPCEntranceSepMinimum(double dtpc)
+{
+  fDTPCMin = dtpc;
+}
index 547e916..064cb79 100644 (file)
@@ -40,10 +40,13 @@ public:
   virtual AliFemtoPairCut* Clone();
   void SetMaxEEMinv(Double_t maxeeminv);
   void SetMaxThetaDiff(Double_t maxdtheta);
+  void SetTPCEntranceSepMinimum(double dtpc);
   
  protected:
   Double_t fMaxEEMinv; // Maximum allowed ee Minv
   Double_t fMaxDTheta; // Maximum polar angle difference
+  Double_t fDTPCMin;          // Minimum allowed pair nominal separation at the entrance to the TPC
+
 
 #ifdef __ROOT__
   ClassDef(AliFemtoPairCutAntiGamma, 0)
index 621f371..49306d5 100644 (file)
@@ -62,5 +62,7 @@
 #pragma link C++ class AliFemtoAODTrackCut+;
 #pragma link C++ class AliAnalysisTaskFemto+;
 #pragma link C++ class AliTwoTrackRes+;
+#pragma link C++ class AliFemtoMCTrackCut+;
+#pragma link C++ class AliFemtoEventReaderKinematicsChain+;
 
 #endif
index 14f500e..72f0949 100644 (file)
@@ -41,3 +41,4 @@
 #pragma link C++ class AliFemtoQinvCorrFctnEMCIC;
 #pragma link C++ class AliFemtoCorrFctn3DSphericalEMCIC;
 #pragma link C++ class AliFemtoBPLCMS3DCorrFctnEMCIC;
+#pragma link C++ class AliFemtoPairCutPt;