Commit of Hypertriton task for Ramona

diff --git a/PWGLF/CMakelibPWGLFSTRANGENESS.pkg b/PWGLF/CMakelibPWGLFSTRANGENESS.pkg
index 06bdfb73d8636c861eab7395dff06477c46b26af..3a19d3a5422b57a6b1d7cc6fdc70e6c6caaf3cf2 100644 (file)
--- a/PWGLF/CMakelibPWGLFSTRANGENESS.pkg
+++ b/PWGLF/CMakelibPWGLFSTRANGENESS.pkg
@@ -52,6 +52,8 @@ set ( SRCS
   STRANGENESS/Cascades/AliAnalysisTaskCheckPerformanceCascadepp276.cxx
   STRANGENESS/Hypernuclei/AliAnalysisTaskLambdaNAOD.cxx
   STRANGENESS/Hypernuclei/AliAnalysisTaskHdibaryonLPpi.cxx
+  STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3Pi.cxx
+  STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3PiMC.cxx
 )
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )

diff --git a/PWGLF/PWGLFSTRANGENESSLinkDef.h b/PWGLF/PWGLFSTRANGENESSLinkDef.h
index 9cae9e98c157efec2f266f265ccbd8e9a8498d5b..2467e5ad59776bf57fc18dd562c22821524c0307 100644 (file)
--- a/PWGLF/PWGLFSTRANGENESSLinkDef.h
+++ b/PWGLF/PWGLFSTRANGENESSLinkDef.h
@@ -31,5 +31,7 @@
 #pragma link C++ class AliAnalysisTaskCheckPerformanceCascadepp276+;
 #pragma link C++ class AliAnalysisTaskLambdaNAOD+;
 #pragma link C++ class AliAnalysisTaskHdibaryonLPpi+;
+#pragma link C++ class AliAnalysisTaskHelium3Pi+;
+#pragma link C++ class AliAnalysisTaskHelium3PiMC+;
 
 #endif

