--- /dev/null
+/**************************************************************************\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
--- /dev/null
+/**************************************************************************\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