]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Adding the feed-down task/class for protons (Marek)
authorpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 12 Aug 2009 10:44:24 +0000 (10:44 +0000)
committerpchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 12 Aug 2009 10:44:24 +0000 (10:44 +0000)
PWG2/SPECTRA/AliProtonFeedDownAnalysis.cxx [new file with mode: 0644]
PWG2/SPECTRA/AliProtonFeedDownAnalysis.h [new file with mode: 0644]
PWG2/SPECTRA/AliProtonFeedDownAnalysisTask.cxx [new file with mode: 0644]
PWG2/SPECTRA/AliProtonFeedDownAnalysisTask.h [new file with mode: 0644]
PWG2/SPECTRA/macros/AddTaskProtonFeedDownAnalysis.C [new file with mode: 0644]
PWG2/SPECTRA/macros/configProtonFeedDownAnalysis.C [new file with mode: 0644]

diff --git a/PWG2/SPECTRA/AliProtonFeedDownAnalysis.cxx b/PWG2/SPECTRA/AliProtonFeedDownAnalysis.cxx
new file mode 100644 (file)
index 0000000..0cecdf9
--- /dev/null
@@ -0,0 +1,382 @@
+#include <Riostream.h>\r
+#include <TFile.h>\r
+#include <TSystem.h>\r
+#include <TF1.h>\r
+#include <TH2D.h>\r
+#include <TH1D.h>\r
+#include <TH1I.h>\r
+#include <TParticle.h>\r
+\r
+#include <AliExternalTrackParam.h>\r
+#include <AliAODEvent.h>\r
+#include <AliESDEvent.h>\r
+//#include <AliLog.h>\r
+#include <AliPID.h>\r
+#include <AliStack.h>\r
+#include <AliCFContainer.h>\r
+#include <AliCFEffGrid.h>\r
+#include <AliCFDataGrid.h>\r
+//#include <AliESDVertex.h>\r
+class AliLog;\r
+class AliESDVertex;\r
+\r
+#include "AliProtonFeedDownAnalysis.h"\r
+#include "AliProtonAnalysisBase.h"\r
+\r
+ClassImp(AliProtonFeedDownAnalysis)\r
+\r
+//____________________________________________________________________//\r
+AliProtonFeedDownAnalysis::AliProtonFeedDownAnalysis() : \r
+  TObject(), fProtonAnalysisBase(0),\r
+  fNBinsY(0), fMinY(0), fMaxY(0),\r
+  fNBinsPt(0), fMinPt(0), fMaxPt(0),\r
+  fProtonContainer(0), fAntiProtonContainer(0), fweightfunction(0),fLambda(0),fLambdaweighted(0),fAntiLambda(0),fAntiLambdaweighted(0)\r
+  {\r
+  //Default constructor\r
+ }\r
+//____________________________________________________________________//\r
+AliProtonFeedDownAnalysis::AliProtonFeedDownAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) : \r
+  TObject(), fProtonAnalysisBase(0),\r
+  fNBinsY(nbinsY), fMinY(fLowY), fMaxY(fHighY),\r
+  fNBinsPt(nbinsPt), fMinPt(fLowPt), fMaxPt(fHighPt),\r
+  fProtonContainer(0), fAntiProtonContainer(0),fweightfunction(0),fLambda(0),fLambdaweighted(0),fAntiLambda(0),fAntiLambdaweighted(0)\r
+  {\r
+       //Default constructor\r
+       \r
+       //setting up the containers\r
+       Int_t iBin[2];\r
+       iBin[0] = nbinsY;\r
+       iBin[1] = nbinsPt;\r
+       Double_t *binLimY = new Double_t[nbinsY+1];\r
+       Double_t *binLimPt = new Double_t[nbinsPt+1];\r
+       //values for bin lower bounds\r
+       for(Int_t i = 0; i <= nbinsY; i++) \r
+               binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;\r
+       for(Int_t i = 0; i <= nbinsPt; i++) \r
+               binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;\r
+       \r
+       fProtonContainer = new AliCFContainer("containerProtons","container for protons",4,2,iBin);\r
+       fProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta\r
+       fProtonContainer->SetBinLimits(1,binLimPt); //pT\r
+       fAntiProtonContainer = new AliCFContainer("containerAntiProtons","container for antiprotons",4,2,iBin);\r
+       fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity or eta\r
+       fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT\r
+       \r
+}\r
+//____________________________________________________________________//\r
+AliProtonFeedDownAnalysis::~AliProtonFeedDownAnalysis() \r
+{\r
+       //Default destructor\r
+       if(fProtonAnalysisBase) delete fProtonAnalysisBase;\r
+       if(fProtonContainer) delete fProtonContainer;\r
+       if(fAntiProtonContainer) delete fAntiProtonContainer;\r
+       if(fweightfunction) delete fweightfunction; \r
+       if(fLambda) delete fLambda;\r
+       if(fLambdaweighted) delete fLambdaweighted;\r
+       if(fAntiLambda) delete fAntiLambda;\r
+       if(fAntiLambdaweighted) delete fAntiLambdaweighted;\r
+}\r
+//____________________________________________________________________//\r
+void AliProtonFeedDownAnalysis::InitAnalysisHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt) \r
+{\r
+       //Initializes the histograms\r
+       fNBinsY = nbinsY;\r
+       fMinY = fLowY;\r
+       fMaxY = fHighY;\r
+       fNBinsPt = nbinsPt;\r
+       fMinPt = fLowPt;\r
+       fMaxPt = fHighPt;\r
+       \r
+       \r
+       //setting up the containers\r
+       Int_t iBin[2];\r
+       iBin[0] = nbinsY;\r
+       iBin[1] = nbinsPt;\r
+       Double_t *binLimY = new Double_t[nbinsY+1];\r
+       Double_t *binLimPt = new Double_t[nbinsPt+1];\r
+       //values for bin lower bounds\r
+       for(Int_t i = 0; i <= nbinsY; i++) \r
+               binLimY[i]=(Double_t)fLowY  + (fHighY - fLowY)  /nbinsY*(Double_t)i;\r
+       for(Int_t i = 0; i <= nbinsPt; i++) \r
+               binLimPt[i]=(Double_t)fLowPt  + (fHighPt - fLowPt)  /nbinsPt*(Double_t)i;\r
+       \r
+       fProtonContainer = new AliCFContainer("containerProtons","container for protons",kNSteps,2,iBin);\r
+       fProtonContainer->SetBinLimits(0,binLimY); //rapidity\r
+       fProtonContainer->SetBinLimits(1,binLimPt); //pT\r
+       fAntiProtonContainer = new AliCFContainer("containerAntiProtons","container for antiprotons",kNSteps,2,iBin);\r
+       fAntiProtonContainer->SetBinLimits(0,binLimY); //rapidity\r
+       fAntiProtonContainer->SetBinLimits(1,binLimPt); //pT\r
+       fLambda=new TH1F("Lambda","Lambda",35,0.5,4.0);\r
+       fLambdaweighted=new TH1F("Lambdaweighted","Lambdaweighted",35,0.5,4.0);\r
+       fAntiLambda=new TH1F("AntiLambda","AntiLambda",35,0.5,4.0);\r
+       fAntiLambdaweighted=new TH1F("AntiLambdaweighted","AntiLambdaweighted",35,0.5,4.0);\r
+}\r
+//_________________________________________________________________________//\r
+void AliProtonFeedDownAnalysis::Analyze(AliESDEvent *esd, const AliESDVertex *vertex,AliStack *stack)\r
+{      \r
+       Int_t nTracks = 0;\r
+       Int_t nIdentifiedProtons = 0, nIdentifiedAntiProtons = 0;\r
+       Int_t nSurvivedProtons = 0, nSurvivedAntiProtons = 0;\r
+       \r
+       //fHistEvents->Fill(0); //number of analyzed events\r
+       Double_t containerInput[2] ;\r
+       Double_t gPt = 0.0, gP = 0.0;\r
+       nTracks = esd->GetNumberOfTracks();\r
+       Float_t weight;\r
+       Int_t trackflag;\r
+       for(Int_t iTracks = 0; iTracks < nTracks; iTracks++) \r
+       {\r
+               \r
+               AliESDtrack* track = esd->GetTrack(iTracks);\r
+               AliESDtrack trackTPC;\r
+               Int_t label= track->GetLabel();\r
+       /*      Int_t trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else \r
+               if (trackflag!=0)\r
+                       weight=GetWeightforProton(label,stack); \r
+               else\r
+                       weight=1.0;     */\r
+       \r
+               if((fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kTPC)||(fProtonAnalysisBase->GetAnalysisMode()==AliProtonAnalysisBase::kHybrid)) \r
+               {\r
+                       AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)track->GetTPCInnerParam();\r
+                       if(!tpcTrack)\r
+                                continue;\r
+                       gPt = tpcTrack->Pt();\r
+                       gP = tpcTrack->P();\r
+                       if(fProtonAnalysisBase->IsProton(track)) \r
+                       {\r
+                               trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else \r
+                               if (trackflag!=0)\r
+                                       weight=GetWeightforProton(label,stack); \r
+                               else\r
+                                       weight=1.0;     \r
+                               if(tpcTrack->Charge() > 0) \r
+                               {\r
+                                       nIdentifiedProtons += 1;\r
+                                       if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) \r
+                                               continue;//track cuts\r
+                                       if(!fProtonAnalysisBase->IsInPhaseSpace(track)) \r
+                                               continue; //track outside the analyzed y-Pt\r
+                                       nSurvivedProtons += 1;\r
+                                       if(fProtonAnalysisBase->GetEtaMode()) \r
+                                       {\r
+                                               containerInput[0] = tpcTrack->Eta();\r
+                                       }\r
+                                       else \r
+                                       {\r
+                                               //fill the container\r
+                                               containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());\r
+                                       }\r
+                                       containerInput[1] = gPt;\r
+                                       fProtonContainer->Fill(containerInput,kSelected,weight);   \r
+                                       //Feed-down check\r
+                                       if (trackflag==1)\r
+                                               fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight); \r
+                                       \r
+                               }//protons\r
+                               else if(tpcTrack->Charge() < 0) \r
+                               {\r
+                                       nIdentifiedAntiProtons += 1;\r
+                                       if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) \r
+                                               continue;//track cuts\r
+                                       if(!fProtonAnalysisBase->IsInPhaseSpace(track)) \r
+                                               continue; //track outside the analyzed y-Pt\r
+                                       nSurvivedAntiProtons += 1;\r
+                                       if(fProtonAnalysisBase->GetEtaMode()) \r
+                                       {\r
+                                               containerInput[0] = tpcTrack->Eta();\r
+                                       }\r
+                                       else \r
+                                       {\r
+                                               //fill the container\r
+                                               containerInput[0] = fProtonAnalysisBase->Rapidity(tpcTrack->Px(),tpcTrack->Py(),tpcTrack->Pz());\r
+                                       }\r
+                                       containerInput[1] = gPt;\r
+                                       fAntiProtonContainer->Fill(containerInput,kSelected,weight);\r
+                                       //Feed-down check                                                       \r
+                                       if(trackflag==-1)\r
+                                               fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);                                  \r
+                               }//antiprotons   \r
+                       }//proton check\r
+               }//TPC only tracks\r
+               else if(fProtonAnalysisBase->GetAnalysisMode() == AliProtonAnalysisBase::kGlobal) \r
+               {\r
+                       gPt = track->Pt();\r
+                       gP = track->P();\r
+                       if(fProtonAnalysisBase->IsProton(track)) \r
+                       {\r
+                               trackflag=LambdaIsMother(label,stack);//1 mother lambda -1 mother anti lambda 0 mother something else \r
+                               if (trackflag!=0)\r
+                                       weight=GetWeightforProton(label,stack); \r
+                               else\r
+                                       weight=1.0;     \r
+                               if(track->Charge() > 0) \r
+                               {\r
+                                       nIdentifiedProtons += 1;\r
+                                       if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) \r
+                                               continue;//track cuts\r
+                                       if(!fProtonAnalysisBase->IsInPhaseSpace(track)) \r
+                                               continue; //track outside the analyzed y-Pt\r
+                                       nSurvivedProtons += 1;\r
+                                       if(fProtonAnalysisBase->GetEtaMode()) \r
+                                       {\r
+                                               containerInput[0] = track->Eta();\r
+                                       }\r
+                                       else \r
+                                       {\r
+                                               //fill the container\r
+                                               containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());\r
+                                       }\r
+                                       containerInput[1] = gPt;\r
+                                       fProtonContainer->Fill(containerInput,kSelected,weight);  \r
+                                        //Feed-down check\r
+                                       if (trackflag==1)\r
+                                               fProtonContainer->Fill(containerInput,kSelectedfromLambda,weight); \r
+                               }//protons\r
+                               else if(track->Charge() < 0) \r
+                               {\r
+                                       nIdentifiedAntiProtons += 1;\r
+                                       if(!fProtonAnalysisBase->IsAccepted(esd,vertex,track)) \r
+                                               continue;//track cuts\r
+                                       if(!fProtonAnalysisBase->IsInPhaseSpace(track))\r
+                                               continue; //track outside the analyzed y-Pt\r
+                                       nSurvivedAntiProtons += 1;\r
+                                       if(fProtonAnalysisBase->GetEtaMode()) \r
+                                       {\r
+                                               containerInput[0] = track->Eta();\r
+                                       }\r
+                                       else \r
+                                       {\r
+                                               //fill the container\r
+                                               containerInput[0] = fProtonAnalysisBase->Rapidity(track->Px(),track->Py(),track->Pz());\r
+                                       }\r
+                                       containerInput[1] = gPt;\r
+                                       fAntiProtonContainer->Fill(containerInput,kSelected,weight);   \r
+                                       if(trackflag==-1)\r
+                                               fAntiProtonContainer->Fill(containerInput,kSelectedfromLambda,weight);\r
+                               }//antiprotons\r
+                       }//proton check \r
+               }//combined tracking\r
+       }//track loop \r
+       \r
+       if(fProtonAnalysisBase->GetDebugMode())\r
+               Printf("Initial number of tracks: %d | Identified (anti)protons: %d - %d | Survived (anti)protons: %d - %d",nTracks,nIdentifiedProtons,nIdentifiedAntiProtons,nSurvivedProtons,nSurvivedAntiProtons);\r
+\r
+}\r
+//_________________________________________________________________________//\r
+void AliProtonFeedDownAnalysis::Analyze(AliAODEvent *fAOD)\r
+{\r
+}\r
+//_________________________________________________________________________//\r
+void AliProtonFeedDownAnalysis::Analyze(AliStack *stack)\r
+{\r
+        Double_t containerInput[2] ;\r
+        Float_t weight;\r
+       for (Int_t ipart=0; ipart<stack->GetNtrack(); ipart++) \r
+       { \r
+               TParticle *particle  = stack->Particle(ipart);\r
+               Int_t code=particle->GetPdgCode();\r
+               if (code==3122)\r
+               {\r
+                       fLambda->Fill(particle->Pt());\r
+                       fLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));                       \r
+               }\r
+               if (code==-3122)\r
+               {\r
+                       fAntiLambda->Fill(particle->Pt());\r
+                       fAntiLambdaweighted->Fill(particle->Pt(),GetWeightforLambda(particle->Pt()));                   \r
+               }\r
+               if (TMath::Abs(code)==2212)\r
+               {\r
+                       Int_t trackflag=LambdaIsMother(ipart,stack);//1 mother lambda -1 mother anti lambda 0 mother something else \r
+                       if (trackflag!=0)\r
+                               weight=GetWeightforProton(ipart,stack);\r
+                       else\r
+                               weight=1.0;     \r
+                       if(particle->GetUniqueID() == 13) //recject hadronic inter.\r
+                               continue; \r
+                       if((particle->Pt() > fMaxPt)||(particle->Pt() < fMinPt)) \r
+                               continue;\r
+                       if(fProtonAnalysisBase->GetEtaMode()) \r
+                       {\r
+                               if((particle->Eta()> fMaxY)||(particle->Eta()< fMinY))\r
+                                       continue; \r
+                       }       \r
+                       else\r
+                       {       \r
+                               if((fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz()) > fMaxY)||(fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz())< fMinY)) \r
+                                       continue;\r
+                       }               \r
+                       if(fProtonAnalysisBase->GetEtaMode()) \r
+                       {\r
+                               containerInput[0] =particle->Eta();\r
+                       }\r
+                       else \r
+                       {\r
+                               containerInput[0] = fProtonAnalysisBase->Rapidity(particle->Px(),particle->Py(),particle->Pz());\r
+                       }\r
+                       //containerInput[0] = particle->Eta() ;\r
+                       containerInput[1] = (Float_t)particle->Pt();\r
+                       if (particle->GetPdgCode()==2212)\r
+                       {\r
+                               fProtonContainer->Fill(containerInput,kAll,weight);\r
+                               if(ipart<stack->GetNprimary())\r
+                               {\r
+                                       fProtonContainer->Fill(containerInput,kPrimary,weight);\r
+                               }\r
+                               else\r
+                               {\r
+                                       if(trackflag==1)\r
+                                               fProtonContainer->Fill(containerInput,kFromLambda,weight);\r
+                               }\r
+                       }\r
+                       if (particle->GetPdgCode()==-2212)\r
+                       {\r
+                               fAntiProtonContainer->Fill(containerInput,kAll,weight);\r
+                               if(ipart<stack->GetNprimary())\r
+                               {\r
+                                       fAntiProtonContainer->Fill(containerInput,kPrimary,weight);\r
+                               }\r
+                               else\r
+                               {                                                               \r
+                                       if(trackflag==-1)\r
+                                               fAntiProtonContainer->Fill(containerInput,kFromLambda,weight);\r
+                               }                       \r
+                       }\r
+               }\r
+       }   \r
+}\r
+//______________________________________________________________________________________________\r
+Int_t AliProtonFeedDownAnalysis::LambdaIsMother(Int_t number, AliStack *stack)\r
+{\r
+        if(number<0)\r
+               return 0;\r
+       TParticle *particle  = stack->Particle(number);\r
+       Int_t motherPDG=0;\r
+       Int_t lmother=-1;\r
+       lmother=particle->GetFirstMother();\r
+       if (lmother<0)          \r
+               return 0;\r
+       TParticle *mparticle=stack->Particle(lmother);\r
+       motherPDG=mparticle->GetPdgCode();                                                                              \r
+       if(motherPDG==3122)\r
+               return 1;\r
+       if(motherPDG==-3122)    \r
+               return -1;\r
+       return 0;       \r
+} \r
+//___________________________________________________________________________________________\r
+Float_t AliProtonFeedDownAnalysis::GetWeightforProton(Int_t number,AliStack *stack)\r
+{\r
+        if(number<0)\r
+               return 1.0;\r
+       TParticle *particle  = stack->Particle(number);\r
+       Int_t lmother=particle->GetFirstMother();\r
+       TParticle *mparticle=stack->Particle(lmother);\r
+       return  fweightfunction->Eval(mparticle->Pt());\r
+}\r
+//______________________________________________________________________________________________\r
+Float_t AliProtonFeedDownAnalysis::GetWeightforLambda(Float_t pt)\r
+{\r
+       return  fweightfunction->Eval(pt);\r
+}\r
diff --git a/PWG2/SPECTRA/AliProtonFeedDownAnalysis.h b/PWG2/SPECTRA/AliProtonFeedDownAnalysis.h
new file mode 100644 (file)
index 0000000..85fdebe
--- /dev/null
@@ -0,0 +1,86 @@
+#ifndef ALIPROTONFEEDDOWNANALYSIS_H\r
+#define ALIPROTONFEEDDOWNANALYSIS_H\r
+\r
+#include "TObject.h"\r
+#include "TH1I.h"\r
+#include "AliCFContainer.h"\r
+class TF1;\r
+class TH2D;\r
+class TH1F;\r
+class TList;\r
+\r
+//#include "AliPID.h"\r
+class AliPID;\r
+\r
+class AliCFDataGrid;\r
+class AliAODEvent;\r
+class AliAODtrack;\r
+class AliESDEvent;\r
+class AliESDtrack;\r
+class AliExternalTrackParam;\r
+class AliStack;\r
+class AliESDVertex;\r
+class AliProtonAnalysisBase;\r
+\r
+class AliProtonFeedDownAnalysis : public TObject \r
+{\r
+       public:\r
+       enum \r
+               {\r
+                       kAll      = 0, //without protons reject the secondaries from hadronic inter. \r
+                       kPrimary = 1,\r
+                       kFromLambda        = 2,\r
+                       kSelected = 3 ,\r
+                       kSelectedfromLambda = 4,\r
+                       kNSteps=5\r
+               };\r
+               AliProtonFeedDownAnalysis();\r
+               AliProtonFeedDownAnalysis(Int_t nbinsY, Float_t fLowY, Float_t fHighY,Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt);\r
+               virtual ~AliProtonFeedDownAnalysis();\r
+               \r
+               void SetBaseAnalysis(AliProtonAnalysisBase * const baseAnalysis) {fProtonAnalysisBase = baseAnalysis;}\r
+               AliProtonAnalysisBase *GetProtonAnalysisBaseObject() const {return fProtonAnalysisBase;}\r
+               \r
+               void InitAnalysisHistograms(Int_t nbinsY, Float_t fLowY, Float_t fHighY, Int_t nbinsPt, Float_t fLowPt, Float_t fHighPt);\r
+               void Analyze(AliESDEvent *fESD, const AliESDVertex *vertex,AliStack *stack);\r
+               void Analyze(AliAODEvent *fAOD);\r
+               void Analyze(AliStack *stack);\r
+               \r
+               AliCFContainer *GetProtonContainer() const {return fProtonContainer;}\r
+               AliCFContainer *GetAntiProtonContainer() const {return fAntiProtonContainer;}\r
+               \r
+               void SetWeightFunction(TF1* weightfunction){fweightfunction=weightfunction;}\r
+               TF1* GetWeightFunction() const{return fweightfunction;} \r
+               TH1F* GetLambdaHist() const{return fLambda;}\r
+               TH1F* GetLambdaweightedHist() const{return fLambdaweighted;}\r
+               TH1F* GetAntiLambdaHist() const{return fAntiLambda;}\r
+               TH1F* GetAntiLambdaweightedHist() const{return fAntiLambdaweighted;}\r
+       private:\r
+               AliProtonFeedDownAnalysis(const AliProtonFeedDownAnalysis&); // Not implemented\r
+               AliProtonFeedDownAnalysis& operator=(const AliProtonFeedDownAnalysis&); // Not implemented\r
+               \r
+               AliProtonAnalysisBase *fProtonAnalysisBase;//base analysis object\r
+               Int_t LambdaIsMother(Int_t numbe, AliStack *stack); \r
+               Float_t GetWeightforProton(Int_t numbe, AliStack *stack);\r
+               Float_t GetWeightforLambda(Float_t pt);\r
+               \r
+               Int_t fNBinsY; //number of bins in y or eta\r
+               Double_t fMinY, fMaxY; //min & max value of y or eta\r
+               Int_t fNBinsPt;  //number of bins in pT\r
+               Double_t fMinPt, fMaxPt; //min & max value of pT\r
+               \r
+               //Analysis containers\r
+               AliCFContainer *fProtonContainer; //container for protons\r
+               AliCFContainer *fAntiProtonContainer; //container for antiprotons\r
+               \r
+               TF1*fweightfunction; \r
+               TH1F *fLambda;\r
+               TH1F *fLambdaweighted;\r
+               TH1F *fAntiLambda;\r
+               TH1F *fAntiLambdaweighted;\r
+               \r
+  ClassDef(AliProtonFeedDownAnalysis,1);\r
+};\r
+\r
+#endif\r
+\r
diff --git a/PWG2/SPECTRA/AliProtonFeedDownAnalysisTask.cxx b/PWG2/SPECTRA/AliProtonFeedDownAnalysisTask.cxx
new file mode 100644 (file)
index 0000000..3fae5d5
--- /dev/null
@@ -0,0 +1,171 @@
+#include "TChain.h"\r
+#include "TTree.h"\r
+#include "TString.h"\r
+#include "TList.h"\r
+#include "TH2F.h"\r
+\r
+#include "AliAnalysisTask.h"\r
+#include "AliAnalysisManager.h"\r
+\r
+#include "AliESDEvent.h"\r
+#include "AliESDInputHandler.h"\r
+#include "AliAODEvent.h"\r
+#include "AliAODInputHandler.h"\r
+#include "AliMCEventHandler.h"\r
+#include "AliMCEvent.h"\r
+#include "AliStack.h"\r
+#include "AliESDVertex.h"\r
+\r
+#include "AliProtonFeedDownAnalysis.h"\r
+#include "AliProtonAnalysisBase.h"\r
+#include "AliProtonFeedDownAnalysisTask.h"\r
+#include <Riostream.h>\r
+\r
+ClassImp(AliProtonFeedDownAnalysisTask)\r
+  \r
+//________________________________________________________________________ \r
+AliProtonFeedDownAnalysisTask::AliProtonFeedDownAnalysisTask()\r
+  : AliAnalysisTask(), fESD(0), fAOD(0), fMC(0),\r
+    fList(0), fProtonAnalysis(0),fStatHist(0)\r
+ {\r
+  //Dummy constructor\r
+  \r
+}\r
+\r
+//________________________________________________________________________\r
+AliProtonFeedDownAnalysisTask::AliProtonFeedDownAnalysisTask(const char *name) \r
+  : AliAnalysisTask(name, ""), fESD(0), fAOD(0), fMC(0),\r
+    fList(0), fProtonAnalysis(0),fStatHist(0)\r
+    {\r
+       // Constructor\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\r
+       DefineOutput(0, TList::Class());\r
+}\r
+\r
+//________________________________________________________________________\r
+void AliProtonFeedDownAnalysisTask::ConnectInputData(Option_t *) \r
+{\r
+       // Connect ESD or AOD here\r
+       // Called once\r
+       TString gAnalysisLevel = dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetAnalysisLevel(); \r
+       \r
+       TTree* tree = dynamic_cast<TTree*> (GetInputData(0));\r
+       if (!tree) \r
+       {\r
+               Printf("ERROR: Could not read chain from input slot 0");\r
+       } \r
+       else \r
+       {\r
+               if(gAnalysisLevel == "ESD") \r
+               {\r
+                       AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());   \r
+                       if (!esdH) \r
+                       {\r
+                               Printf("ERROR: Could not get ESDInputHandler");\r
+                       } \r
+                       else\r
+                               fESD = esdH->GetEvent();\r
+                       AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());\r
+                       if (!mcH) \r
+                       {\r
+                               Printf("ERROR: Could not retrieve MC event handler");\r
+                       }\r
+                       else\r
+                               fMC = mcH->MCEvent();\r
+               }\r
+               else if(gAnalysisLevel == "AOD") \r
+               {\r
+                       AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());        \r
+                       if (!aodH) \r
+                       {\r
+                               Printf("ERROR: Could not get AODInputHandler");\r
+                       }\r
+                        else\r
+                               fAOD = aodH->GetEvent();\r
+               }\r
+               else\r
+                       Printf("Wrong analysis type: Only ESD, AOD  types are allowed!");\r
+       }\r
+}\r
+//________________________________________________________________________\r
+void AliProtonFeedDownAnalysisTask::CreateOutputObjects() \r
+{\r
+       // Create output objects\r
+       // Called once\r
+       fList = new TList();\r
+       fList->Add(fProtonAnalysis->GetProtonContainer());\r
+       fList->Add(fProtonAnalysis->GetAntiProtonContainer());\r
+       fList->Add(fProtonAnalysis->GetLambdaHist());\r
+       fList->Add(fProtonAnalysis->GetLambdaweightedHist());\r
+       fList->Add(fProtonAnalysis->GetAntiLambdaHist());\r
+       fList->Add(fProtonAnalysis->GetAntiLambdaweightedHist());\r
+       fStatHist=new TH1F("StatsHist","StatsHist",10,-0.5,9.5);\r
+       fList->Add(fStatHist);\r
+       fStatHist->GetXaxis()->SetBinLabel(1,"level1cut");\r
+       fStatHist->GetXaxis()->SetBinLabel(2,"level2cut");\r
+       fStatHist->GetXaxis()->SetBinLabel(3,"level3cut");\r
+       fStatHist->GetXaxis()->SetBinLabel(4,"level4cut");\r
+}\r
+//________________________________________________________________________\r
+void AliProtonFeedDownAnalysisTask::Exec(Option_t *) \r
+{\r
+       // Main loop\r
+       // Called for each event\r
+       TString gAnalysisLevel = dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetAnalysisLevel(); \r
+       //TString gAnalysisLevel = (fProtonAnalysis->GetProtonAnalysisBaseObject())->GetAnalysisLevel(); \r
+       if(gAnalysisLevel == "ESD") \r
+       {\r
+               if (!fESD) \r
+               {\r
+                       Printf("ERROR: fESD not available");\r
+                       return;\r
+               }\r
+               fStatHist->Fill(0);\r
+               if(dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->IsEventTriggered(fESD,dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetTriggerMode())) \r
+               {\r
+                       fStatHist->Fill(1);\r
+                       const AliESDVertex *vertex = dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetVertex(fESD,dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetAnalysisMode(),dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetVxMax(),dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetVyMax(),dynamic_cast<AliProtonAnalysisBase*>(fProtonAnalysis->GetProtonAnalysisBaseObject())->GetVzMax());\r
+                       if(vertex) \r
+                       {\r
+                               fStatHist->Fill(2);\r
+                               Printf("Proton ESD analysis task: There are %d tracks in this event", fESD->GetNumberOfTracks());\r
+                               AliStack* stack1=0x0;\r
+                               if(!fMC)\r
+                                       return;                 \r
+                               stack1 = fMC->Stack();\r
+                               if(!stack1)\r
+                                       return;\r
+                               fStatHist->Fill(3);     \r
+                               fProtonAnalysis->Analyze(fESD,vertex,stack1);\r
+                               fProtonAnalysis->Analyze(stack1);\r
+\r
+                       }//reconstructed vertex\r
+               }//triggered event\r
+       }//ESD analysis              \r
+       \r
+       else if(gAnalysisLevel == "AOD") \r
+       {\r
+               if (!fAOD) \r
+               {\r
+                       Printf("ERROR: fAOD not available");\r
+                       return;\r
+               }\r
+               Printf("Proton AOD analysis task: There are %d tracks in this event", fAOD->GetNumberOfTracks());\r
+               fProtonAnalysis->Analyze(fAOD);\r
+       }//AOD analysis\r
+                        \r
+       \r
+       // Post output data.\r
+       PostData(0, fList);\r
+}    \r
+//__________________________________________________________________________________________\r
+void AliProtonFeedDownAnalysisTask::Terminate(Option_t *) \r
+{\r
+\r
+}\r
+\r
+\r
diff --git a/PWG2/SPECTRA/AliProtonFeedDownAnalysisTask.h b/PWG2/SPECTRA/AliProtonFeedDownAnalysisTask.h
new file mode 100644 (file)
index 0000000..c9de5c0
--- /dev/null
@@ -0,0 +1,46 @@
+#ifndef ALIPROTONFEEDDOWNANALYSISTASK_H\r
+#define ALIPROTONFEEDDOWNANALYSISTASK_H\r
+#include "AliAnalysisTask.h"\r
+\r
+class TList;\r
+class AliESDEvent;\r
+class AliAODEvent;\r
+class AliMCEvent;\r
+class AliProtonFeedDownAnalysis;\r
+\r
+\r
+\r
+class AliProtonFeedDownAnalysisTask : public AliAnalysisTask {\r
+ public:\r
+       AliProtonFeedDownAnalysisTask();\r
+       AliProtonFeedDownAnalysisTask(const char *name);\r
+       virtual ~AliProtonFeedDownAnalysisTask() {}\r
+       \r
+       virtual void   ConnectInputData(Option_t *);\r
+       virtual void   CreateOutputObjects();\r
+       virtual void   Exec(Option_t *option);\r
+       virtual void   Terminate(Option_t *);\r
+       \r
+       void SetAnalysisObject(AliProtonFeedDownAnalysis *const analysis) {fProtonAnalysis = analysis;}\r
+  \r
+ private:\r
+       AliESDEvent *fESD;    //ESD object \r
+       AliAODEvent *fAOD;    //AOD object\r
+       AliMCEvent  *fMC;     //MC object \r
+       \r
+       TList  *fList; //TList output object \r
+       \r
+       AliProtonFeedDownAnalysis *fProtonAnalysis; //analysis object \r
+       \r
+       TH1F *fStatHist;\r
+       \r
+       AliProtonFeedDownAnalysisTask(const AliProtonFeedDownAnalysisTask&); // not implemented\r
+       AliProtonFeedDownAnalysisTask& operator=(const AliProtonFeedDownAnalysisTask&); // not implemented\r
+       \r
+ClassDef(AliProtonFeedDownAnalysisTask, 1); // example of analysis\r
+};\r
+\r
+#endif\r
+\r
+\r
+\r
diff --git a/PWG2/SPECTRA/macros/AddTaskProtonFeedDownAnalysis.C b/PWG2/SPECTRA/macros/AddTaskProtonFeedDownAnalysis.C
new file mode 100644 (file)
index 0000000..1dba0cd
--- /dev/null
@@ -0,0 +1,41 @@
+AliProtonFeedDownAnalysisTask* AddTaskProtonFeedDownAnalysis(const char *analysisType="Hybrid",const char *pidMode="Bayesian")
+{
+  // Creates a proton analysis task and adds it to the analysis manager.
+  
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskProtons", "No analysis manager to connect to.");
+    return NULL;
+  }   
+  
+  // Check the analysis type using the event handlers connected to the analysis manager.
+  //==============================================================================
+  if (!mgr->GetInputEventHandler()) {
+    ::Error("AddTaskProtons", "This task requires an input event handler");
+    return NULL;
+  }   
+  TString type = mgr->GetInputEventHandler()->GetDataType(); // can be "ESD" or "AOD"
+  gROOT->LoadMacro("$ALICE_ROOT/PWG2/SPECTRA/macros/configProtonFeedDownAnalysis.C");
+  AliProtonFeedDownAnalysis *pa = 0;
+  if (type=="ESD") pa = GetProtonFeedDownAnalysisObject("ESD", analysisType, pidMode);
+  else if (type=="AOD") pa = GetProtonFeedDownAnalysisObject("AOD", analysisType, pidMode);
+  else return NULL;
+
+  // Create the task, add it to manager and configure it.
+  //===========================================================================
+  AliProtonFeedDownAnalysisTask *taskproton = new  AliProtonFeedDownAnalysisTask("TaskProtons");
+  mgr->AddTask(taskproton);
+  taskproton->SetAnalysisObject(pa);
+
+  // Create ONLY the output containers for the data produced by the task.
+  // Get and connect other common input/output containers via the manager as below
+  //==============================================================================
+  AliAnalysisDataContainer *cout_proton = mgr->CreateContainer("protonFeedDown", TList::Class(),AliAnalysisManager::kOutputContainer,"outputAliProtonFeedDownAnalysisTask.root");                              
+  mgr->ConnectInput(taskproton, 0, mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(taskproton, 0, cout_proton);
+  
+  // Return task pointer at the end
+  return taskproton;
+}
diff --git a/PWG2/SPECTRA/macros/configProtonFeedDownAnalysis.C b/PWG2/SPECTRA/macros/configProtonFeedDownAnalysis.C
new file mode 100644 (file)
index 0000000..d993341
--- /dev/null
@@ -0,0 +1,154 @@
+//__________________________________________________________
+
+Double_t test1(Double_t *x,Double_t *p)
+{
+       /*if(x[0]<1.5) 
+               return -1.0*x[0]+2.5; 
+       else 
+               return 1.0;*/
+       if(x[0]<3.0) 
+               return -0.2*x[0]+1.6; 
+       else 
+               return 1.0;     
+}
+//_____________________________________________________________
+Double_t test2(Double_t *x,Double_t *p)
+{
+       if(x[0]<1.5) 
+               return -1.0*x[0]+2.5; 
+       else 
+               return 1.0;
+}
+
+//__________________________________________________//
+AliProtonFeedDownAnalysis *GetProtonFeedDownAnalysisObject(const char* analysisLevel = "ESD", 
+                                          const char* esdAnalysisType = "TPC", 
+                                          const char* pidMode = "Bayesian") {
+  //Function to setup the AliProtonAnalysis object and return it
+  AliProtonAnalysisBase *baseAnalysis = GetProtonAnalysisBaseObject(analysisLevel,esdAnalysisType,pidMode);
+
+  AliProtonFeedDownAnalysis *analysis = new AliProtonFeedDownAnalysis();
+  analysis->SetBaseAnalysis(baseAnalysis);
+  //if(analysisBase->GetEtaMode()) analysis->SetEtaMode();
+  analysis->InitAnalysisHistograms(baseAnalysis->GetNBinsX(),
+                                  baseAnalysis->GetMinX(),
+                                  baseAnalysis->GetMaxX(),
+                                  baseAnalysis->GetNBinsY(),
+                                  baseAnalysis->GetMinY(),
+                                  baseAnalysis->GetMaxY());
+TF1* weightfunction=new TF1("weightfunction","1");
+//TF1* weightfunction=new TF1("weightfunction",test1,0.5,4.0,0);       
+//TF1* weightfunction=new TF1("weightfunction",test2,0.5,4.0,0);                          
+ analysis->SetWeightFunction(weightfunction);     
+    
+  return analysis;
+}
+
+
+//__________________________________________________//
+AliProtonAnalysisBase *GetProtonAnalysisBaseObject(const char* analysisLevel = "ESD",
+                                                  const char* esdAnalysisType = "TPC",
+                                                  const char* pidMode = "Bayesian") {
+  //Function to setup the AliProtonAnalysisBase object and return it
+  AliProtonAnalysisBase *baseAnalysis = new AliProtonAnalysisBase();
+  //baseAnalysis->SetDebugMode();
+  baseAnalysis->SetAnalysisLevel(analysisLevel);
+  if(analysisLevel == "ESD") {  
+    baseAnalysis->SetTriggerMode(AliProtonAnalysisBase::kMB2);
+    baseAnalysis->SetMinTPCClusters(110);
+    baseAnalysis->SetMaxChi2PerTPCCluster(2.2);
+    baseAnalysis->SetMaxCov11(0.5);
+    baseAnalysis->SetMaxCov22(0.5);
+    baseAnalysis->SetMaxCov33(0.5);
+    baseAnalysis->SetMaxCov44(0.5);
+    baseAnalysis->SetMaxCov55(0.5);
+    baseAnalysis->SetMinTPCdEdxPoints(80);
+    switch(esdAnalysisType) {
+    case "TPC":
+      baseAnalysis->SetAnalysisMode(AliProtonAnalysisBase::kTPC);
+      baseAnalysis->SetPhaseSpace(10, -0.5, 0.5, 16, 0.5, 0.9);
+      baseAnalysis->SetTPCpid();
+      baseAnalysis->SetMaxSigmaToVertexTPC(2.0);
+      //baseAnalysis->SetMaxDCAXYTPC(1.5);
+      //baseAnalysis->SetMaxDCAZTPC(1.5);
+      break;
+    case "Hybrid":
+      baseAnalysis->SetAnalysisMode(AliProtonAnalysisBase::kHybrid);
+      baseAnalysis->SetPhaseSpace(10, -0.5, 0.5, 16, 0.5, 0.9);
+      baseAnalysis->SetTPCpid();
+      baseAnalysis->SetMaxSigmaToVertex(2.0);
+      /*baseAnalysis->SetMaxDCAXY(1.5);
+       baseAnalysis->SetMaxDCAZ(1.5);*/
+      baseAnalysis->SetPointOnITSLayer6();
+      baseAnalysis->SetPointOnITSLayer5();
+      //baseAnalysis->SetPointOnITSLayer4();
+      //baseAnalysis->SetPointOnITSLayer3();
+      baseAnalysis->SetPointOnITSLayer2();
+      baseAnalysis->SetPointOnITSLayer1();
+      baseAnalysis->SetMinITSClusters(4);
+      break;
+    case "Global":
+      baseAnalysis->SetAnalysisMode(AliProtonAnalysisBase::kGlobal);
+      baseAnalysis->SetPhaseSpace(20, -1.0, 1.0, 48, 0.3, 1.5);
+      baseAnalysis->SetMaxSigmaToVertex(2.0);
+      //baseAnalysis->SetMaxDCAXY(2.0);
+      //baseAnalysis->SetMaxDCAZ(2.0);
+      baseAnalysis->SetTPCRefit();
+      baseAnalysis->SetPointOnITSLayer1();
+      baseAnalysis->SetPointOnITSLayer2();
+      //baseAnalysis->SetPointOnITSLayer3();
+      //baseAnalysis->SetPointOnITSLayer4();
+      baseAnalysis->SetPointOnITSLayer5();
+      baseAnalysis->SetPointOnITSLayer6();
+      baseAnalysis->SetMinITSClusters(5);
+      baseAnalysis->SetITSRefit();
+      baseAnalysis->SetESDpid();
+      baseAnalysis->SetTOFpid();
+      break;
+    default:
+      break;
+    }
+    baseAnalysis->SetAcceptedVertexDiamond(5.,5.,15.);
+    baseAnalysis->SetEtaMode();
+    switch(pidMode) {
+    case "Bayesian":
+      baseAnalysis->SetPIDMode(AliProtonAnalysisBase::kBayesian);
+      //Momentum dependent priors
+      /*TFile *f = TFile::Open("$ALICE_ROOT/PWG2/data/PriorProbabilities.root ");
+       TF1 *fitElectrons = (TF1 *)f->Get("fitElectrons");
+       TF1 *fitMuons = (TF1 *)f->Get("fitMuons");
+       TF1 *fitPions = (TF1 *)f->Get("fitPions");
+       TF1 *fitKaons = (TF1 *)f->Get("fitKaons");
+       TF1 *fitProtons = (TF1 *)f->Get("fitProtons");
+       baseAnalysis->SetPriorProbabilityFunctions(fitElectrons,
+       fitMuons,
+       fitPions,
+       fitKaons,
+       fitProtons);*/
+      //Fixed prior probabilities
+      Double_t partFrac[10] = {0.01, 0.01, 0.85, 0.10, 0.05,0.0,0.0,0.0,0.0,0.0};
+      if(!baseAnalysis->IsPriorProbabilityFunctionUsed())
+       baseAnalysis->SetPriorProbabilities(partFrac);
+      break;
+    case "Ratio":
+      baseAnalysis->SetPIDMode(AliProtonAnalysisBase::kRatio);
+      break;
+    case "Sigma1":
+      baseAnalysis->SetPIDMode(AliProtonAnalysisBase::kSigma1);
+      baseAnalysis->SetNSigma(3);
+      baseAnalysis->SetdEdxBandInfo("$ALICE_ROOT/PWG2/data/protonsdEdxInfo.dat");
+      break;
+    case "Sigma2":
+      baseAnalysis->SetPIDMode(AliProtonAnalysisBase::kSigma2);
+      baseAnalysis->SetNSigma(3);
+      baseAnalysis->SetdEdxBandInfo("$ALICE_ROOT/PWG2/data/protonsdEdxInfo.dat");
+      break;
+    default:
+      break;
+    }//PID mode
+  }//ESD
+  if(analysisLevel == "MC") 
+    baseAnalysis->SetPhaseSpace(10, -0.5, 0.5, 16, 0.5, 0.9);
+
+  return baseAnalysis;
+}