diff --git a/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3Pi.cxx b/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3Pi.cxx
new file mode 100644 (file)
index 0000000..98bb836
--- /dev/null
+++ b/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3Pi.cxx
@@ -0,0 +1,971 @@
+/**************************************************************************\r
+ * Contributors are not mentioned at all.                                 *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+//\r
+//-----------------------------------------------------------------\r
+//                 AliAnalysisTaskHelium3Pi class\r
+//-----------------------------------------------------------------\r
+\r
+\r
+class TTree;\r
+class TParticle;\r
+class TVector3;\r
+\r
+#include "AliAnalysisManager.h"\r
+#include <AliMCEventHandler.h>\r
+#include <AliMCEvent.h>\r
+#include <AliStack.h>\r
+\r
+class AliESDVertex;\r
+class AliAODVertex;\r
+class AliESDv0;\r
+class AliAODv0; \r
+class AliCascadeVertexer;\r
+\r
+#include <iostream>\r
+#include "AliAnalysisTaskSE.h"\r
+#include "TList.h"\r
+#include "TH1.h"\r
+#include "TH2.h"\r
+#include "TH3.h"\r
+#include "TNtuple.h"\r
+#include "TGraph.h"\r
+#include "TCutG.h"\r
+#include "TF1.h"\r
+#include "TCanvas.h"\r
+#include "TMath.h"\r
+#include "TChain.h"\r
+#include "Riostream.h"\r
+#include "AliLog.h"\r
+#include "AliCascadeVertexer.h"\r
+#include "AliESDEvent.h"\r
+#include "AliESDtrack.h"\r
+#include "AliExternalTrackParam.h"\r
+#include "AliAODEvent.h"\r
+#include "AliInputEventHandler.h"\r
+#include "AliESDcascade.h"\r
+#include "AliAODcascade.h"\r
+#include "AliAnalysisTaskHelium3Pi.h"\r
+#include "AliESDtrackCuts.h"\r
+#include "AliCentrality.h"\r
+#include "TString.h"\r
+#include <TDatime.h>\r
+#include <TRandom3.h>\r
+#include <TLorentzVector.h>\r
+\r
+const Int_t AliAnalysisTaskHelium3Pi::fgNrot = 15;\r
+\r
+\r
+ClassImp(AliAnalysisTaskHelium3Pi)\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskHelium3Pi::AliAnalysisTaskHelium3Pi() \r
+: AliAnalysisTaskSE(),\r
+  fAnalysisType("ESD"), \r
+  fCollidingSystems(0), \r
+  fDataType("REAL"),\r
+  fListHistCascade(0), \r
+  fHistEventMultiplicity(0),         \r
+  fHistTrackMultiplicity(0),      \r
+  fHistTrackMultiplicityCent(0),      \r
+  fHistTrackMultiplicitySemiCent(0),  \r
+  fHistTrackMultiplicityMB(0),        \r
+  fHistTrackMultiplicityPVCent(0),      \r
+  fHistTrackMultiplicityPVSemiCent(0),  \r
+  fHistTrackMultiplicityPVMB(0),        \r
+  fHistMult(0),\r
+  fhBB(0),    \r
+  fhTOF(0),   \r
+  fhMassTOF(0),\r
+  fhBBPions(0),\r
+  fhBBHe(0),   \r
+  fhNaPos(0),  \r
+  fhNaNeg(0),  \r
+  fBetavsTPCsignalPos(0),  \r
+  fBetavsTPCsignalNeg(0),  \r
+  fHelium3TOF(0),                       \r
+  fhTreeESD(0)  \r
+\r
+{\r
+  // Dummy Constructor \r
+}\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskHelium3Pi::AliAnalysisTaskHelium3Pi(const char *name) \r
+  : AliAnalysisTaskSE(name), \r
+    fAnalysisType("ESD"), \r
+    fCollidingSystems(0), \r
+    fDataType("REAL"),\r
+    fListHistCascade(0), \r
+    fHistEventMultiplicity(0),    \r
+    fHistTrackMultiplicity(0),            \r
+    fHistTrackMultiplicityCent(0),      \r
+    fHistTrackMultiplicitySemiCent(0),  \r
+    fHistTrackMultiplicityMB(0),        \r
+    fHistTrackMultiplicityPVCent(0),      \r
+    fHistTrackMultiplicityPVSemiCent(0),  \r
+    fHistTrackMultiplicityPVMB(0),        \r
+    fHistMult(0),\r
+    fhBB(0),    \r
+    fhTOF(0),   \r
+    fhMassTOF(0),\r
+    fhBBPions(0),\r
+    fhBBHe(0),   \r
+    fhNaPos(0),  \r
+    fhNaNeg(0),  \r
+    fBetavsTPCsignalPos(0),  \r
+    fBetavsTPCsignalNeg(0),  \r
+    fHelium3TOF(0),                       \r
+    fhTreeESD(0)  \r
+  \r
+{\r
+  // Define input and output slots here\r
+  // Input slot #0 works with a TChain\r
+  //DefineInput(0, TChain::Class());\r
+  // Output slot #0 writes into a TList container (Cascade)\r
+  DefineOutput(1, TList::Class());\r
+}\r
+//_______________________________________________________\r
+AliAnalysisTaskHelium3Pi::~AliAnalysisTaskHelium3Pi() \r
+{ \r
+  // Destructor\r
+  if (fListHistCascade) {\r
+    delete fListHistCascade;\r
+    fListHistCascade = 0;\r
+  }\r
+\r
+}\r
+//=================DEFINITION BETHE BLOCH==============================\r
+\r
+Double_t AliAnalysisTaskHelium3Pi::BetheBloch(Double_t betaGamma,Double_t charge,Bool_t isPbPb) {\r
+\r
+  Double_t kp1, kp2, kp3, kp4, kp5;\r
+  \r
+  if(isPbPb){\r
+\r
+    //    pass1 2011\r
+    kp1 = 4.7*charge*charge;\r
+    kp2 = 8.98482806165147636e+00;\r
+    kp3 = 1.54000000000000005e-05;\r
+    kp4 = 2.30445734159456084e+00;\r
+    kp5 = 2.25624744086878559e+00;\r
+\r
+  }\r
+  \r
+  else{\r
+\r
+    // to be defined\r
+    //pass1 2011\r
+    kp1 = 4.7*charge*charge;\r
+    kp2 = 8.98482806165147636e+00;\r
+    kp3 = 1.54000000000000005e-05;\r
+    kp4 = 2.30445734159456084e+00;\r
+    kp5 = 2.25624744086878559e+00;\r
+\r
+  }\r
+\r
+  Double_t beta = betaGamma / TMath::Sqrt(1.0 + betaGamma * betaGamma);\r
+  \r
+  Double_t aa = TMath::Power(beta, kp4);\r
+  Double_t bb = TMath::Power(1.0 / betaGamma, kp5);\r
+  \r
+  bb = TMath::Log(kp3 + bb);\r
+  \r
+  Double_t out = (kp2 - aa - bb) * kp1 / aa;\r
+\r
+  return out;\r
+ \r
+}\r
+\r
+//==================DEFINITION OF OUTPUT OBJECTS==============================\r
+\r
+void AliAnalysisTaskHelium3Pi::UserCreateOutputObjects()\r
+{\r
+  fListHistCascade = new TList();\r
+  fListHistCascade->SetOwner();  // IMPORTANT!\r
+\r
+  if(! fHistEventMultiplicity ){\r
+    fHistEventMultiplicity   = new TH1F( "fHistEventMultiplicity" , "Nb of Events" , 10 , 0, 10);\r
+    //    fHistEventMultiplicity->GetXaxis()->SetTitle("Event Type");\r
+    fHistEventMultiplicity->GetXaxis()->SetBinLabel(1,"All Events");\r
+    fHistEventMultiplicity->GetXaxis()->SetBinLabel(2,"Events w/PV");\r
+    fHistEventMultiplicity->GetXaxis()->SetBinLabel(3,"Events w/|Vz|<10cm");\r
+    fHistEventMultiplicity->GetXaxis()->SetBinLabel(4,"Central Events");\r
+    fHistEventMultiplicity->GetXaxis()->SetBinLabel(5,"SemiCentral Events");\r
+    fHistEventMultiplicity->GetXaxis()->SetBinLabel(6,"MB Events");\r
+    fHistEventMultiplicity->GetXaxis()->SetBinLabel(7,"Central Events  w/|Vz|<10cm");\r
+    fHistEventMultiplicity->GetXaxis()->SetBinLabel(8,"SemiCentral Events  w/|Vz|<10cm");\r
+    fHistEventMultiplicity->GetXaxis()->SetBinLabel(9,"MB Events   w/|Vz|<10cm");\r
+\r
+    fListHistCascade->Add(fHistEventMultiplicity);\r
+  }\r
+\r
+  if(! fHistTrackMultiplicity ){\r
+    fHistTrackMultiplicity   = new TH2F( "fHistTrackMultiplicity" , "Nb of Tracks", 25000,0, 25000,105,-1,104);\r
+    fHistTrackMultiplicity->GetXaxis()->SetTitle("Number of tracks");\r
+    fHistTrackMultiplicity->GetYaxis()->SetTitle("Percentile");\r
+    fListHistCascade->Add(fHistTrackMultiplicity);\r
+  } \r
+\r
+  if(! fHistTrackMultiplicityCent ){\r
+    fHistTrackMultiplicityCent   = new TH2F( "fHistTrackMultiplicityCent" , "Nb of Tracks Central Events", 25000,0, 25000,105,-1,104 );\r
+    fHistTrackMultiplicityCent->GetXaxis()->SetTitle("Number of tracks");\r
+    fHistTrackMultiplicityCent->GetYaxis()->SetTitle("Percentile");\r
+    fListHistCascade->Add(fHistTrackMultiplicityCent);\r
+  } \r
+\r
+  if(! fHistTrackMultiplicitySemiCent ){\r
+    fHistTrackMultiplicitySemiCent   = new TH2F( "fHistTrackMultiplicitySemiCent" , "Nb of Tracks SemiCentral Events", 25000,0, 25000 ,105,-1,104);\r
+    fHistTrackMultiplicitySemiCent->GetXaxis()->SetTitle("Number of tracks");\r
+    fHistTrackMultiplicitySemiCent->GetYaxis()->SetTitle("Percentile");\r
+    fListHistCascade->Add(fHistTrackMultiplicitySemiCent);\r
+  } \r
+ \r
+  if(! fHistTrackMultiplicityMB ){\r
+    fHistTrackMultiplicityMB   = new TH2F( "fHistTrackMultiplicityMB" , "Nb of Tracks MBral Events", 25000,0, 25000,105,-1,104 );\r
+    fHistTrackMultiplicityMB->GetXaxis()->SetTitle("Number of tracks");\r
+    fHistTrackMultiplicityMB->GetYaxis()->SetTitle("Percentile");\r
+    fListHistCascade->Add(fHistTrackMultiplicityMB);\r
+  } \r
+\r
+  if(! fHistTrackMultiplicityPVCent ){\r
+    fHistTrackMultiplicityPVCent   = new TH2F( "fHistTrackMultiplicityPVCent" , "Nb of Tracks Central Events", 25000,0, 25000,105,-1,104 );\r
+    fHistTrackMultiplicityPVCent->GetXaxis()->SetTitle("Number of tracks");\r
+    fHistTrackMultiplicityPVCent->GetYaxis()->SetTitle("Percentile");\r
+    fListHistCascade->Add(fHistTrackMultiplicityPVCent);\r
+  } \r
+\r
+  if(! fHistTrackMultiplicityPVSemiCent ){\r
+    fHistTrackMultiplicityPVSemiCent   = new TH2F( "fHistTrackMultiplicityPVSemiCent" , "Nb of Tracks SemiCentral Events", 25000,0, 25000 ,105,-1,104);\r
+    fHistTrackMultiplicityPVSemiCent->GetXaxis()->SetTitle("Number of tracks");\r
+    fHistTrackMultiplicityPVSemiCent->GetYaxis()->SetTitle("Percentile");\r
+    fListHistCascade->Add(fHistTrackMultiplicityPVSemiCent);\r
+  } \r
+ \r
+  if(! fHistTrackMultiplicityPVMB ){\r
+    fHistTrackMultiplicityPVMB   = new TH2F( "fHistTrackMultiplicityPVMB" , "Nb of Tracks MBral Events", 25000,0, 25000,105,-1,104 );\r
+    fHistTrackMultiplicityPVMB->GetXaxis()->SetTitle("Number of tracks");\r
+    fHistTrackMultiplicityPVMB->GetYaxis()->SetTitle("Percentile");\r
+    fListHistCascade->Add(fHistTrackMultiplicityPVMB);\r
+  } \r
+\r
+  \r
+\r
+\r
+\r
+  if(! fHistMult){\r
+    fHistMult=new TH1F ("fHistMult","Number neg-pos", 10, -1,9);\r
+    fHistMult->GetXaxis()->SetTitle("Type of tracks");\r
+    fListHistCascade->Add(fHistMult);\r
+  }\r
+ \r
+  if(! fhBB ){\r
+    fhBB = new TH2F( "fhBB" , "BetheBlochTPC" , 1000,-6,6,1500,0,5000);\r
+    fhBB->GetXaxis()->SetTitle("p/z [GeV/c]");\r
+    fhBB->GetYaxis()->SetTitle("TPC Signal");\r
+    fListHistCascade->Add(fhBB);\r
+  }\r
+\r
+  if(! fhTOF ){\r
+    fhTOF = new TH2F( "fhTOF" , "Scatter Plot TOF" , 1000,-6,6,1000,0,1.2);\r
+    fhTOF->GetXaxis()->SetTitle("p/z [GeV/c]");\r
+    fhTOF->GetYaxis()->SetTitle("#beta");\r
+    fListHistCascade->Add(fhTOF);\r
+  }\r
+\r
+  if(! fhMassTOF){\r
+    fhMassTOF=new TH1F ("fhMassTOF","Particle Mass - TOF", 300,0 ,5);\r
+    fhMassTOF->GetXaxis()->SetTitle("Mass [GeV/c^{2}]");\r
+    fListHistCascade->Add(fhMassTOF);\r
+  }\r
+\r
+  if(! fhBBPions ){\r
+    fhBBPions = new TH2F( "fhBBPions" , "Bethe-Bloch TPC Pions" , 1000,-6,6,1500,0,5000);\r
+    fhBBPions->GetXaxis()->SetTitle("p/z [GeV/c]");\r
+    fhBBPions->GetYaxis()->SetTitle("TPC Signal");\r
+    fListHistCascade->Add(fhBBPions);\r
+  }\r
+  \r
+  if(! fhBBHe ){\r
+    fhBBHe = new TH2F( "fhBBHe" , "Bethe-Bloch TPC He" , 1000,-6,6,1500,0,5000);\r
+    fhBBHe->GetXaxis()->SetTitle("p/z [GeV/c]");\r
+    fhBBHe->GetYaxis()->SetTitle("TPC Signal");\r
+    fListHistCascade->Add(fhBBHe);\r
+  }\r
+  \r
+  if(! fhNaPos ){\r
+    fhNaPos = new TH2F( "fhNaPos" , "Distribution Pos" , 500,0,5,500,-2,2);\r
+    fhNaPos->GetXaxis()->SetTitle("p/z [GeV/c]");\r
+    fhNaPos->GetYaxis()->SetTitle("(TPCSignal-bbtheo)/bbtheo (He)");\r
+    fListHistCascade->Add(fhNaPos);\r
+  }\r
+  \r
+  if(! fhNaNeg ){\r
+    fhNaNeg = new TH2F( "fhNaNeg" , "Distribution Neg" , 500,0,5,500,-2,2);\r
+    fhNaNeg->GetXaxis()->SetTitle("p/z [GeV/c]");\r
+    fhNaNeg->GetYaxis()->SetTitle("(TPCSignal-bbtheo)/bbtheo (He)");\r
+    fListHistCascade->Add(fhNaNeg);\r
+  }\r
+\r
+  if(! fBetavsTPCsignalPos ){\r
+    fBetavsTPCsignalPos = new TH2F("fBetavsTPCsignalPos","fBetavsTPCsignalPos",1000,0,1.2,1500,0,5000);\r
+    fBetavsTPCsignalPos->GetXaxis()->SetTitle("#beta");\r
+    fBetavsTPCsignalPos->GetYaxis()->SetTitle("TPC Signal");\r
+    fListHistCascade->Add(fBetavsTPCsignalPos);\r
+  }\r
+  \r
+  if(! fBetavsTPCsignalNeg ){\r
+    fBetavsTPCsignalNeg = new TH2F("fBetavsTPCsignalNeg","fBetavsTPCsignalNeg",1000,0,1.2,1500,0,5000);\r
+    fBetavsTPCsignalNeg->GetXaxis()->SetTitle("#beta");\r
+    fBetavsTPCsignalNeg->GetYaxis()->SetTitle("TPC Signal");\r
+    fListHistCascade->Add(fBetavsTPCsignalNeg);\r
+  }\r
+  \r
+  if(! fHelium3TOF){\r
+    fHelium3TOF = new TH2F("fHelium3massTOF","Helium3 beta vs p/z",1000,0,6,1000,0,1.2);\r
+    fHelium3TOF->GetXaxis()->SetTitle("p/z [GeV/c]");\r
+    fHelium3TOF->GetYaxis()->SetTitle("#beta");\r
+    fListHistCascade->Add(fHelium3TOF);\r
+  }\r
+  \r
+  if(! fNtuple1 ) {\r
+    fNtuple1 = new TNtuple("fNtuple1","Ntuple1","runNumber:bunchcross:orbit:period:eventtype:TrackNumber:percentile:xPrimaryVertex:yPrimaryVertex:zPrimaryVertex:xSecondaryVertex:ySecondaryVertex:zSecondaryVertex:dcaTracks:CosPointingAngle:DCAV0toPrimaryVertex:HeSign:HepInTPC:HeTPCsignal:DcaHeToPrimVertex:HeEta:momHex:momHey:momHez:momHeAtSVx:momHeAtSVy:momHeAtSVz:HeTPCNcls:HeimpactXY:HeimpactZ:HeITSClusterMap:IsHeITSRefit:PionSign:PionpInTPC:PionTPCsignal:DcaPionToPrimVertex:PionEta:momPionx:momPiony:momPionz:momNegPionAtSVx:momNegPionAtSVy:momNegPionAtSVz:PionTPCNcls:PionimpactXY:PionimpactZ:PionITSClusterMap:IsPiITSRefit:xn:xp");\r
+    \r
+    fListHistCascade->Add(fNtuple1);\r
+  }\r
+  \r
+  if(! fNtuple4 ) {\r
+    fNtuple4 = new TNtuple("fNtuple4","Ntuple4","runNumber:BCNumber:OrbitNumber:PeriodNumber:eventtype:isHeITSrefit:percentile:Sign:pinTPC:GetTPCsignal:Px:Py:Pz:Eta:isTOF:poutTPC:timeTOF:trackLenghtTOF:impactXY:impactZ:mapITS:TPCNcls:TRDsignal:xPrimaryVertex:yPrimaryVertex:zPrimaryVertex:chi2PerClusterTPC");\r
+    fListHistCascade->Add(fNtuple4);\r
+  } \r
+\r
+  if(! fhTreeESD){\r
+    fhTreeESD = new TTree("fhTreeESD","Lista ESD file");\r
+    fhTreeESD->Branch("nameESD", &nameESD,25000,0); \r
+    fhTreeESD->Branch("nEvent", &nEv);\r
+    fListHistCascade->Add(fhTreeESD);\r
+  }\r
+  PostData(1,  fListHistCascade);\r
+\r
+}// end UserCreateOutputObjects\r
+\r
+\r
+\r
+//====================== USER EXEC ========================\r
+\r
+void AliAnalysisTaskHelium3Pi::UserExec(Option_t *) \r
+{\r
+  //_______________________________________________________________________\r
+  \r
+  //!*********************!//\r
+  //!  Define variables   !//\r
+  //!*********************!//\r
+  Float_t vett1[50];\r
+  for(Int_t i=0;i<50;i++) vett1[i]=0;\r
+\r
+  Float_t vett4[40];\r
+  for(Int_t i=0;i<40;i++) vett4[i]=0;\r
+  \r
+  Double_t ITSsample[4];\r
+  for(Int_t i=0;i<4;i++)ITSsample[i]=0;\r
+\r
+  Double_t  pinTPC=0.,poutTPC=0.,TPCSignal=0.;\r
+  Double_t xPrimaryVertex=0.,yPrimaryVertex=0.,zPrimaryVertex=0.;\r
+  Double_t massTOF=0.,timeTOF=0.,trackLenghtTOF=0.,betaTOF=0.;\r
+\r
+  ULong_t  status;\r
+  ULong_t  statusT;\r
+  ULong_t  statusPi;\r
+\r
+  Bool_t   isTPC,isTOF,isTOFHe,IsHeITSRefit,isTOFPi,IsPiITSRefit ;\r
+\r
+  Float_t nSigmaNegPion=0.;\r
+  Float_t nSigmaNegPion1=0.;\r
+  Float_t nSigmaNegPion2=0.;\r
+\r
+  Double_t cutNSigma = 3;\r
+  Double_t bbtheoM=0.,bbtheoP=0.,bbtheo=0.;\r
+  Double_t zNathashaNeg=0;\r
+  Double_t zNathashaPos=0;\r
+  Double_t fPos[3]={0.,0.,0.};\r
+  Double_t runNumber=0.;\r
+  Double_t evNumber=0.;\r
+\r
+  Double_t BCNumber=0.;\r
+  Double_t OrbitNumber=0.;\r
+  Double_t PeriodNumber=0.;\r
+\r
+\r
+  Double_t        Helium3Mass = 2.80839; //!OK\r
+  Double_t        PionMass = 0.13957;     //!OK\r
+\r
+  // TLORENTZ vectors\r
+  \r
+  TLorentzVector  vPion,vHelium,vSum;\r
+\r
+  //!----------------------------------------------------------------\r
+\r
+  //! A set of very loose parameters for cuts \r
+  \r
+  Double_t fgChi2max=33.;     //! max chi2\r
+  Double_t fgDNmin=0.05;      //! min imp parameter for the 1st daughter = 500um\r
+  //  Double_t fgDPmin=0.05;      //! min imp parameter for the 2nd daughter = 500um\r
+  //  Double_t fgDCAmax=1.5;  //! max DCA between the daughter tracks in cm\r
+  Double_t fgDCAmax=1.0;      //! max DCA between the daughter tracks in cm\r
+  Double_t fgCPAmin=0.99;     //! min cosine of V0's pointing angle  //qui\r
+  //Double_t fgCPAmin=-1.;    //! min cosine of V0's pointing angle\r
+  //  Double_t fgRmin=0.2;    //! min radius of the fiducial volume //original\r
+  Double_t fgRmin=0.1;        //! min radius of the fiducial volume = 1 mm \r
+  Double_t fgRmax=200.;       //! max radius of the fiducial volume = 2 m\r
+\r
+  //------------------------------------------\r
+\r
+  // Main loop\r
+  // Called for EACH event\r
+\r
+  AliVEvent *event = InputEvent();\r
+  if (!event) { Printf("ERROR: Could not retrieve event"); return; }\r
+    \r
+  Info("AliAnalysisTaskHelium3Pi","Starting UserExec");  \r
+\r
+  SetDataType("REAL");\r
+  \r
+  // create pointer to event\r
+  AliESDEvent* lESDevent = dynamic_cast<AliESDEvent*>(event);\r
+  if (!lESDevent) {\r
+    AliError("Cannot get the ESD event");\r
+    return;\r
+  }  \r
+\r
+  fHistEventMultiplicity->Fill(0);\r
+\r
+  Double_t lMagneticField=lESDevent->GetMagneticField();\r
+  Int_t TrackNumber = -1;\r
+\r
+   \r
+  //*****************//  \r
+  //*   Centrality  *//\r
+  //*****************//\r
+ \r
+  AliCentrality *centrality = lESDevent->GetCentrality();\r
+  Float_t percentile=centrality->GetCentralityPercentile("V0M");\r
+\r
+  TrackNumber = lESDevent->GetNumberOfTracks();\r
+  if (TrackNumber<2) return;  \r
+\r
+  fHistTrackMultiplicity->Fill(TrackNumber,percentile); //tracce per evento\r
+\r
+  //****************************************\r
+\r
+  Int_t eventtype=-99;\r
+  \r
+  Bool_t isSelectedCentral     = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kCentral);\r
+  Bool_t isSelectedSemiCentral = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kSemiCentral);\r
+  Bool_t isSelectedMB          = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);\r
+ \r
+  if(isSelectedCentral){\r
+    fHistEventMultiplicity->Fill(3);\r
+    fHistTrackMultiplicityCent->Fill(TrackNumber,percentile); \r
+    eventtype=1;\r
+  }\r
+\r
+  if(isSelectedSemiCentral){\r
+    fHistEventMultiplicity->Fill(4);\r
+    fHistTrackMultiplicitySemiCent->Fill(TrackNumber,percentile); \r
+    eventtype=2;\r
+  }\r
+\r
+  if(isSelectedMB){\r
+    fHistEventMultiplicity->Fill(5);\r
+    fHistTrackMultiplicityMB->Fill(TrackNumber,percentile); \r
+    eventtype=3;\r
+  }\r
+\r
+  if(isSelectedCentral || isSelectedSemiCentral || isSelectedMB){\r
+  // ANALISYS\r
+  \r
+  // Primary vertex cut\r
+  \r
+  const AliESDVertex *vtx = lESDevent->GetPrimaryVertexTracks();\r
+  \r
+  if(vtx->GetNContributors()<1) {\r
+    \r
+    // SPD vertex cut\r
+    vtx = lESDevent->GetPrimaryVertexSPD();\r
+    \r
+    if(vtx->GetNContributors()<1) {\r
+      Info("AliAnalysisTaskHelium3Pi","No good vertex, skip event");\r
+      return; // NO GOOD VERTEX, SKIP EVENT \r
+    }\r
+  }\r
+\r
+  fHistEventMultiplicity->Fill(1); // analyzed events with PV\r
+ \r
+  xPrimaryVertex=vtx->GetXv();\r
+  yPrimaryVertex=vtx->GetYv();\r
+  zPrimaryVertex=vtx->GetZv();  \r
+\r
+  if(TMath::Abs(zPrimaryVertex)>10) return;\r
+\r
+  if(eventtype==1){\r
+    fHistTrackMultiplicityPVCent->Fill(TrackNumber,percentile); \r
+    fHistEventMultiplicity->Fill(6); \r
+  }\r
+\r
+  if(eventtype==2){\r
+    fHistTrackMultiplicityPVSemiCent->Fill(TrackNumber,percentile); \r
+    fHistEventMultiplicity->Fill(7); \r
+  }\r
+\r
+  if(eventtype==3){\r
+    fHistTrackMultiplicityPVMB->Fill(TrackNumber,percentile); \r
+    fHistEventMultiplicity->Fill(8); \r
+  }\r
+\r
+\r
+  fHistEventMultiplicity->Fill(2);\r
+\r
+  //Find Pair candidates\r
+  \r
+  TArrayI PionsTPC(TrackNumber);        //Neg pions\r
+  Int_t nPionsTPC=0;\r
+  \r
+  TArrayI HeTPC(TrackNumber);        //helium3\r
+  Int_t nHeTPC=0;\r
+\r
+  const Double_t speedOfLight =  TMath::C()*1E2*1E-12; // cm/ps\r
+  \r
+  Float_t impactXY=-999, impactZ=-999;\r
+  Float_t impactXYpi=-999, impactZpi=-999;\r
+\r
+  //!   SELECTIONS:\r
+  //! - No ITSpureSA\r
+  //! - ITSrefit\r
+  //! - TPCrefit\r
+  //! - ITSmap !=0\r
+\r
+  // ******************* Track Cuts Definitions ********************//\r
+  \r
+  //! ITS\r
+\r
+  AliESDtrackCuts* esdtrackCutsITS = new AliESDtrackCuts("esdtrackCutsITS");\r
+  esdtrackCutsITS->SetRequireITSStandAlone(kFALSE);\r
+  esdtrackCutsITS->SetRequireITSPureStandAlone(kFALSE);\r
+\r
+  //! TPC\r
+  \r
+  Int_t    minclsTPC=60;\r
+  Double_t maxchi2perTPCcl=5.;\r
+  \r
+  AliESDtrackCuts* esdtrackCutsTPC = new AliESDtrackCuts("esdtrackCutsTPC");\r
+  esdtrackCutsTPC->SetRequireTPCRefit(kTRUE);\r
+  esdtrackCutsTPC->SetAcceptKinkDaughters(kFALSE);\r
+  esdtrackCutsTPC->SetMinNClustersTPC(minclsTPC);\r
+  esdtrackCutsTPC->SetMaxChi2PerClusterTPC(maxchi2perTPCcl);\r
+  \r
+  //*********************************************+\r
+ \r
+  runNumber = lESDevent->GetRunNumber();\r
+  evNumber  = lESDevent->GetEventNumberInFile();\r
+      \r
+  BCNumber    = lESDevent->GetBunchCrossNumber();\r
+  OrbitNumber = lESDevent->GetOrbitNumber();\r
+  PeriodNumber= lESDevent->GetPeriodNumber();\r
+\r
+  //*************************************************************\r
+  \r
+  TF1 *foPion=new TF1("foPion", "[0]*([1]*TMath::Power(TMath::Sqrt(1 + (x/0.13957)*(x/0.13957))/(x/0.13957) , [3]) - 1 - TMath::Power(TMath::Sqrt(1 + (x/0.13957)*(x/0.13957))/(x/0.13957) , [3])*TMath::Log([2] + 1/TMath::Power((x/0.13957), [4])))",0.01,20);\r
+  foPion->SetParameters(4.1,8.98482806165147636e+00,1.54000000000000005e-05,2.30445734159456084e+00,2.25624744086878559e+00);\r
+\r
+  //--------------------------------------------\r
+\r
+  for (Int_t j=0; j<TrackNumber; j++) { //loop on tracks\r
+    \r
+    AliESDtrack *esdtrack=lESDevent->GetTrack(j);\r
+   \r
+    if(!esdtrack) { \r
+      AliError(Form("ERROR: Could not retrieve esdtrack %d",j)); \r
+      continue; \r
+    }\r
+\r
+    // ************** Track cuts ****************\r
+    \r
+    status  = (ULong_t)esdtrack->GetStatus();\r
+    \r
+    isTPC   = ((status & AliESDtrack::kTPCin)  != 0);\r
+    isTOF   = (((status & AliESDtrack::kTOFout) != 0) && ((status & AliESDtrack::kTIME) != 0));\r
+    \r
+    Bool_t IsTrackAcceptedTPC =  esdtrackCutsTPC->AcceptTrack(esdtrack);\r
+    Bool_t IsTrackAcceptedITS =  esdtrackCutsITS->AcceptTrack(esdtrack);\r
+\r
+    if (!(IsTrackAcceptedTPC && IsTrackAcceptedITS)) continue;\r
+\r
+    UInt_t mapITS=esdtrack->GetITSClusterMap();\r
+    //    if(mapITS==0) continue;\r
+\r
+    //----------------------------------------------\r
+    \r
+    //****** Cuts from  AliV0Vertex.cxx *************\r
+     \r
+    Double_t d=esdtrack->GetD(xPrimaryVertex,yPrimaryVertex,lMagneticField);\r
+    //    if (TMath::Abs(d)<fgDPmin) continue;\r
+    if (TMath::Abs(d)>fgRmax) continue;\r
+    \r
+    //---- (Usefull) Stuff\r
+    \r
+    TPCSignal=esdtrack->GetTPCsignal(); \r
+    \r
+    if (TPCSignal<10)continue;\r
+    \r
+    if(!isTPC)continue;\r
+\r
+    if(!esdtrack->GetTPCInnerParam())continue;\r
+    \r
+    AliExternalTrackParam trackIn(*esdtrack->GetInnerParam()); \r
+    pinTPC = trackIn.GetP(); \r
+        \r
+    poutTPC=pinTPC;\r
+    \r
+    fHistMult->Fill(0);\r
+   \r
+    if(status & AliESDtrack::kITSrefit!=0){\r
+      fHistMult->Fill(1);\r
+      fhBB->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);\r
+    }\r
+    \r
+    timeTOF=esdtrack->GetTOFsignal();                 // ps\r
+    trackLenghtTOF= esdtrack->GetIntegratedLength();  // cm\r
+\r
+    if(isTOF){\r
+      \r
+      if(!esdtrack->GetOuterParam())continue;    \r
+      \r
+      AliExternalTrackParam trackOut(*esdtrack->GetOuterParam()); \r
+      \r
+      poutTPC = trackOut.GetP(); \r
+      \r
+      betaTOF= (trackLenghtTOF/timeTOF)/2.99792458e-2;\r
+      \r
+      fhTOF->Fill(poutTPC*esdtrack->GetSign(),betaTOF);\r
+\r
+      Double_t mass2=(poutTPC*poutTPC)*((((speedOfLight*speedOfLight)*(timeTOF*timeTOF))-(trackLenghtTOF*trackLenghtTOF))/(trackLenghtTOF*trackLenghtTOF));\r
+      if(mass2>0) massTOF=TMath::Sqrt(mass2);\r
+      fhMassTOF->Fill(massTOF);\r
+\r
+      if(esdtrack->GetSign() < 0.)fBetavsTPCsignalNeg->Fill(betaTOF,TPCSignal);\r
+      if(esdtrack->GetSign() > 0.)fBetavsTPCsignalPos->Fill(betaTOF,TPCSignal);\r
+         \r
+    }\r
+\r
+    //pass2\r
+    \r
+    bbtheo =BetheBloch((2*pinTPC)/3.,2,kTRUE);    //! OK\r
+    bbtheoM=(1 - 0.08*5)*bbtheo;                  //! OK \r
+    bbtheoP=(1 + 0.08*5)*bbtheo;                  //! OK\r
+    \r
+    fHistMult->Fill(2);\r
+    \r
+    if(esdtrack->GetSign()<0){\r
+      zNathashaNeg=(TPCSignal-bbtheo)/bbtheo;\r
+      fhNaNeg->Fill(pinTPC,zNathashaNeg); \r
+    }\r
+\r
+    if(esdtrack->GetSign() > 0.){\r
+      zNathashaPos=(TPCSignal-bbtheo)/bbtheo;\r
+      fhNaPos->Fill(pinTPC,zNathashaPos); \r
+    }\r
+    \r
+    nSigmaNegPion1 = TMath::Abs((TPCSignal - foPion->Eval(pinTPC))/foPion->Eval(pinTPC))/0.07; \r
+    nSigmaNegPion2 = TMath::Abs((TPCSignal - foPion->Eval(pinTPC))/foPion->Eval(pinTPC))/0.04; \r
+    \r
+    if(pinTPC<0.6)\r
+      nSigmaNegPion=nSigmaNegPion1;\r
+    if(pinTPC>=0.6)\r
+      nSigmaNegPion=nSigmaNegPion2;\r
+    \r
+    if ( (nSigmaNegPion < cutNSigma)){ \r
+      \r
+      fhBBPions->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);\r
+\r
+      if(pinTPC<3){\r
+       PionsTPC[nPionsTPC++]=j;\r
+      }\r
+    }\r
+    \r
+    if( TPCSignal > bbtheoM ) {\r
+ \r
+      if(pinTPC>0.6){\r
+\r
+       fhBBHe->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);\r
+       HeTPC[nHeTPC++]=j;\r
+       \r
+       nameESD=(TString)AliAnalysisTaskSE::CurrentFileName();\r
+       nameESD.Data();\r
+       nEv=(Int_t)AliAnalysisTaskSE::Entry();\r
+       fhTreeESD->Fill();\r
+       \r
+       Bool_t isHeITSrefit=(status & AliESDtrack::kITSrefit);\r
+       \r
+       esdtrack->GetImpactParameters(impactXY, impactZ);\r
+       esdtrack->GetITSdEdxSamples(ITSsample);\r
+       \r
+       \r
+       Int_t  fIdxInt[200]; //dummy array\r
+       Int_t nClustersTPC = esdtrack->GetTPCclusters(fIdxInt);\r
+       \r
+       Float_t chi2PerClusterTPC = esdtrack->GetTPCchi2()/(Float_t)(nClustersTPC);\r
+       \r
+       vett4[0] =(Float_t)runNumber;\r
+       vett4[1] =(Float_t)BCNumber;\r
+       vett4[2] =(Float_t)OrbitNumber;\r
+       vett4[3] =(Float_t)PeriodNumber;\r
+       vett4[4] =(Float_t)eventtype;\r
+       vett4[5] =(Float_t)isHeITSrefit;\r
+       vett4[6] =(Float_t)percentile;\r
+       vett4[7] =(Float_t)esdtrack->GetSign();\r
+       vett4[8] =(Float_t)pinTPC;\r
+       vett4[9] =(Float_t)esdtrack->GetTPCsignal();\r
+       vett4[10]=(Float_t)esdtrack->Px();\r
+       vett4[11]=(Float_t)esdtrack->Py();\r
+       vett4[12]=(Float_t)esdtrack->Pz();\r
+       vett4[13]=(Float_t)esdtrack->Eta();\r
+       vett4[14]=(Float_t)isTOF;\r
+       vett4[15]=(Float_t)poutTPC;\r
+       vett4[16]=(Float_t)timeTOF;\r
+       vett4[17]=(Float_t)trackLenghtTOF;\r
+       vett4[18]=(Float_t)impactXY;\r
+       vett4[19]=(Float_t)impactZ;\r
+       vett4[20]=(Float_t)mapITS;\r
+       vett4[21]=(Float_t)esdtrack->GetTPCNcls();\r
+       vett4[22]=(Float_t)esdtrack->GetTRDsignal();\r
+       vett4[23]=(Float_t)xPrimaryVertex;\r
+       vett4[24]=(Float_t)yPrimaryVertex;\r
+       vett4[25]=(Float_t)zPrimaryVertex;\r
+       vett4[26]=(Float_t)chi2PerClusterTPC;\r
+       \r
+       fNtuple4->Fill(vett4);\r
+      }\r
+    }\r
+  }       //! track\r
+  \r
+  \r
+  Double_t        DcaHeToPrimVertex=0;\r
+  Double_t        DcaPionToPrimVertex=0;\r
+  \r
+  impactXY=-999, impactZ=-999;\r
+  impactXYpi=-999, impactZpi=-999;\r
+  \r
+  // Track \r
+\r
+  AliESDtrack  *PionTrack = 0x0;\r
+  AliESDtrack  *HeTrack = 0x0;\r
+\r
+  // Vettors for il PxPyPz\r
+\r
+  Double_t momPionVett[3];\r
+  for(Int_t i=0;i<3;i++)momPionVett[i]=0;\r
+  \r
+  Double_t momHeVett[3];\r
+  for(Int_t i=0;i<3;i++)momHeVett[i]=0;\r
+\r
+  //At SV\r
+\r
+  Double_t momPionVettAt[3];\r
+  for(Int_t i=0;i<3;i++)momPionVettAt[i]=0;\r
+  \r
+  Double_t momHeVettAt[3];\r
+  for(Int_t i=0;i<3;i++)momHeVettAt[i]=0;\r
+\r
+  //---------------   LOOP PAIRS   ----------------\r
+\r
+  for (Int_t k=0; k < nPionsTPC; k++) {                           //! Pions Loop\r
+\r
+    DcaPionToPrimVertex=0.;\r
+    DcaHeToPrimVertex=0;\r
+\r
+    Int_t PionIdx=PionsTPC[k];\r
+    \r
+    PionTrack=lESDevent->GetTrack(PionIdx);\r
+\r
+    statusPi = (ULong_t)PionTrack->GetStatus();\r
+    isTOFPi  = (((statusPi & AliESDtrack::kTOFout) != 0) && ((statusPi & AliESDtrack::kTIME) != 0));\r
+    IsPiITSRefit = (statusPi & AliESDtrack::kITSrefit); \r
+\r
+    if (PionTrack) \r
+      DcaPionToPrimVertex = TMath::Abs(PionTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK\r
+   \r
+    if(DcaPionToPrimVertex<0.2)continue; //qui\r
+\r
+    AliExternalTrackParam trackInPion(*PionTrack);  \r
+  \r
+    for (Int_t i=0; i<nHeTPC; i++){                               //! Helium Loop\r
+\r
+      Int_t HeIdx=HeTPC[i];\r
+      \r
+      HeTrack=lESDevent->GetTrack(HeIdx);\r
+     \r
+      statusT= (ULong_t)HeTrack->GetStatus();\r
+      isTOFHe   = (((statusT & AliESDtrack::kTOFout) != 0) && ((statusT & AliESDtrack::kTIME) != 0));\r
+      IsHeITSRefit = (status & AliESDtrack::kITSrefit); \r
+ \r
+     if (HeTrack) \r
+       DcaHeToPrimVertex = TMath::Abs(HeTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK\r
+    \r
+      AliExternalTrackParam trackInHe(*HeTrack); \r
+    \r
+      if ( DcaPionToPrimVertex < fgDNmin)                //OK\r
+       if ( DcaHeToPrimVertex < fgDNmin) continue;    //OK\r
+      \r
+      Double_t xn, xp;\r
+      Double_t dca=0.;\r
+      \r
+      dca= PionTrack->GetDCA(HeTrack,lMagneticField,xn,xp); //!dca (Neg to Pos)\r
+\r
+      if (dca > fgDCAmax) continue;\r
+      if ((xn+xp) > 2*fgRmax) continue;\r
+      if ((xn+xp) < 2*fgRmin) continue;\r
+\r
+      //CORREZIONE da AliV0Vertex\r
+      \r
+      Bool_t corrected=kFALSE;\r
+      if ((trackInPion.GetX() > 3.) && (xn < 3.)) {\r
+       //correct for the beam pipe material\r
+       corrected=kTRUE;\r
+      }\r
+      if ((trackInHe.GetX() > 3.) && (xp < 3.)) {\r
+       //correct for the beam pipe material\r
+       corrected=kTRUE;\r
+      }\r
+      if (corrected) {\r
+       dca=trackInPion.GetDCA(&trackInHe,lMagneticField,xn,xp);\r
+       if (dca > fgDCAmax) continue;\r
+       if ((xn+xp) > 2*fgRmax) continue;\r
+       if ((xn+xp) < 2*fgRmin) continue;\r
+      }\r
+ \r
+      //=============================================//\r
+      // Faccio un "V0" con le tracce che ho trovato //\r
+      //=============================================//\r
+\r
+      trackInPion.PropagateTo(xn,lMagneticField); \r
+      trackInHe.PropagateTo(xp,lMagneticField);\r
+            \r
+      AliESDv0 vertex(trackInPion,PionIdx,trackInHe,HeIdx);\r
+      if (vertex.GetChi2V0() > fgChi2max) continue;\r
+      \r
+      Float_t CosPointingAngle=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex); //PointingAngle\r
+      if (CosPointingAngle < fgCPAmin) continue;\r
+\r
+      vertex.SetDcaV0Daughters(dca);\r
+      vertex.SetV0CosineOfPointingAngle(CosPointingAngle);\r
+\r
+      fPos[0]=vertex.Xv();\r
+      fPos[1]=vertex.Yv(); \r
+      fPos[2]=vertex.Zv(); \r
+      \r
+      HeTrack->PxPyPz(momHeVett);\r
+      PionTrack->PxPyPz(momPionVett); \r
+      \r
+      Double_t raggio=TMath::Sqrt(fPos[0]*fPos[0]+fPos[1]*fPos[1]+fPos[2]*fPos[2]);\r
+      HeTrack->GetPxPyPzAt(raggio,lMagneticField,momHeVettAt);\r
+      PionTrack->GetPxPyPzAt(raggio,lMagneticField,momPionVettAt); \r
+\r
+      //------------------------------------------------------------------------//\r
+\r
+      HeTrack->GetImpactParameters(impactXY, impactZ);\r
+      \r
+      PionTrack->GetImpactParameters(impactXYpi, impactZpi);\r
+      \r
+      if(vertex.GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex)>3) continue;\r
+\r
+      //salvo solo fino a 3.1 GeV/c2\r
+\r
+      vHelium.SetXYZM(2*momHeVettAt[0],2*momHeVettAt[1],2*momHeVettAt[2],Helium3Mass); \r
+      vPion.SetXYZM(momPionVettAt[0],momPionVettAt[1],momPionVettAt[2],PionMass);       \r
+      vSum=vHelium+vPion;\r
+      \r
+      if(vSum.M()>3.1)continue;\r
+\r
+      //----------------------------------------------------------------------//\r
+      \r
+      vett1[0]=(Float_t)runNumber;\r
+      vett1[1]=(Float_t)BCNumber;\r
+      vett1[2]=(Float_t)OrbitNumber;\r
+      vett1[3]=(Float_t)PeriodNumber;\r
+      vett1[4]=(Float_t)eventtype;\r
+      vett1[5]=(Float_t)TrackNumber;\r
+      vett1[6]=(Float_t)percentile;\r
+      vett1[7]=(Float_t)xPrimaryVertex; //PRIMARY\r
+      vett1[8]=(Float_t)yPrimaryVertex;\r
+      vett1[9]=(Float_t)zPrimaryVertex;\r
+      vett1[10]=(Float_t)fPos[0]; //SECONDARY\r
+      vett1[11]=(Float_t)fPos[1];\r
+      vett1[12]=(Float_t)fPos[2];\r
+      vett1[13]=(Float_t)dca;           //between 2 tracks\r
+      vett1[14]=(Float_t)CosPointingAngle;          //cosPointingAngle da V0\r
+      vett1[15]=(Float_t)vertex.GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);\r
+      vett1[16]=(Float_t)HeTrack->GetSign(); //He\r
+      vett1[17]=(Float_t)trackInHe.GetP();\r
+      vett1[18]=(Float_t)HeTrack->GetTPCsignal();\r
+      vett1[19]=(Float_t)DcaHeToPrimVertex;\r
+      vett1[20]=(Float_t)HeTrack->Eta();\r
+      vett1[21]=(Float_t)momHeVett[0];\r
+      vett1[22]=(Float_t)momHeVett[1];\r
+      vett1[23]=(Float_t)momHeVett[2];\r
+      vett1[24]=(Float_t)momHeVettAt[0];\r
+      vett1[25]=(Float_t)momHeVettAt[1];\r
+      vett1[26]=(Float_t)momHeVettAt[2];\r
+      vett1[27]=(Float_t)HeTrack->GetTPCNcls();\r
+      vett1[28]=(Float_t)impactXY;\r
+      vett1[29]=(Float_t)impactZ;\r
+      vett1[30]=(Float_t)HeTrack->GetITSClusterMap();\r
+      vett1[31]=(Float_t)IsHeITSRefit;\r
+      vett1[32]=(Float_t)PionTrack->GetSign(); //Pion\r
+      vett1[33]=(Float_t)trackInPion.GetP();\r
+      vett1[34]=(Float_t)PionTrack->GetTPCsignal();\r
+      vett1[35]=(Float_t)DcaPionToPrimVertex;\r
+      vett1[36]=(Float_t)PionTrack->Eta();\r
+      vett1[37]=(Float_t)momPionVett[0];\r
+      vett1[38]=(Float_t)momPionVett[1];\r
+      vett1[39]=(Float_t)momPionVett[2];\r
+      vett1[40]=(Float_t)momPionVettAt[0];\r
+      vett1[41]=(Float_t)momPionVettAt[1];\r
+      vett1[42]=(Float_t)momPionVettAt[2];\r
+      vett1[43]=(Float_t)PionTrack->GetTPCNcls();\r
+      vett1[44]=(Float_t)impactXYpi;\r
+      vett1[45]=(Float_t)impactZpi;\r
+      vett1[46]=(Float_t)PionTrack->GetITSClusterMap();\r
+      vett1[47]=(Float_t)IsPiITSRefit;\r
+      vett1[48]=(Float_t)xn;\r
+      vett1[49]=(Float_t)xp;\r
+\r
+      fNtuple1->Fill(vett1);  \r
+\r
+    }// positive TPC\r
+    \r
+  } //negative tpc\r
+  \r
+  }\r
+\r
+  PostData(1,fListHistCascade);\r
+    \r
+} //end userexec\r
+\r
+\r
+//________________________________________________________________________\r
+\r
+void AliAnalysisTaskHelium3Pi::Terminate(Option_t *) \r
+{\r
+  // Draw result to the screen\r
+  // Called once at the end of the query\r
+}\r
+\r
+\r

diff --git a/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3Pi.h b/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3Pi.h
new file mode 100644 (file)
index 0000000..0491e05
--- /dev/null
+++ b/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3Pi.h
@@ -0,0 +1,94 @@
+#ifndef ALIANALYSISTASKHELIUM3PI_H
+#define ALIANALYSISTASKHELIUM3PI_H
+
+/*  See cxx source for full Copyright notice */
+
+//-----------------------------------------------------------------
+//                 AliAnalysisTaskHelium3Pion class
+//-----------------------------------------------------------------
+
+class TList;
+class TH1F;
+class TH2F;
+class TH3F;
+class TNtuple;
+class AliESDcascade;
+//class AliCascadeVertexer; 
+
+#include "TString.h"
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskHelium3Pi : public AliAnalysisTaskSE {
+ public:
+  AliAnalysisTaskHelium3Pi();
+  AliAnalysisTaskHelium3Pi(const char *name);
+  virtual ~AliAnalysisTaskHelium3Pi();
+  
+  virtual void  UserCreateOutputObjects();
+  virtual void  UserExec(Option_t *option);
+  virtual void  Terminate(Option_t *);
+  
+  void SetCollidingSystems(Short_t collidingSystems = 0)     {fCollidingSystems = collidingSystems;}
+  void SetAnalysisType    (const char* analysisType = "ESD") {fAnalysisType = analysisType;}
+  void SetDataType    (const char* dataType = "REAL") {fDataType = dataType;}
+  void SetMinNSigma(Double_t ns=3.){
+    fMinNSigma=ns;
+  }
+
+  void SetMinSPDPoints(Int_t np=1){
+    fMinSPDPts=np;
+  }
+
+  //Double_t BetheBloch(Double_t bg,Double_t Charge,Bool_t optMC);
+  Double_t BetheBloch(Double_t bg,Double_t Charge,Bool_t isPbPb);
+
+  
+ private:
+
+  TString fAnalysisType;            //! "ESD" or "AOD" analysis type   
+
+  Short_t fCollidingSystems;        //! 0 = pp collisions or 1 = AA collisions
+  TString fDataType;                //! "REAL" or "SIM" data type      
+  TList        *fListHistCascade;           //! List of Cascade histograms
+  TH1F *fHistEventMultiplicity;
+  TH2F *fHistTrackMultiplicity;
+  TH2F *fHistTrackMultiplicityCent;
+  TH2F *fHistTrackMultiplicitySemiCent;
+  TH2F *fHistTrackMultiplicityMB;
+  TH2F *fHistTrackMultiplicityPVCent;
+  TH2F *fHistTrackMultiplicityPVSemiCent;
+  TH2F *fHistTrackMultiplicityPVMB;
+  TH1F *fHistMult;
+  TH2F *fhBB;
+  TH2F *fhTOF;
+  TH1F *fhMassTOF;
+  TH2F *fhBBPions;
+  TH2F *fhBBHe;
+  TH2F *fhNaPos;
+  TH2F *fhNaNeg;
+  TH2F *fBetavsTPCsignalPos;
+  TH2F *fBetavsTPCsignalNeg;
+  TH2F *fHelium3TOF;
+   
+  TNtuple *fNtuple1;                  //! Ntupla Pairs Pi/Proton "standard"
+  TNtuple *fNtuple4;                  //! He carateristiche
+
+  TTree *fhTreeESD;                  //! Tree di Tstring
+  TString nameESD ; 
+  Int_t nEv;
+
+  Int_t    fMinSPDPts;       // minimum number of SPD Points
+  Double_t fMinNSigma;       // minimum number of sigmas
+  
+
+ static const Int_t fgNrot;
+ 
+
+  AliAnalysisTaskHelium3Pi(const AliAnalysisTaskHelium3Pi&);            // not implemented
+  AliAnalysisTaskHelium3Pi& operator=(const AliAnalysisTaskHelium3Pi&); // not implemented
+  
+  ClassDef(AliAnalysisTaskHelium3Pi, 0);
+};
+
+#endif

diff --git a/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3PiMC.cxx b/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3PiMC.cxx
new file mode 100644 (file)
index 0000000..273b0e1
--- /dev/null
+++ b/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3PiMC.cxx
@@ -0,0 +1,1413 @@
+/**************************************************************************\r
+ * Contributors are not mentioned at all.                                 *\r
+ *                                                                        *\r
+ * Permission to use, copy, modify and distribute this software and its   *\r
+ * documentation strictly for non-commercial purposes is hereby granted   *\r
+ * without fee, provided that the above copyright notice appears in all   *\r
+ * copies and that both the copyright notice and this permission notice   *\r
+ * appear in the supporting documentation. The authors make no claims     *\r
+ * about the suitability of this software for any purpose. It is          *\r
+ * provided "as is" without express or implied warranty.                  *\r
+ **************************************************************************/\r
+//\r
+//-----------------------------------------------------------------\r
+//                 AliAnalysisTaskHelium3Pi class\r
+//-----------------------------------------------------------------\r
+\r
+\r
+class TTree;\r
+class TParticle;\r
+class TVector3;\r
+\r
+#include "AliAnalysisManager.h"\r
+#include <AliMCEventHandler.h>\r
+#include <AliMCEvent.h>\r
+#include <AliStack.h>\r
+\r
+class AliESDVertex;\r
+class AliAODVertex;\r
+class AliESDv0;\r
+class AliAODv0; \r
+class AliCascadeVertexer;\r
+\r
+#include <iostream>\r
+#include "AliAnalysisTaskSE.h"\r
+#include "TList.h"\r
+#include "TH1.h"\r
+#include "TH2.h"\r
+#include "TH3.h"\r
+#include "TNtuple.h"\r
+#include "TGraph.h"\r
+#include "TCutG.h"\r
+#include "TF1.h"\r
+#include "TCanvas.h"\r
+#include "TMath.h"\r
+#include "TChain.h"\r
+#include "Riostream.h"\r
+#include "AliLog.h"\r
+#include "AliCascadeVertexer.h"\r
+#include "AliESDEvent.h"\r
+#include "AliESDtrack.h"\r
+#include "AliExternalTrackParam.h"\r
+#include "AliAODEvent.h"\r
+#include "AliInputEventHandler.h"\r
+#include "AliESDcascade.h"\r
+#include "AliAODcascade.h"\r
+#include "AliAnalysisTaskHelium3PiMC.h"\r
+#include "AliESDtrackCuts.h"\r
+#include "AliCentrality.h"\r
+#include "TString.h"\r
+#include <TDatime.h>\r
+#include <TRandom3.h>\r
+\r
+const Int_t AliAnalysisTaskHelium3PiMC::fgNrot = 15;\r
+\r
+\r
+ClassImp(AliAnalysisTaskHelium3PiMC)\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskHelium3PiMC::AliAnalysisTaskHelium3PiMC() \r
+: AliAnalysisTaskSE(),\r
+  fAnalysisType("ESD"), \r
+  fCollidingSystems(0), \r
+  fDataType("SIM"),\r
+  fListHistCascade(0), \r
+  fHistEventMultiplicity(0), \r
+  fHistTrackMultiplicity(0),\r
+  fHistMCMultiplicityTracks(0),\r
+  fHistMCEta(0), \r
+  fHistMCPt(0), \r
+  fHistMCTheta(0), \r
+  fHistMCDecayPosition(0),\r
+  fHistMCDecayRho(0), \r
+  fhRigidityHevsMomPiMC(0),\r
+  fhRigidityHevsMomPiRec(0),\r
+  fhInvMassMC(0),\r
+  fhInvMassMum(0),\r
+  fhInvMassRec(0),\r
+  fhInvMassRec1(0),\r
+  fhInvMassRec2(0), \r
+  fhInvMassRec3(0),\r
+  fhInvMassRec4(0),\r
+  fhInvMassRec5(0),\r
+  fhInvMassRec6(0),\r
+  fhInvMassRec7(0),\r
+  fhHeMCRigidity(0),\r
+  fhPioneMC(0),\r
+  hBBTPCnoCuts(0),\r
+  fhBBTPC(0),\r
+  fhBBTPCNegativePions(0),\r
+  fhBBTPCPositivePions(0),\r
+  fhBBTPCHe3(0),\r
+  fHistProvaDCA(0),\r
+  fHistPercentileVsTrackNumber(0),\r
+  fhHeDCAXY(0),\r
+  fhHeDCAZ(0),\r
+  fhPiDCAXY(0),\r
+  fhPiDCAZ(0),\r
+  hITSClusterMap(0)\r
+\r
+{\r
+  // Dummy Constructor(0); \r
+}\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskHelium3PiMC::AliAnalysisTaskHelium3PiMC(const char *name) \r
+  : AliAnalysisTaskSE(name), \r
+    fAnalysisType("ESD"), \r
+    fCollidingSystems(0), \r
+    fDataType("SIM"),\r
+    fListHistCascade(0), \r
+    fHistEventMultiplicity(0), \r
+    fHistTrackMultiplicity(0),\r
+    fHistMCMultiplicityTracks(0),\r
+    fHistMCEta(0), \r
+    fHistMCPt(0), \r
+    fHistMCTheta(0), \r
+    fHistMCDecayPosition(0),\r
+    fHistMCDecayRho(0), \r
+    fhRigidityHevsMomPiMC(0),\r
+    fhRigidityHevsMomPiRec(0),\r
+    fhInvMassMC(0),\r
+    fhInvMassMum(0),\r
+    fhInvMassRec(0),\r
+    fhInvMassRec1(0),\r
+    fhInvMassRec2(0), \r
+    fhInvMassRec3(0),\r
+    fhInvMassRec4(0),\r
+    fhInvMassRec5(0),\r
+    fhInvMassRec6(0),\r
+    fhInvMassRec7(0),\r
+    fhHeMCRigidity(0),\r
+    fhPioneMC(0),\r
+    hBBTPCnoCuts(0),\r
+    fhBBTPC(0),\r
+    fhBBTPCNegativePions(0),\r
+    fhBBTPCPositivePions(0),\r
+    fhBBTPCHe3(0),\r
+    fHistProvaDCA(0),\r
+    fHistPercentileVsTrackNumber(0),\r
+    fhHeDCAXY(0),\r
+    fhHeDCAZ(0),\r
+    fhPiDCAXY(0),\r
+    fhPiDCAZ(0),\r
+    hITSClusterMap(0)\r
+\r
+\r
+{\r
+  // Define input and output slots here\r
+  // Input slot #0 works with a TChain\r
+  //DefineInput(0, TChain::Class());\r
+  // Output slot #0 writes into a TList container (Cascade)\r
+  DefineOutput(1, TList::Class());\r
+}\r
+//_______________________________________________________\r
+AliAnalysisTaskHelium3PiMC::~AliAnalysisTaskHelium3PiMC() \r
+{ \r
+  // Destructor\r
+  if (fListHistCascade) {\r
+    delete fListHistCascade;\r
+    fListHistCascade = 0;\r
+  }\r
+\r
+}\r
+//=================DEFINITION BETHE BLOCH==============================\r
+\r
+Double_t AliAnalysisTaskHelium3PiMC::BetheBloch(Double_t betaGamma,Double_t charge,Bool_t isPbPb) {\r
+\r
+  Double_t kp1, kp2, kp3, kp4, kp5;\r
+  \r
+  if(isPbPb){\r
+\r
+    //pass2\r
+    kp1 = 5.2*charge*charge;\r
+    kp2 = 8.98482806165147636e+00;\r
+    kp3 = 1.54000000000000005e-05;\r
+    kp4 = 2.30445734159456084e+00;\r
+    kp5 = 2.25624744086878559e+00;\r
+\r
+//     //pass1\r
+//     kp1 = 1.25202*charge*charge;///50.; \r
+//     kp2 = 2.74992e+01; \r
+//     kp3 = TMath::Exp(-3.31517e+01); \r
+//     kp4 = 2.46246; \r
+//     kp5 = 6.78938;\r
+\r
+  }\r
+  \r
+  else{\r
+\r
+    //pass2\r
+    kp1 = 5.2*charge*charge;\r
+    kp2 = 8.98482806165147636e+00;\r
+    kp3 = 1.54000000000000005e-05;\r
+    kp4 = 2.30445734159456084e+00;\r
+    kp5 = 2.25624744086878559e+00;\r
+\r
+    //     //pass1\r
+    //     kp1 = 1.25202*charge*charge;///50.; \r
+    //     kp2 = 2.74992e+01; \r
+    //     kp3 = TMath::Exp(-3.31517e+01); \r
+    //     kp4 = 2.46246; \r
+    //     kp5 = 6.78938;\r
+    \r
+    //     kp1 = 4.23232575531564326e+00*charge*charge;  //50*0.76176e-1;\r
+    //     kp2 = 8.68482806165147636e+00;                //10.632; \r
+    //     kp3 = 1.34000000000000005e-05;                //0.13279e-4;\r
+    //     kp4 = 2.30445734159456084e+00;                //1.8631;\r
+    //     kp5 = 2.25624744086878559e+00;                //1.9479;\r
+  }\r
+\r
+  Double_t beta = betaGamma / TMath::Sqrt(1.0 + betaGamma * betaGamma);\r
+  \r
+  Double_t aa = TMath::Power(beta, kp4);\r
+  Double_t bb = TMath::Power(1.0 / betaGamma, kp5);\r
+  \r
+  bb = TMath::Log(kp3 + bb);\r
+  \r
+  Double_t out = (kp2 - aa - bb) * kp1 / aa;\r
+\r
+  return out;\r
+ \r
+}\r
+\r
+//==================DEFINITION OF OUTPUT OBJECTS==============================\r
+\r
+void AliAnalysisTaskHelium3PiMC::UserCreateOutputObjects()\r
+{\r
+  fListHistCascade = new TList();\r
+  fListHistCascade->SetOwner();  // IMPORTANT!\r
+\r
+  if(! fHistEventMultiplicity ){\r
+    fHistEventMultiplicity   = new TH1F( "fHistEventMultiplicity" , "Nb of Events" , 6 , -1, 5 );\r
+    fHistEventMultiplicity->GetXaxis()->SetTitle("Event Type");\r
+    fListHistCascade->Add(fHistEventMultiplicity);\r
+  }\r
+\r
+  if(! fHistTrackMultiplicity ){\r
+\r
+    fHistTrackMultiplicity   = new TH1F( "fHistTrackMultiplicity" , "Nb of Tracks" , 25000,0, 25000 );\r
+    fHistTrackMultiplicity->GetXaxis()->SetTitle("Number of tracks");\r
+    fListHistCascade->Add(fHistTrackMultiplicity);\r
+  } \r
+ \r
+  if(! fHistMCMultiplicityTracks){ \r
+    fHistMCMultiplicityTracks =new TH1F("fHistMCMultiplicityTracks","MC Multiplicity Tracks",1000,0,1000); \r
+    fHistMCMultiplicityTracks->GetXaxis()->SetTitle("MC Number of tracks");\r
+    fListHistCascade->Add(fHistMCMultiplicityTracks); \r
+  }\r
+  if(!fHistMCEta ){ \r
+    fHistMCEta=new TH1F("fHistMCEta","MC eta",1000,-3,3);                                     \r
+    fHistMCEta->GetXaxis()->SetTitle("Injected Eta");\r
+    fListHistCascade->Add(fHistMCEta);\r
+  } \r
+  if(! fHistMCPt){ \r
+    fHistMCPt =new TH1F("fHistMCPt","MC pt",1000,0,20);     \r
+    fHistMCPt->GetXaxis()->SetTitle("Injected Pt");\r
+    fListHistCascade->Add(fHistMCPt); \r
+  }  \r
+  if(!fHistMCTheta ){ \r
+    fHistMCTheta=new TH1F("fHistMCTheta","MC theta",1000,-6,6);                                   \r
+    fHistMCTheta->GetXaxis()->SetTitle("Injected Theta");\r
+    fListHistCascade->Add(fHistMCTheta); \r
+  }\r
+  if(!fHistMCDecayPosition){ \r
+    fHistMCDecayPosition =new TH1F("fHistMCDecayPosition","MC Decay Position",10000,0,1000);  \r
+    fHistMCDecayPosition->GetXaxis()->SetTitle("Decay Position");\r
+    fListHistCascade->Add(fHistMCDecayPosition); \r
+  }    \r
+  if(!fHistMCDecayRho ){ \r
+    fHistMCDecayRho  =new TH1F("fHistMCDecayRho","MC decay position 3d",10000,0,1000);  \r
+    fHistMCDecayRho->GetXaxis()->SetTitle("Decay rho");\r
+    fListHistCascade->Add(fHistMCDecayRho); \r
+  } \r
+\r
+  if(!fhRigidityHevsMomPiMC ){ \r
+    fhRigidityHevsMomPiMC=new TH2F("fhRigidityHevsMomPiMC","Rigidity He vs Mom Pi MC",20,0,10,300,0,30);\r
+    fhRigidityHevsMomPiMC->GetXaxis()->SetTitle("He3 Rigidity");\r
+    fhRigidityHevsMomPiMC->GetYaxis()->SetTitle("Pi momentum");\r
+    fListHistCascade->Add(fhRigidityHevsMomPiMC); \r
+  }\r
+\r
+  if(! fhRigidityHevsMomPiRec){ \r
+    fhRigidityHevsMomPiRec=new TH2F("fhRigidityHevsMomPiRec","Rigidity He vs Mom Pi Rec",20,0,10,300,0,30);\r
+    fhRigidityHevsMomPiRec->GetXaxis()->SetTitle("He3 Rigidity");\r
+    fhRigidityHevsMomPiRec->GetYaxis()->SetTitle("Pi momentum");\r
+    fListHistCascade->Add(fhRigidityHevsMomPiRec); \r
+  }\r
+  \r
+  if(!fhInvMassMC){\r
+    fhInvMassMC=new TH1F("fhInvMassMC","fhInvMassMC",800,2.,6.);\r
+    fhInvMassMC->GetXaxis()->SetTitle("(He3,#pi) InvMass");\r
+    fListHistCascade->Add(fhInvMassMC); \r
+  }\r
+  \r
+  if(!fhInvMassMum){\r
+    fhInvMassMum=new TH1F("fhInvMassMum","fhInvMassMum",800,2.,6.);\r
+    fhInvMassMum->GetXaxis()->SetTitle("(He3,#pi) InvMass");\r
+    fListHistCascade->Add(fhInvMassMum); \r
+  }\r
+  \r
+  if(!fhInvMassRec){\r
+    fhInvMassRec=new TH1F("fhInvMassRec","fhInvMassRec",800,2.,6.);\r
+    fhInvMassRec->GetXaxis()->SetTitle("(He3,#pi) InvMass");\r
+    fListHistCascade->Add(fhInvMassRec);\r
+  }\r
+\r
+  if(!fhInvMassRec1){\r
+    fhInvMassRec1=new TH1F("fhInvMassRec1","No Altri tagli",800,2.,6.);\r
+    fhInvMassRec1->GetXaxis()->SetTitle("(He3,#pi) InvMass");\r
+    fListHistCascade->Add(fhInvMassRec1);\r
+  }\r
+  if(!fhInvMassRec2){\r
+    fhInvMassRec2=new TH1F("fhInvMassRec2","DCA pi > 0.1",800,2.,6.);\r
+    fhInvMassRec2->GetXaxis()->SetTitle("(He3,#pi) InvMass");\r
+    fListHistCascade->Add(fhInvMassRec2);\r
+  }\r
+  if(!fhInvMassRec3){\r
+    fhInvMassRec3=new TH1F("fhInvMassRec3","DCA He > 0.05",800,2.,6.);\r
+    fhInvMassRec3->GetXaxis()->SetTitle("(He3,#pi) InvMass");\r
+    fListHistCascade->Add(fhInvMassRec3);\r
+  }\r
+  if(!fhInvMassRec4){\r
+    fhInvMassRec4=new TH1F("fhInvMassRec4","DCA tracks < 1 cm",800,2.,6.);\r
+    fhInvMassRec4->GetXaxis()->SetTitle("(He3,#pi) InvMass");\r
+    fListHistCascade->Add(fhInvMassRec4);\r
+  }\r
+  if(!fhInvMassRec5){\r
+    fhInvMassRec5=new TH1F("fhInvMassRec5","Condizione xn+xp",800,2.,6.);\r
+    fhInvMassRec5->GetXaxis()->SetTitle("(He3,#pi) InvMass");\r
+    fListHistCascade->Add(fhInvMassRec5);\r
+  }\r
+  if(!fhInvMassRec6){\r
+    fhInvMassRec6=new TH1F("fhInvMassRec6","Ho fatto V0 ",800,2.,6.);\r
+    fhInvMassRec6->GetXaxis()->SetTitle("(He3,#pi) InvMass");\r
+    fListHistCascade->Add(fhInvMassRec6);\r
+  }\r
+  if(!fhInvMassRec7){\r
+    fhInvMassRec7=new TH1F("fhInvMassRec7","V0+Taglio CPA",800,2.,6.);\r
+    fhInvMassRec7->GetXaxis()->SetTitle("(He3,#pi) InvMass");\r
+    fListHistCascade->Add(fhInvMassRec7);\r
+  }\r
+\r
+  if(!fhHeMCRigidity ){ \r
+    fhHeMCRigidity=new TH1F("fhHeMCRigidity","He3 rigidity distribution",200,0,20);\r
+    fhHeMCRigidity->GetXaxis()->SetTitle("He3 rigidity");\r
+    fListHistCascade->Add(fhHeMCRigidity); \r
+  }\r
+  if(!fhPioneMC ){ \r
+    fhPioneMC=new TH1F("hPioneMC","Pion mom distribution",200,0,50);\r
+    fhPioneMC->GetXaxis()->SetTitle("Pion momentum");\r
+    fListHistCascade->Add(fhPioneMC); \r
+  }\r
+  \r
+  if(!hBBTPCnoCuts ){ \r
+    hBBTPCnoCuts=new TH2F("hBBTPCnoCuts","scatterPlot TPC no cuts",2000,-10,10,1000,0,3000);\r
+    hBBTPCnoCuts->GetXaxis()->SetTitle("p/Z (GeV/c)");\r
+    hBBTPCnoCuts->GetYaxis()->SetTitle("TPC Signal (a.u)");\r
+    fListHistCascade->Add(hBBTPCnoCuts); \r
+  }\r
+  if(!fhBBTPC ){ \r
+    fhBBTPC=new TH2F("fhBBTPC","scatterPlot TPC",2000,-10,10,1000,0,3000);\r
+    fhBBTPC->GetXaxis()->SetTitle("p/Z (GeV/c)");\r
+    fhBBTPC->GetYaxis()->SetTitle("TPC Signal (a.u)");\r
+    fListHistCascade->Add(fhBBTPC); \r
+  }\r
+  if(!fhBBTPCNegativePions ){ \r
+    fhBBTPCNegativePions=new TH2F("fhBBTPCNegativePions","scatterPlot Neg Pions",2000,-10,10,1000,0,3000);\r
+    fhBBTPCNegativePions->GetXaxis()->SetTitle("p/Z (GeV/c)");\r
+    fhBBTPCNegativePions->GetYaxis()->SetTitle("TPC Signal (a.u)");\r
+    fListHistCascade->Add(fhBBTPCNegativePions); \r
+  }\r
+  if(!fhBBTPCPositivePions ){ \r
+    fhBBTPCPositivePions=new TH2F("fhBBTPCPositivePions","scatterPlot Pos Pions",2000,-10,10,1000,0,3000);\r
+    fhBBTPCPositivePions->GetXaxis()->SetTitle("p/Z (GeV/c)");\r
+    fhBBTPCPositivePions->GetYaxis()->SetTitle("TPC Signal (a.u)");\r
+    fListHistCascade->Add(fhBBTPCPositivePions); \r
+  }\r
+  if(!fhBBTPCHe3 ){ \r
+    fhBBTPCHe3=new TH2F("fhBBTPCHe3","scatterPlot TPC -  He3",2000,-10,10,1000,0,3000);\r
+    fhBBTPCHe3->GetXaxis()->SetTitle("p/Z (GeV/c)");\r
+    fhBBTPCHe3->GetYaxis()->SetTitle("TPC Signal (a.u)");\r
+    fListHistCascade->Add(fhBBTPCHe3); \r
+  }\r
+  if(!fHistProvaDCA ){ \r
+    fHistProvaDCA=new TH2F("fHistProvaDCA","fHistProvaDCA",1000,-50,50,1000,0,100);\r
+    fHistProvaDCA->GetXaxis()->SetTitle("xn+xp");\r
+    fHistProvaDCA->GetYaxis()->SetTitle("dca tracks");\r
+    fListHistCascade->Add(fHistProvaDCA); \r
+  }\r
+  \r
+  if(!hITSClusterMap ){ \r
+    hITSClusterMap=new TH1F("hITSClusterMap","hITSClusterMap",65,-1,64);\r
+    fListHistCascade->Add(hITSClusterMap); \r
+  }\r
+\r
+  if(!fHistPercentileVsTrackNumber){\r
+    fHistPercentileVsTrackNumber=new TH2F("fHistPercentileVsTrackNumber","fHistPercentileVsTrackNumber",120,-3,117,2500,0,25000);\r
+    fHistPercentileVsTrackNumber->GetXaxis()->SetTitle("Percentile");\r
+    fHistPercentileVsTrackNumber->GetYaxis()->SetTitle("Tracks Number");\r
+    fListHistCascade->Add(fHistPercentileVsTrackNumber); \r
+  }\r
+\r
+  if(!fhHeDCAXY){\r
+    fhHeDCAXY=new TH1F("fhHeDCAXY","fhHeDCAXY",800,-4,4);\r
+    fListHistCascade->Add(fhHeDCAXY); \r
+  }\r
+  if(!fhHeDCAZ){\r
+    fhHeDCAZ=new TH1F("fhHeDCAZ","fhHeDCAZ",800,-30,30);\r
+    fListHistCascade->Add(fhHeDCAZ); \r
+  }\r
+  if(!fhPiDCAXY){\r
+    fhPiDCAXY=new TH1F("fhPiDCAXY","fhPiDCAXY",800,-4,4);\r
+    fListHistCascade->Add(fhPiDCAXY); \r
+  }\r
+  if(!fhPiDCAZ){\r
+    fhPiDCAZ=new TH1F("fhPiDCAZ","fhPiDCAZ",800,-30,30);\r
+    fListHistCascade->Add(fhPiDCAZ); \r
+  }\r
+\r
+  if(! fNtuple1 ) {\r
+    fNtuple1 = new TNtuple("fNtuple1","Ntuple1","runNumber:evNumber:TrackNumber:percentile:xPrimaryVertex:yPrimaryVertex:zPrimaryVertex:xSecondaryVertex:ySecondaryVertex:zSecondaryVertex:dcaTracks:CosPointingAngle:DCAV0toPrimaryVertex:HeSign:HepInTPC:HeTPCsignal:DcaHeToPrimVertex:HeEta:momHex:momHey:momHez:momHeAtSVx:momHeAtSVy:momHeAtSVz:HeTPCNcls:HeimpactXY:HeimpactZ:isTOFHe:HeBeta:HeITSClusterMap:IsHeITSRefit:PionSign:PionpInTPC:PionTPCsignal:DcaPionToPrimVertex:PionEta:momPionx:momPiony:momPionz:momNegPionAtSVx:momNegPionAtSVy:momNegPionAtSVz:PionTPCNcls:PionimpactXY:PionimpactZ:isTOFPion:PionBeta:PionITSClusterMap:IsPiITSRefit:PDGCodeNeg:PDCCodePos:motherPDGNeg:motherPDGPos:labelPi:labelHe:mumidNeg:mumidPos");\r
+  \r
+    //    fNtuple1->SetDirectory(0);\r
+    fListHistCascade->Add(fNtuple1);\r
+  }\r
+  \r
+  if(! fNtuple2 ) {\r
+    \r
+    fNtuple2 = new TNtuple("fNtuple2","Ntuple2","run:event:iMC:Centrality:PVx:PVy:PVz:PDGcodeMum:MotherIndex:SVxD0:SVyD0:SVzD0:SVxD1:SVyD1:SVzD1:SV3d:EtaMum:YMum:ThetaMum:PhiMum:PxMum:PyMum:PzMum:PdgDaughter0:PdgDaughter1:PxD0:PyD0:PzD0:PxD1:PyD1:PzD1");\r
+    \r
+    //fNtuple2 = new TNtuple("fNtuple2","Ntuple2","run:event:iMC:Centrality:PxHeMc:PyHeMc:PzHeMc:PxPionMc:PyPionMc:PzPionMc:fhInvMassMC");\r
+    //fNtuple2->SetDirectory(0);\r
+    fListHistCascade->Add(fNtuple2);\r
+  } \r
+    \r
+  PostData(1,fListHistCascade);\r
+\r
+}// end UserCreateOutputObjects\r
+\r
+\r
+\r
+//====================== USER EXEC ========================\r
+\r
+void AliAnalysisTaskHelium3PiMC::UserExec(Option_t *) \r
+{\r
+  //_______________________________________________________________________\r
+  \r
+  //  cout<<"Inside the event loop"<<endl;\r
+\r
+\r
+  //!*********************!//\r
+  //!  Define variables   !//\r
+  //!*********************!//\r
+  Float_t vett1[60];\r
+  for(Int_t i=0;i<60;i++) vett1[i]=0;\r
+\r
+  Float_t vett2[40];\r
+  for(Int_t i=0;i<40;i++) vett2[i]=0;\r
+\r
+  Float_t vett3[40];\r
+  for(Int_t i=0;i<40;i++) vett3[i]=0;\r
+\r
+  Float_t vett4[40];\r
+  for(Int_t i=0;i<40;i++) vett4[i]=0;\r
+  \r
+  Double_t ITSsample[4];\r
+  for(Int_t i=0;i<4;i++)ITSsample[i]=0;\r
+\r
+  Double_t ITSsamplePos[4];\r
+  for(Int_t i=0;i<4;i++)ITSsamplePos[i]=0;\r
+\r
+  Double_t ITSsampleNeg[4];\r
+  for(Int_t i=0;i<4;i++)ITSsampleNeg[i]=0;\r
+\r
+  Double_t  pinTPC=0.,poutTPC=0.,TPCSignal=0.;\r
+  Double_t xPrimaryVertex=0.,yPrimaryVertex=0.,zPrimaryVertex=0.;\r
+\r
+  ULong_t  status;\r
+  ULong_t  statusT;\r
+  ULong_t  statusPi;\r
+\r
+  Bool_t   isTPC,isTOF,isTOFHe3,isTOFPi;\r
+  //  Bool_t   isTOFPos,isTOFNeg,isTOFPosPion;\r
+\r
+  Double_t fPos[3]={0.,0.,0.};\r
+  Double_t runNumber=0.;\r
+  Double_t evNumber=0.;\r
+ \r
+  Int_t id0           = 0, id1          = 0;\r
+  Double_t mcDecayPosXD0 = 0, mcDecayPosYD0 = 0, mcDecayPosR = 0, mcDecayPosZD0 =0, mcDecayPosRho=0.;\r
+  Double_t mcDecayPosXD1 = 0, mcDecayPosYD1 = 0, mcDecayPosZD1 =0;\r
+\r
+  Double_t lEtaCurrentPart =0., lPtCurrentPart = 0.,lThetaCurrentPart = 0., lPhiCurrentPart = 0.;\r
+  Int_t iCurrentMother = 0;\r
+  Double_t mcPosX = 0.,        mcPosY = 0.,mcPosZ = 0., mcPosR = 0.;\r
+\r
+  Double_t lPdgCurrentDaughter0 = 0, lPdgCurrentDaughter1= 0., lPdgCurrentMother=0.,lPdgCurrentDaughter =0;\r
+\r
+  Double_t PxD0 = 0, PyD0 = 0,PzD0 = 0;\r
+  Double_t PxD1 = 0, PyD1 = 0,PzD1 = 0;\r
+\r
+  Int_t lNbMCPrimary        = 0;\r
+  Int_t lNbMCPart           = 0;\r
+  Int_t lPdgcodeCurrentPart = 0;\r
+  //!----------------------------------------------------------------\r
+\r
+  //! A set of very loose parameters for cuts \r
+  \r
+  Double_t fgChi2max=33.;     //! max chi2\r
+  Double_t fgDNmin=0.05;      //! min imp parameter for the 1st daughter = 500um\r
+  //  Double_t fgDPmin=0.05;      //! min imp parameter for the 2nd daughter = 500um\r
+  //  Double_t fgDCAmax=1.5;  //! max DCA between the daughter tracks in cm\r
+  Double_t fgDCAmax=1.;       //! max DCA between the daughter tracks in cm\r
+  Double_t fgCPAmin=0.9;      //! min cosine of V0's pointing angle\r
+  //  Double_t fgCPAmin=-1.;      //! min cosine of V0's pointing angle\r
+  //  Double_t fgRmin=0.2;    //! min radius of the fiducial volume //original\r
+  Double_t fgRmin=0.1;        //! min radius of the fiducial volume = 1 mm \r
+  Double_t fgRmax=200.;       //! max radius of the fiducial volume = 2 m\r
+\r
+  //------------------------------------------\r
+  // create pointer to event\r
+\r
+  AliVEvent *event = InputEvent();\r
+  if (!event) { Printf("ERROR: Could not retrieve event"); return; }\r
+\r
+\r
+\r
+//   AliVEvent *lESDevent = InputEvent();\r
+//   if (!lESDevent) { \r
+//     Printf("ERROR: Could not retrieve event"); \r
+//     return; \r
+//  }\r
+    \r
+  Info("AliAnalysisTaskHelium3PiMC","Starting UserExec");  \r
+\r
+  SetDataType("SIM");\r
+  AliStack *stack;\r
+  if(fDataType == "SIM") {\r
+  // Main loop\r
+  // Called for EACH event\r
+  AliMCEvent *mcEvent = MCEvent();\r
+  if (!mcEvent) { \r
+    Printf("ERROR: Could not retrieve MC event"); \r
+    return; \r
+  }\r
+  \r
+  Printf("MC particles: %d", mcEvent->GetNumberOfTracks());\r
+  \r
+  // set up a stack for use in check for primary/stable particles\r
+  stack = mcEvent->Stack();\r
+  if( !stack ) { Printf( "Stack not available"); return; }\r
+  }\r
+ \r
+\r
+  AliESDEvent *lESDevent = 0x0;\r
+\r
+  //********************************** Connect to the InputEvent ******//\r
+\r
+  //Int_t TrackNumber = 0;\r
+  if(fAnalysisType == "ESD"){\r
+    lESDevent = dynamic_cast<AliESDEvent*>(event);\r
+    if (!lESDevent) {\r
+      Printf("ERROR: lESDevent not available \n");\r
+      return;\r
+    }\r
+    //TrackNumber = lESDevent->GetNumberOfTracks();\r
+  }\r
+\r
+  //*****************//  \r
+  //*   Centrality  *//\r
+  //*****************//\r
+\r
+  AliCentrality *centrality = lESDevent->GetCentrality();\r
+ \r
+  Float_t percentile=centrality->GetCentralityPercentile("V0M");\r
+  \r
+  //------------------------------\r
+\r
+  runNumber = lESDevent->GetRunNumber();\r
+  evNumber =lESDevent->GetEventNumberInFile();\r
+\r
+  //---------------------\r
+\r
+  //  Int_t primary = stack->GetNprimary();\r
+  Int_t label =-1;\r
+\r
+  lNbMCPrimary = stack->GetNprimary();\r
+  lNbMCPart    = stack->GetNtrack();\r
+      \r
+  fHistMCMultiplicityTracks->Fill(lNbMCPart);     //histo\r
+\r
+  TArrayD MomPionsMC(lNbMCPart);        //Neg pions\r
+  Int_t nPionsMC=0;\r
+  TArrayD MomHeMC(lNbMCPart);        //helium3\r
+  Int_t nHeMC=0;\r
+\r
+  //------ Trimomento pion\r
+  TArrayD PxPionsMC(lNbMCPart);       \r
+  Int_t nPxPionsMC=0;\r
+  TArrayD PyPionsMC(lNbMCPart);       \r
+  Int_t nPyPionsMC=0;\r
+  TArrayD PzPionsMC(lNbMCPart);       \r
+  Int_t nPzPionsMC=0;\r
+  //------ Trimomento He\r
+  TArrayD PxHeMC(lNbMCPart);       \r
+  Int_t nPxHeMC=0;\r
+  TArrayD PyHeMC(lNbMCPart);       \r
+  Int_t nPyHeMC=0;\r
+  TArrayD PzHeMC(lNbMCPart);       \r
+  Int_t nPzHeMC=0;\r
+\r
+  //Mass Definition\r
+\r
+  Double_t        Helium3Mass = 2.80839; //!OK\r
+  Double_t        PionMass = 0.13957;     //!OK\r
+  \r
+  TLorentzVector  vPionMC,vHeliumMC,vSumMC;\r
+  TLorentzVector  vPionMum,vHeliumMum,vSumMum;\r
+  TLorentzVector  vPionRec,vHeliumRec,vSumRec;\r
+  Bool_t isTwoBody=kFALSE;\r
+  for (Int_t iMC=0; iMC<stack->GetNtrack(); iMC++)\r
+    {\r
+      TParticle *p0 = stack->Particle(iMC);\r
+\r
+      if (!p0) {\r
+       //Printf("ERROR: particle with label %d not found in stack (mc loop)", iMc);\r
+       continue;\r
+      }\r
+      \r
+      \r
+      lPdgcodeCurrentPart  = p0->GetPdgCode(); \r
+     \r
+      if(lPdgcodeCurrentPart == 1000020030 || lPdgcodeCurrentPart == -1000020030 ){\r
+\r
+      // if(lPdgcodeCurrentPart == 1000020030){\r
+\r
+//     lPoscodePart == lPdgcodeCurrentPart;\r
+\r
+//     Int_t lmumidPos = lPdgcodeCurrentPart->GetFirstMother();\r
+//     if(lmumidPos>-1){\r
+//       TParticle *pmotherPos=(TParticle*)stack->Particle(lmumidPos);\r
+//       lmotherPDGPos = pmotherPos->GetPdgCode();\r
+//     }\r
+\r
+       MomHeMC[nHeMC++]=p0->P();\r
+       \r
+        PxHeMC[nPxHeMC++]=p0->Px();\r
+       PyHeMC[nPyHeMC++]=p0->Py();\r
+       PzHeMC[nPzHeMC++]=p0->Pz();\r
+       \r
+       fhHeMCRigidity->Fill(p0->P()/2);\r
+      }\r
+\r
+      if(lPdgcodeCurrentPart == 211 || lPdgcodeCurrentPart == -211 ){\r
+\r
+      //      if(lPdgcodeCurrentPart == -211 ){\r
+\r
+//     lNegcodePart = lPdgcodeCurrentPart;\r
+\r
+//     Int_t lmumidNeg = lPdgcodeCurrentPart->GetFirstMother();\r
+//     if(lmumidNeg>-1){\r
+//       TParticle *pmotherNeg=(TParticle*)stack->Particle(lmumidNeg);\r
+//       lmotherPDGNeg = pmotherNeg->GetPdgCode();\r
+//     }\r
+\r
+       MomPionsMC[nPionsMC++]=p0->P();\r
+\r
+       PxPionsMC[nPxPionsMC++]=p0->Px();\r
+       PyPionsMC[nPyPionsMC++]=p0->Py();\r
+       PzPionsMC[nPzPionsMC++]=p0->Pz();\r
+\r
+       fhPioneMC->Fill(p0->P());\r
+      }\r
+\r
+      if ( lPdgcodeCurrentPart == 1010010030 || lPdgcodeCurrentPart == -1010010030 ){\r
+\r
+      //   if (lPdgcodeCurrentPart == -1010010030){\r
+       \r
+       //      fhSimulatedMultiplicity->Fill(1);\r
+       lEtaCurrentPart   = p0->Eta();\r
+       lPtCurrentPart    = p0->Pt();\r
+       lThetaCurrentPart = p0->Theta();\r
+       lPhiCurrentPart   = p0->Phi();\r
+       iCurrentMother    = p0->GetFirstMother();\r
+       \r
+       fHistMCEta->Fill(lEtaCurrentPart);    //histo\r
+       fHistMCPt->Fill(lPtCurrentPart);      //histo\r
+       fHistMCTheta->Fill(lThetaCurrentPart);//histo\r
+       \r
+       //        lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();\r
+       if (iCurrentMother == -1){lPdgCurrentMother=0; } else {lPdgCurrentMother = stack->Particle(iCurrentMother)->GetPdgCode();}\r
+        \r
+       mcPosX = p0->Vx();\r
+       mcPosY = p0->Vy();\r
+       mcPosZ = p0->Vz();\r
+       mcPosR = TMath::Sqrt(mcPosX*mcPosX+mcPosY*mcPosY);\r
+\r
+//     cout<<"ndaughters: "<<p0->GetNDaughters()<<endl;\r
+//     cout<<"First daughters: "<<p0->GetFirstDaughter()<<endl;\r
+//     cout<<"Last daughters: "<<p0->GetLastDaughter()<<endl;\r
+\r
+       isTwoBody=kFALSE;\r
+\r
+       for(Int_t i=p0->GetFirstDaughter(); i<= p0->GetLastDaughter(); i++){\r
+         TParticle *pDaughter = stack->Particle(i);\r
+         lPdgCurrentDaughter= pDaughter->GetPdgCode();\r
+         cout<<lPdgCurrentDaughter<<endl;\r
+         if(lPdgCurrentDaughter == 1000020030 || lPdgCurrentDaughter ==-1000020030 ){\r
+           //  if(lPdgCurrentDaughter == 1000020030){\r
+           //  cout<<"\nHelio 3!!!\n"<<endl;\r
+           isTwoBody=kTRUE;\r
+           //  cout<<"Is 2 body? inside "<<isTwoBody<<endl;\r
+         }\r
+       }\r
+       \r
+       if(isTwoBody){\r
+\r
+         for(Int_t i=p0->GetFirstDaughter(); i<= p0->GetLastDaughter(); i++){\r
+         \r
+           TParticle *pDaughter = stack->Particle(i);\r
+         \r
+           lPdgCurrentDaughter= pDaughter->GetPdgCode();\r
+           //    cout<<"i "<<i<<" "<<lPdgCurrentDaughter<<endl;\r
+           if(lPdgCurrentDaughter == 211 || lPdgCurrentDaughter == -211 ){\r
+             // if(lPdgCurrentDaughter == 211 ){\r
+             id0 = i;\r
+           }\r
+           \r
+           if(lPdgCurrentDaughter == 1000020030 || lPdgCurrentDaughter == -1000020030 ){\r
+             //if(lPdgCurrentDaughter == -1000020030){\r
+             id1 = i;\r
+           }\r
+         }\r
+\r
+         TParticle *pDaughter0 = stack->Particle(id0);\r
+         TParticle *pDaughter1 = stack->Particle(id1);\r
+         lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();\r
+         lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();\r
+         \r
+         // Decay Radius and Production Radius\r
+         \r
+         if ( id0 <= lNbMCPart && id0 > 0 && id1 <= lNbMCPart && id1 > 0) {\r
+           \r
+           lPdgCurrentDaughter0 = pDaughter0->GetPdgCode();\r
+           lPdgCurrentDaughter1 = pDaughter1->GetPdgCode();\r
+           \r
+           PxD0 = pDaughter0->Px();\r
+           PyD0 = pDaughter0->Py();\r
+           PzD0 = pDaughter0->Pz();\r
+           \r
+           PxD1 = pDaughter1->Px();\r
+           PyD1 = pDaughter1->Py();\r
+           PzD1 = pDaughter1->Pz();\r
+           \r
+           mcDecayPosXD0 = pDaughter0->Vx();\r
+           mcDecayPosYD0 = pDaughter0->Vy();\r
+           mcDecayPosZD0 = pDaughter0->Vz();\r
+           \r
+           mcDecayPosXD1 = pDaughter0->Vx();\r
+           mcDecayPosYD1 = pDaughter0->Vy();\r
+           mcDecayPosZD1 = pDaughter0->Vz();\r
+           \r
+           mcDecayPosR = TMath::Sqrt(mcDecayPosXD0*mcDecayPosXD0+mcDecayPosYD0*mcDecayPosYD0);\r
+           fHistMCDecayPosition->Fill(mcDecayPosR);  //histo\r
+           \r
+           //mcDecayPosZ = pDaughter0->Vz();\r
+           mcDecayPosRho = TMath::Sqrt(mcDecayPosXD0*mcDecayPosXD0+mcDecayPosYD0*mcDecayPosYD0+mcDecayPosZD0*mcDecayPosZD0);\r
+           fHistMCDecayRho->Fill(mcDecayPosRho);    //histo\r
+           \r
+           //---- Massa iniziale prova\r
+           \r
+           vHeliumMum.SetXYZM(PxD1,PyD1,PzD1,Helium3Mass); \r
+           vPionMum.SetXYZM(PxD0,PyD0,PzD0,PionMass);       \r
+           vSumMum=vHeliumMum+vPionMum;\r
+           \r
+           fhInvMassMum->Fill(vSumMum.M());\r
+           \r
+           //Ntupla hyper triton\r
+           \r
+           vett2[0]=(Float_t)lESDevent->GetRunNumber();\r
+           vett2[1]=(Float_t)lESDevent->GetEventNumberInFile();\r
+           vett2[2]=(Float_t)iMC;\r
+           vett2[3]=(Float_t)percentile;\r
+           vett2[4]=(Float_t)mcPosX;\r
+           vett2[5]=(Float_t)mcPosY;\r
+           vett2[6]=(Float_t)mcPosZ;\r
+           vett2[7]=(Float_t)lPdgcodeCurrentPart;\r
+           vett2[8]=(Float_t)iCurrentMother;\r
+           vett2[9]=(Float_t)mcDecayPosXD0;\r
+           vett2[10]=(Float_t)mcDecayPosYD0;\r
+           vett2[11]=(Float_t)mcDecayPosZD0;\r
+           vett2[12]=(Float_t)mcDecayPosXD1;\r
+           vett2[13]=(Float_t)mcDecayPosYD1;\r
+           vett2[14]=(Float_t)mcDecayPosZD1;\r
+           vett2[15]=(Float_t)mcDecayPosRho;\r
+           vett2[16]=(Float_t)lEtaCurrentPart;\r
+           vett2[17]=(Float_t)p0->Y();\r
+           vett2[18]=(Float_t)lThetaCurrentPart;\r
+           vett2[19]=(Float_t)lPhiCurrentPart;\r
+           vett2[20]=(Float_t)p0->Px();\r
+           vett2[21]=(Float_t)p0->Py();\r
+           vett2[22]=(Float_t)p0->Pz();\r
+           vett2[23]=(Float_t)lPdgCurrentDaughter0;\r
+           vett2[24]=(Float_t)lPdgCurrentDaughter1;\r
+           vett2[25]=(Float_t)PxD0; //pion\r
+           vett2[26]=(Float_t)PyD0;\r
+           vett2[27]=(Float_t)PzD0;\r
+           vett2[28]=(Float_t)PxD1; //He3\r
+           vett2[29]=(Float_t)PyD1;\r
+           vett2[30]=(Float_t)PzD1;\r
+           \r
+           fNtuple2->Fill(vett2);\r
+           \r
+         }//if check daughters index\r
+       }//is twobody\r
+      } // if mother\r
+    } // Kinetic Track loop ends here \r
+    \r
+  // Loop phase - space\r
+\r
+  // NB ora e' solo He3+pi meno\r
+\r
+  Double_t HeMomMC =0;\r
+  Double_t PionMomMC=0;\r
+  Double_t PxHeMc=0, PyHeMc=0,PzHeMc=0;\r
+  Double_t PxPionMc=0, PyPionMc=0,PzPionMc=0;\r
+\r
+  for(Int_t l=0; l < nHeMC; l++){\r
+      \r
+    HeMomMC=MomHeMC[l];\r
+\r
+    PxHeMc=PxHeMC[l];\r
+    PyHeMc=PyHeMC[l];\r
+    PzHeMc=PzHeMC[l];\r
+\r
+    for(Int_t k=0; k < nPionsMC; k++){\r
+   \r
+      PionMomMC=MomPionsMC[k];\r
+      \r
+      PxPionMc=PxPionsMC[k];\r
+      PyPionMc=PyPionsMC[k];\r
+      PzPionMc=PzPionsMC[k];\r
+      \r
+      fhRigidityHevsMomPiMC->Fill(HeMomMC/2,PionMomMC);\r
+\r
+      vHeliumMC.SetXYZM(PxHeMc,PyHeMc,PzHeMc,Helium3Mass); \r
+      vPionMC.SetXYZM(PxPionMc,PyPionMc,PzPionMc,PionMass);       \r
+      vSumMC=vHeliumMC+vPionMC;\r
+      \r
+      fhInvMassMC->Fill(vSumMC.M());\r
+\r
+    }\r
+    \r
+  } // fine loop phase space\r
+\r
+  //-------------- Inizia la parte ricostruita -------------------\r
+\r
+  fHistEventMultiplicity->Fill(0);\r
+  \r
+  Double_t lMagneticField=lESDevent->GetMagneticField();\r
+  //  cout<<"lMagneticField: "<<lMagneticField<<endl;\r
+  Int_t TrackNumber = -1;\r
+\r
+  // ANALISYS reconstructed tracks\r
+  \r
+  // Primary vertex cut\r
+  \r
+  const AliESDVertex *vtx = lESDevent->GetPrimaryVertexTracks();\r
+  \r
+  if(vtx->GetNContributors()<1) {\r
+    \r
+    // SPD vertex cut\r
+    vtx = lESDevent->GetPrimaryVertexSPD();\r
+    \r
+    if(vtx->GetNContributors()<1) {\r
+      Info("AliAnalysisTaskHelium3PiMC","No good vertex, skip event");\r
+      return; // NO GOOD VERTEX, SKIP EVENT \r
+    }\r
+  }\r
+\r
+  fHistEventMultiplicity->Fill(1); // analyzed events with PV\r
+ \r
+  xPrimaryVertex=vtx->GetXv();\r
+  yPrimaryVertex=vtx->GetYv();\r
+  zPrimaryVertex=vtx->GetZv();  \r
+\r
+  TrackNumber = lESDevent->GetNumberOfTracks();\r
+  fHistTrackMultiplicity->Fill(TrackNumber); //tracce per evento\r
+\r
+  fHistPercentileVsTrackNumber->Fill(percentile,TrackNumber);\r
+\r
+  if (TrackNumber<2) return;  \r
+\r
+  fHistEventMultiplicity->Fill(2);\r
+    \r
+  //Find Pair candidates\r
+    \r
+  TArrayI PionsTPC(TrackNumber);        //Neg pions\r
+  Int_t nPionsTPC=0;\r
+  \r
+  TArrayI HeTPC(TrackNumber);        //helium3\r
+  Int_t nHeTPC=0;\r
+ \r
+  // find pairs candidates phase daughter tracks rec\r
+\r
+  TArrayD MomPionsRec(TrackNumber);        //Neg pions\r
+  Int_t nPionsRec=0;\r
+  \r
+  TArrayD MomHeRec(TrackNumber);        //helium3\r
+  Int_t nHeRec=0;\r
+\r
+  //------ Trimomento pion\r
+  TArrayD PxPionsRec(TrackNumber);       \r
+  Int_t nPxPionsRec=0;\r
+  TArrayD PyPionsRec(TrackNumber);       \r
+  Int_t nPyPionsRec=0;\r
+  TArrayD PzPionsRec(TrackNumber);       \r
+  Int_t nPzPionsRec=0;\r
+\r
+  //------ Trimomento He\r
+  TArrayD PxHeRec(TrackNumber);       \r
+  Int_t nPxHeRec=0;\r
+  TArrayD PyHeRec(TrackNumber);       \r
+  Int_t nPyHeRec=0;\r
+  TArrayD PzHeRec(TrackNumber);       \r
+  Int_t nPzHeRec=0;\r
+\r
+  Float_t impactXY=-999, impactZ=-999;\r
+  Float_t impactXYpi=-999, impactZpi=-999;\r
+  \r
+  Int_t PDGCodePos;\r
+  Int_t PDGCodeNeg;\r
+  Int_t motherPDGNeg;\r
+  Int_t motherPDGPos;\r
+\r
+  Int_t mumpdg=-100;\r
+  \r
+  //!   SELECTIONS:\r
+  //! - No ITSpureSA\r
+  //! - ITSrefit\r
+  //! - TPCrefit\r
+  //! - ITSmap !=0\r
+\r
+  // ******************* Track Cuts Definitions ********************//\r
+  \r
+  //! ITS\r
+\r
+  AliESDtrackCuts* esdtrackCutsITS = new AliESDtrackCuts("esdtrackCutsITS");\r
+  //  esdtrackCutsITS->SetRequireITSRefit(kTRUE);\r
+  esdtrackCutsITS->SetRequireITSStandAlone(kFALSE);\r
+  esdtrackCutsITS->SetRequireITSPureStandAlone(kFALSE);\r
+\r
+  //! TPC\r
+  \r
+  Int_t    minclsTPC=60;\r
+  Double_t maxchi2perTPCcl=5.;\r
+  \r
+  AliESDtrackCuts* esdtrackCutsTPC = new AliESDtrackCuts("esdtrackCutsTPC");\r
+  esdtrackCutsTPC->SetRequireTPCRefit(kTRUE);\r
+  esdtrackCutsTPC->SetAcceptKinkDaughters(kFALSE);\r
+  esdtrackCutsTPC->SetMinNClustersTPC(minclsTPC);\r
+  esdtrackCutsTPC->SetMaxChi2PerClusterTPC(maxchi2perTPCcl);\r
+        \r
+  //*************************************************************\r
+ \r
+\r
+  //--------------------------------------------\r
+  //cout<<"TrackNumber: "<<TrackNumber<<endl;\r
+\r
+  for (Int_t j=0; j<TrackNumber; j++) { //loop on tracks\r
+  \r
+    AliESDtrack *esdtrack=lESDevent->GetTrack(j);\r
+   \r
+    if(!esdtrack) { \r
+      AliError(Form("ERROR: Could not retrieve esdtrack %d",j)); \r
+      continue; \r
+    }\r
+\r
+    hBBTPCnoCuts->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());\r
+\r
+    // ************** Track cuts ****************\r
+    \r
+    status  = (ULong_t)esdtrack->GetStatus();\r
+    \r
+    isTPC   = ((status & AliESDtrack::kTPCin)  != 0);\r
+    isTOF   = (((status & AliESDtrack::kTOFout) != 0) && ((status & AliESDtrack::kTIME) != 0));\r
+    \r
+    Bool_t IsTrackAcceptedTPC =  esdtrackCutsTPC->AcceptTrack(esdtrack);\r
+    Bool_t IsTrackAcceptedITS =  esdtrackCutsITS->AcceptTrack(esdtrack);\r
+\r
+    if (!(IsTrackAcceptedTPC && IsTrackAcceptedITS)) continue;\r
+\r
+    //----------------------------------------------\r
+    \r
+    //****** Cuts from  AliV0Vertex.cxx *************\r
+     \r
+    Double_t d=esdtrack->GetD(xPrimaryVertex,yPrimaryVertex,lMagneticField);\r
+    //    if (TMath::Abs(d)<fgDPmin) continue;\r
+    if (TMath::Abs(d)>fgRmax) continue;\r
+    \r
+    //---- (Usefull) Stuff\r
+    \r
+    TPCSignal=esdtrack->GetTPCsignal(); \r
+    \r
+    if (TPCSignal<10)continue;\r
+    \r
+    if(!isTPC)continue;\r
+\r
+    if(!esdtrack->GetTPCInnerParam())continue;\r
+    \r
+    AliExternalTrackParam trackIn(*esdtrack->GetInnerParam()); \r
+    pinTPC = trackIn.GetP(); \r
+        \r
+    poutTPC=pinTPC;\r
+  \r
+    fhBBTPC->Fill(pinTPC*esdtrack->GetSign(),TPCSignal);\r
+    \r
+    d=esdtrack->GetD(xPrimaryVertex,yPrimaryVertex,lMagneticField);\r
+    //    if (TMath::Abs(d)<fgDPmin) continue;\r
+    if (TMath::Abs(d)>fgRmax) continue;\r
+\r
+    label = TMath::Abs(esdtrack->GetLabel());\r
+    \r
+    if (label>=10000000) {\r
+      // Underlying event. 10000000 is the\r
+      // value of fkMASKSTEP in AliRunDigitizer\r
+      cout <<"Strange, there should be no underlying event"<<endl;\r
+    }\r
+    \r
+    else {\r
+      if (label>=lNbMCPart) {\r
+       cout <<"Strange, label outside the range"<< endl;\r
+       continue;\r
+      }\r
+    }\r
+    \r
+    TParticle * part = stack->Particle(label);\r
+    //if(part)\r
+    //         part->Print();\r
+      \r
+    //      cout<<"PDG CODE: "<<part->GetPdgCode()<<endl;\r
+       \r
+    Int_t PDGCode=part->GetPdgCode();\r
+    Int_t mumid = part->GetFirstMother();\r
+\r
+    //      cout<<"Traccia ricostruita: \nID:"<<PDGCode<<" label "<<label<<" Theta Mum "<<thetaMumRec<<" mum id "<<mumid<<endl;      \r
+    if(mumid>-1){\r
+      TParticle *mother=(TParticle*)stack->Particle(mumid);\r
+      mumpdg = mother->GetPdgCode();\r
+    }\r
+\r
+    //     if (track->GetSign() < 0.){\r
+       \r
+    //    if(PDGCode==-211 || PDGCode==+211){\r
+    \r
+    if(PDGCode==-211)\r
+      fhBBTPCNegativePions->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());\r
+    \r
+    if(PDGCode==+211)\r
+      fhBBTPCPositivePions->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());\r
+\r
+\r
+    //    if(PDGCode == 211){\r
+         \r
+    if(PDGCode==-211 || PDGCode==+211){\r
+      //if(esdtrack->GetSign()>0)continue;\r
+\r
+      PionsTPC[nPionsTPC++]=j;\r
+    \r
+      esdtrack->GetImpactParameters(impactXY, impactZ);\r
+      fhPiDCAXY->Fill(impactXY);\r
+      fhPiDCAZ->Fill(impactZ);\r
+      \r
+      MomPionsRec[nPionsRec++]=esdtrack->P();\r
+\r
+      PxPionsRec[nPxPionsRec++]=esdtrack->Px();\r
+      PyPionsRec[nPyPionsRec++]=esdtrack->Py();\r
+      PzPionsRec[nPzPionsRec++]=esdtrack->Pz();\r
+       \r
+    }\r
+       \r
+    if(PDGCode==1000020030 ||PDGCode==-1000020030 ){\r
+\r
+    //    if(PDGCode==-1000020030){\r
+       \r
+      //cout<<"He3"<<endl;\r
+      //  if(esdtrack->GetSign()<0)continue;\r
+\r
+      HeTPC[nHeTPC++]=j;\r
+\r
+      fhBBTPCHe3->Fill(esdtrack->GetSign()*esdtrack->P(),esdtrack->GetTPCsignal());\r
+\r
+      esdtrack->GetImpactParameters(impactXY, impactZ);\r
+      fhHeDCAXY->Fill(impactXY);\r
+      fhHeDCAZ->Fill(impactZ);\r
+      \r
+      MomHeRec[nHeRec++]=esdtrack->P();\r
+\r
+      PxHeRec[nPxHeRec++]=esdtrack->Px();\r
+      PyHeRec[nPyHeRec++]=esdtrack->Py();\r
+      PzHeRec[nPzHeRec++]=esdtrack->Pz();\r
+       \r
+    }  \r
+      \r
+  }  //     ! track\r
+\r
+  //-------------- LOOP pairs 1 -------------\r
+  // Fill phase space and inva mass before cuts\r
+  \r
+  Double_t HeMomRec =0;\r
+  Double_t PionMomRec=0;\r
+  Double_t PxHeReco=0, PyHeReco=0,PzHeReco=0;\r
+  Double_t PxPionReco=0, PyPionReco=0,PzPionReco=0;\r
+\r
+  for(Int_t l=0; l < nHeRec; l++){\r
+      \r
+    HeMomRec=MomHeRec[l];\r
+\r
+    PxHeReco=PxHeRec[l];\r
+    PyHeReco=PyHeRec[l];\r
+    PzHeReco=PzHeRec[l];\r
+\r
+    for(Int_t k=0; k < nPionsRec; k++){\r
+   \r
+      PionMomRec=MomPionsRec[k];\r
+      \r
+      PxPionReco=PxPionsRec[k];\r
+      PyPionReco=PyPionsRec[k];\r
+      PzPionReco=PzPionsRec[k];\r
+      \r
+      fhRigidityHevsMomPiRec->Fill(HeMomRec,PionMomRec);\r
+\r
+      vHeliumRec.SetXYZM(2*PxHeReco,2*PyHeReco,2*PzHeReco,Helium3Mass); \r
+      vPionRec.SetXYZM(PxPionReco,PyPionReco,PzPionReco,PionMass);       \r
+      vSumRec=vHeliumRec+vPionRec;\r
+      \r
+      fhInvMassRec->Fill(vSumRec.M());\r
+    }\r
+    \r
+  } // fine loop phase space\r
+\r
+  //--------------- LOOP PAIRS ----------------------//\r
+  \r
+  //    cout<<"nPosTrackHe3 "<<nHeTPC<<endl;\r
+  \r
+  Double_t        DcaHeToPrimVertex=0;\r
+  Double_t        DcaPionToPrimVertex=0;\r
+  \r
+  impactXY=-999, impactZ=-999;\r
+  impactXYpi=-999, impactZpi=-999;\r
+  \r
+  // Track \r
+  \r
+  AliESDtrack  *PionTrack = 0x0;\r
+  AliESDtrack  *HeTrack = 0x0;\r
+  \r
+  // Vettori per il PxPyPz\r
+  \r
+  Double_t momPionVett[3];\r
+  for(Int_t i=0;i<3;i++)momPionVett[i]=0;\r
+    \r
+  Double_t momHeVett[3];\r
+  for(Int_t i=0;i<3;i++)momHeVett[i]=0;\r
+  \r
+  //At SV\r
+  \r
+  Double_t momPionVettAt[3];\r
+  for(Int_t i=0;i<3;i++)momPionVettAt[i]=0;\r
+  \r
+  Double_t momHeVettAt[3];\r
+  for(Int_t i=0;i<3;i++)momHeVettAt[i]=0;\r
+  \r
+  //    Double_t fPos[3]={0.,0.,0.};\r
+  \r
+  Bool_t   IsHeITSRefit,IsPiITSRefit ;\r
+  \r
+  //----------- Mio 2nd Vertex Finder\r
+    \r
+  for (Int_t k=0; k < nPionsTPC; k++) {                           //! Pions Loop\r
+    \r
+    DcaPionToPrimVertex=0.;\r
+    DcaHeToPrimVertex=0;\r
+    \r
+    Int_t PionIdx=PionsTPC[k];\r
+    \r
+    PionTrack=lESDevent->GetTrack(PionIdx);\r
+    \r
+    statusPi = (ULong_t)PionTrack->GetStatus();\r
+    IsPiITSRefit = (statusPi & AliESDtrack::kITSrefit); \r
+    \r
+    Int_t labelPi = TMath::Abs(PionTrack->GetLabel());\r
+    TParticle * partNeg = stack->Particle(labelPi);\r
+    PDGCodeNeg=partNeg->GetPdgCode();\r
+    \r
+    Int_t mumidNeg = partNeg->GetFirstMother();\r
+    if(mumidNeg>-1){\r
+      TParticle *motherNeg=(TParticle*)stack->Particle(mumidNeg);\r
+      motherPDGNeg = motherNeg->GetPdgCode();\r
+    }\r
+    \r
+    if (PionTrack) \r
+      DcaPionToPrimVertex = TMath::Abs(PionTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK\r
+\r
+ \r
+  \r
+    if(DcaPionToPrimVertex<0.1)continue; //qui\r
+    \r
+    AliExternalTrackParam trackInPion(*PionTrack);  \r
+    \r
+    for (Int_t i=0; i<nHeTPC; i++){                               //! Helium Loop\r
+      \r
+      Int_t HeIdx=HeTPC[i];\r
+      \r
+      HeTrack=lESDevent->GetTrack(HeIdx);\r
+      \r
+      statusT= (ULong_t)HeTrack->GetStatus();\r
+      IsHeITSRefit = (statusT & AliESDtrack::kITSrefit); \r
+      \r
+      Int_t labelHe = TMath::Abs(HeTrack->GetLabel());\r
+      TParticle * partPos = stack->Particle(labelHe);\r
+      PDGCodePos=partPos->GetPdgCode();\r
+      \r
+      Int_t mumidPos = partPos->GetFirstMother();\r
+      if(mumidPos>-1){\r
+       TParticle *motherPos=(TParticle*)stack->Particle(mumidPos);\r
+       motherPDGPos = motherPos->GetPdgCode();\r
+      }\r
+      \r
+      if (HeTrack) \r
+       DcaHeToPrimVertex = TMath::Abs(HeTrack->GetD(xPrimaryVertex, yPrimaryVertex,lMagneticField)); //OK\r
+      \r
+      AliExternalTrackParam trackInHe(*HeTrack); \r
+     \r
+      HeTrack->PxPyPz(momHeVett);\r
+      PionTrack->PxPyPz(momPionVett);   \r
+     \r
+      vHeliumRec.SetXYZM(2*momHeVett[0],2*momHeVett[1],2*momHeVett[2],Helium3Mass); \r
+      vPionRec.SetXYZM(momPionVett[0],momPionVett[1],momPionVett[2],PionMass);       \r
+      vSumRec=vHeliumRec+vPionRec;\r
+      \r
+      fhInvMassRec1->Fill(vSumRec.M());\r
+      \r
+      //      if(DcaPionToPrimVertex<0.1)continue; \r
+      fhInvMassRec2->Fill(vSumRec.M());\r
+      \r
+      if ( DcaPionToPrimVertex < fgDNmin)                //OK\r
+       if ( DcaHeToPrimVertex < fgDNmin) continue;    //OK\r
+\r
+      fhInvMassRec3->Fill(vSumRec.M());\r
+\r
+      Double_t xn, xp;\r
+      Double_t dca=0.;\r
+      \r
+      dca= PionTrack->GetDCA(HeTrack,lMagneticField,xn,xp); //!distanza tra le due tracce (Neg to Pos)\r
+      fHistProvaDCA->Fill(xn-xp,dca);\r
+      if (dca > fgDCAmax) continue;\r
+\r
+      fhInvMassRec4->Fill(vSumRec.M());\r
+       \r
+      if ((xn+xp) > 2*fgRmax) continue;\r
+      if ((xn+xp) < 2*fgRmin) continue;\r
+      fhInvMassRec5->Fill(vSumRec.M());\r
+\r
+     \r
+      \r
+      //CORREZIONE da AliV0Vertex\r
+      \r
+      Bool_t corrected=kFALSE;\r
+      if ((trackInPion.GetX() > 3.) && (xn < 3.)) {\r
+       //correct for the beam pipe material\r
+       corrected=kTRUE;\r
+      }\r
+      if ((trackInHe.GetX() > 3.) && (xp < 3.)) {\r
+       //correct for the beam pipe material\r
+       corrected=kTRUE;\r
+      }\r
+      if (corrected) {\r
+       dca=trackInPion.GetDCA(&trackInHe,lMagneticField,xn,xp);\r
+       if (dca > fgDCAmax) continue;\r
+       if ((xn+xp) > 2*fgRmax) continue;\r
+       if ((xn+xp) < 2*fgRmin) continue;\r
+      }\r
+   \r
+      //=============================================//\r
+      // Faccio un "V0" con le tracce che ho trovato //\r
+      //=============================================//\r
+      \r
+      trackInPion.PropagateTo(xn,lMagneticField); \r
+      trackInHe.PropagateTo(xp,lMagneticField);\r
+            \r
+      AliESDv0 vertex(trackInPion,PionIdx,trackInHe,HeIdx);\r
+      if (vertex.GetChi2V0() > fgChi2max) continue;\r
+      fhInvMassRec6->Fill(vSumRec.M());\r
+\r
+      Float_t CosPointingAngle=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex); //PointingAngle\r
+      if (CosPointingAngle < fgCPAmin) continue;\r
+    \r
+      fhInvMassRec7->Fill(vSumRec.M());\r
+\r
+      vertex.SetDcaV0Daughters(dca);\r
+      vertex.SetV0CosineOfPointingAngle(CosPointingAngle);\r
+\r
+      fPos[0]=vertex.Xv();\r
+      fPos[1]=vertex.Yv(); \r
+      fPos[2]=vertex.Zv(); \r
+      \r
+\r
+      \r
+      Double_t raggio=TMath::Sqrt(fPos[0]*fPos[0]+fPos[1]*fPos[1]+fPos[2]*fPos[2]);\r
+      HeTrack->GetPxPyPzAt(raggio,lMagneticField,momHeVettAt);\r
+      PionTrack->GetPxPyPzAt(raggio,lMagneticField,momPionVettAt); \r
+      \r
+      //------------------------------------------------------------------------//\r
+\r
+      HeTrack->GetImpactParameters(impactXY, impactZ);\r
+      \r
+      PionTrack->GetImpactParameters(impactXYpi, impactZpi);\r
+      \r
+      Float_t timeTOFHe= HeTrack->GetTOFsignal();                 // ps\r
+      Float_t trackLenghtTOFHe= HeTrack->GetIntegratedLength();  // cm\r
+\r
+      Float_t timeTOFPi= PionTrack->GetTOFsignal();                 // ps\r
+      Float_t trackLenghtTOFPi= PionTrack->GetIntegratedLength();  // cm\r
+\r
+      //----------------------------------------------------------------------//\r
+\r
+      vett1[0]=(Float_t)runNumber;\r
+      vett1[1]=(Float_t)evNumber;\r
+      vett1[2]=(Float_t)lNbMCPart;\r
+      vett1[3]=(Float_t)percentile;\r
+      vett1[4]=(Float_t)xPrimaryVertex; //PRIMARY\r
+      vett1[5]=(Float_t)yPrimaryVertex;\r
+      vett1[6]=(Float_t)zPrimaryVertex;\r
+      vett1[7]=(Float_t)fPos[0]; //SECONDARY\r
+      vett1[8]=(Float_t)fPos[1];\r
+      vett1[9]=(Float_t)fPos[2];\r
+      vett1[10]=(Float_t)dca;           //between 2 tracks\r
+      vett1[11]=(Float_t)CosPointingAngle;          //cosPointingAngle da V0\r
+      vett1[12]=(Float_t)vertex.GetD(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);\r
+      vett1[13]=(Float_t)HeTrack->GetSign(); //He\r
+      vett1[14]=(Float_t)trackInHe.GetP();\r
+      vett1[15]=(Float_t)HeTrack->GetTPCsignal();\r
+      vett1[16]=(Float_t)DcaHeToPrimVertex;\r
+      vett1[17]=(Float_t)HeTrack->Eta();\r
+      vett1[18]=(Float_t)momHeVett[0];\r
+      vett1[19]=(Float_t)momHeVett[1];\r
+      vett1[20]=(Float_t)momHeVett[2];\r
+      vett1[21]=(Float_t)momHeVettAt[0];\r
+      vett1[22]=(Float_t)momHeVettAt[1];\r
+      vett1[23]=(Float_t)momHeVettAt[2];\r
+      vett1[24]=(Float_t)HeTrack->GetTPCNcls();\r
+      vett1[25]=(Float_t)impactXY;\r
+      vett1[26]=(Float_t)impactZ;\r
+      vett1[27]=(Float_t)isTOFHe3;\r
+      vett1[28]=(Float_t)(trackLenghtTOFHe/timeTOFHe)/2.99792458e-2;\r
+      vett1[29]=(Float_t)HeTrack->GetITSClusterMap();\r
+      vett1[30]=(Float_t)IsHeITSRefit;\r
+      vett1[31]=(Float_t)PionTrack->GetSign(); //Pion\r
+      vett1[32]=(Float_t)trackInPion.GetP();\r
+      vett1[33]=(Float_t)PionTrack->GetTPCsignal();\r
+      vett1[34]=(Float_t)DcaPionToPrimVertex;\r
+      vett1[35]=(Float_t)PionTrack->Eta();\r
+      vett1[36]=(Float_t)momPionVett[0];\r
+      vett1[37]=(Float_t)momPionVett[1];\r
+      vett1[38]=(Float_t)momPionVett[2];\r
+      vett1[39]=(Float_t)momPionVettAt[0];\r
+      vett1[40]=(Float_t)momPionVettAt[1];\r
+      vett1[41]=(Float_t)momPionVettAt[2];\r
+      vett1[42]=(Float_t)PionTrack->GetTPCNcls();\r
+      vett1[43]=(Float_t)impactXYpi;\r
+      vett1[44]=(Float_t)impactZpi;\r
+      vett1[45]=(Float_t)isTOFPi;\r
+      vett1[46]=(Float_t)(trackLenghtTOFPi/timeTOFPi)/2.99792458e-2;\r
+      vett1[47]=(Float_t)PionTrack->GetITSClusterMap();\r
+      vett1[48]=(Float_t)IsPiITSRefit;\r
+      vett1[49]=(Float_t)PDGCodeNeg;\r
+      vett1[50]=(Float_t)PDGCodePos;\r
+      vett1[51]=(Float_t)motherPDGNeg;\r
+      vett1[52]=(Float_t)motherPDGPos;\r
+      vett1[53]=(Float_t)labelPi;\r
+      vett1[54]=(Float_t)labelHe;\r
+      vett1[55]=(Float_t)mumidNeg;\r
+      vett1[56]=(Float_t)mumidPos;\r
+\r
+      fNtuple1->Fill(vett1);  \r
+\r
+    }// positive TPC\r
+    \r
+  } //negative tpc\r
+\r
+  PostData(1,fListHistCascade);\r
+  \r
+}// end userexec\r
+\r
+\r
+//________________________________________________________________________\r
+//\r
+void AliAnalysisTaskHelium3PiMC::Terminate(Option_t *) \r
+{\r
+  //  Draw result to the screen\r
+  //   Called once at the end of the query\r
+}\r
+\r

diff --git a/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3PiMC.h b/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3PiMC.h
new file mode 100644 (file)
index 0000000..da41f40
--- /dev/null
+++ b/PWGLF/STRANGENESS/Hypernuclei/AliAnalysisTaskHelium3PiMC.h
@@ -0,0 +1,113 @@
+#ifndef ALIANALYSISTASKHELIUM3PIMC_H
+#define ALIANALYSISTASKHELIUM3PIMC_H
+
+/*  See cxx source for full Copyright notice */
+
+//-----------------------------------------------------------------
+//                 AliAnalysisTaskHelium3Pion class
+//-----------------------------------------------------------------
+
+class TList;
+class TH1F;
+class TH2F;
+class TH3F;
+class TNtuple;
+class AliESDcascade;
+//class AliCascadeVertexer; 
+
+#include "TString.h"
+
+#include "AliAnalysisTaskSE.h"
+
+class AliAnalysisTaskHelium3PiMC : public AliAnalysisTaskSE {
+ public:
+  AliAnalysisTaskHelium3PiMC();
+  AliAnalysisTaskHelium3PiMC(const char *name);
+  virtual ~AliAnalysisTaskHelium3PiMC();
+  
+  virtual void  UserCreateOutputObjects();
+  virtual void  UserExec(Option_t *option);
+  virtual void  Terminate(Option_t *);
+  
+  void SetCollidingSystems(Short_t collidingSystems = 0)     {fCollidingSystems = collidingSystems;}
+  void SetAnalysisType    (const char* analysisType = "ESD") {fAnalysisType = analysisType;}
+  void SetDataType    (const char* dataType = "SIM") {fDataType = dataType;}
+  void SetMinNSigma(Double_t ns=3.){
+    fMinNSigma=ns;
+  }
+
+  void SetMinSPDPoints(Int_t np=1){
+    fMinSPDPts=np;
+  }
+
+  //Double_t BetheBloch(Double_t bg,Double_t Charge,Bool_t optMC);
+  Double_t BetheBloch(Double_t bg,Double_t Charge,Bool_t isPbPb);
+
+  
+ private:
+
+  TString fAnalysisType;            //! "ESD" or "AOD" analysis type   
+
+  Short_t fCollidingSystems;        //! 0 = pp collisions or 1 = AA collisions
+  TString fDataType;                //! "REAL" or "SIM" data type      
+  TList        *fListHistCascade;           //! List of Cascade histograms
+
+  TH1F *fHistEventMultiplicity;      //! event multiplicity
+  TH1F *fHistTrackMultiplicity;      //! track multiplicity  
+  TH1F *fHistMCMultiplicityTracks;
+  TH1F *fHistMCEta; 
+  TH1F *fHistMCPt; 
+  TH1F *fHistMCTheta; 
+  TH1F *fHistMCDecayPosition;
+  TH1F *fHistMCDecayRho; 
+  TH2F *fhRigidityHevsMomPiMC;
+  TH2F *fhRigidityHevsMomPiRec;
+  TH1F *fhInvMassMC;
+  TH1F *fhInvMassMum;
+  TH1F *fhInvMassRec;
+
+  TH1F *fhInvMassRec1;
+  TH1F *fhInvMassRec2;
+  TH1F *fhInvMassRec3;
+  TH1F *fhInvMassRec4;
+  TH1F *fhInvMassRec5;
+  TH1F *fhInvMassRec6;
+  TH1F *fhInvMassRec7;
+
+  TH1F *fhHeMCRigidity;
+  TH1F *fhPioneMC;
+  TH2F *hBBTPCnoCuts;
+  TH2F *fhBBTPC;
+  TH2F *fhBBTPCNegativePions;
+  TH2F *fhBBTPCPositivePions;
+  TH2F *fhBBTPCHe3;
+  TH2F *fHistProvaDCA;
+  TH2F *fHistPercentileVsTrackNumber;
+  TH1F *fhHeDCAXY; 
+  TH1F *fhHeDCAZ; 
+  TH1F *fhPiDCAXY; 
+  TH1F *fhPiDCAZ; 
+  TH1F *hITSClusterMap;
+  
+ 
+  TNtuple *fNtuple1;                  //! Ntupla Pairs Pi/Proton "standard"
+  TNtuple *fNtuple2;                  //! Ntupla Pairs PiPos/Proton "background"
+
+  TTree *fhTreeESD;                  //! Tree di Tstring
+  TString nameESD ; 
+  Int_t nEv;
+
+  Int_t    fMinSPDPts;       // minimum number of SPD Points
+  Double_t fMinNSigma;       // minimum number of sigmas
+  
+
+ static const Int_t fgNrot;
+ 
+
+  AliAnalysisTaskHelium3PiMC(const AliAnalysisTaskHelium3PiMC&);            // not implemented
+  AliAnalysisTaskHelium3PiMC& operator=(const AliAnalysisTaskHelium3PiMC&); // not implemented
+  
+  ClassDef(AliAnalysisTaskHelium3PiMC, 0);
+};
+
+#endif