SPECTRA/PiKaPr/TestAOD/AliSpectraAODPID.cxx
SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraAOD.cxx
SPECTRA/Kinks/AliAnalysisKinkESDat.cxx
+ SPECTRA/PiKaPr/TestAOD/AliSpectraBothTrackCuts.cxx
+ SPECTRA/PiKaPr/TestAOD/AliSpectraBothEventCuts.cxx
+ SPECTRA/PiKaPr/TestAOD/AliSpectraBothHistoManager.cxx
+ SPECTRA/PiKaPr/TestAOD/AliSpectraBothPID.cxx
+ SPECTRA/PiKaPr/TestAOD/AliAnalysisTaskSpectraBoth.cxx
+
)
string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
#pragma link C++ class AliSpectraAODTrackCuts+;
#pragma link C++ class AliAnalysisKinkESDat+;
+#pragma link C++ class AliAnalysisTaskSpectraBoth+;
+#pragma link C++ class AliSpectraBothEventCuts+;
+#pragma link C++ class AliSpectraBothHistoManager+;
+#pragma link C++ class AliSpectraBothPID+;
+#pragma link C++ class AliSpectraBothTrackCuts+;
+
#endif
--- /dev/null
+/**************************************************************************\r
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *\r
+ * *\r
+ * Author: The ALICE Off-line Project. *\r
+ * Contributors are mentioned in the code where appropriate. *\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
+// AliAnalysisTaskSpectraBoth class\r
+//-----------------------------------------------------------------\r
+\r
+#include "TChain.h"\r
+#include "TTree.h"\r
+#include "TLegend.h"\r
+#include "TH1F.h"\r
+#include "TH2F.h"\r
+#include "TCanvas.h"\r
+#include "AliAnalysisTask.h"\r
+#include "AliAnalysisManager.h"\r
+#include "AliAODTrack.h"\r
+#include "AliAODMCParticle.h"\r
+#include "AliVParticle.h"\r
+#include "AliAODEvent.h"\r
+#include "AliAODInputHandler.h"\r
+#include "AliAnalysisTaskSpectraBoth.h"\r
+#include "AliAnalysisTaskESDfilter.h"\r
+#include "AliAnalysisDataContainer.h"\r
+#include "AliSpectraBothHistoManager.h"\r
+#include "AliSpectraBothTrackCuts.h"\r
+#include "AliSpectraBothEventCuts.h"\r
+#include "AliCentrality.h"\r
+#include "TProof.h"\r
+#include "AliPID.h"\r
+#include "AliVEvent.h"\r
+#include "AliESDEvent.h"\r
+#include "AliPIDResponse.h"\r
+#include "AliStack.h"\r
+#include "AliSpectraBothPID.h"\r
+#include <TMCProcess.h>\r
+\r
+#include <iostream>\r
+\r
+\r
+\r
+\r
+using namespace AliSpectraNameSpaceBoth;\r
+using namespace std;\r
+\r
+ClassImp(AliAnalysisTaskSpectraBoth)\r
+\r
+//________________________________________________________________________\r
+AliAnalysisTaskSpectraBoth::AliAnalysisTaskSpectraBoth(const char *name) : AliAnalysisTaskSE(name), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0), fPID(0), fIsMC(0), fNRebin(0),fUseMinSigma(0),fCuts(0),fdotheMCLoopAfterEventCuts(0)\r
+{\r
+ // Default constructor\r
+ \r
+ DefineInput(0, TChain::Class());\r
+ DefineOutput(1, AliSpectraBothHistoManager::Class());\r
+ DefineOutput(2, AliSpectraBothEventCuts::Class());\r
+ DefineOutput(3, AliSpectraBothTrackCuts::Class());\r
+ DefineOutput(4, AliSpectraBothPID::Class());\r
+ fNRebin=0;\r
+ \r
+}\r
+//________________________________________________________________________\r
+//________________________________________________________________________\r
+void AliAnalysisTaskSpectraBoth::UserCreateOutputObjects()\r
+{\r
+ // create output objects\r
+ fHistMan = new AliSpectraBothHistoManager("SpectraHistos",fNRebin);\r
+\r
+ if (!fTrackCuts) AliFatal("Track Cuts should be set in the steering macro");\r
+ if (!fEventCuts) AliFatal("Event Cuts should be set in the steering macro");\r
+ if (!fPID) AliFatal("PID object should be set in the steering macro");\r
+ fTrackCuts->SetAliESDtrackCuts(fCuts);\r
+ PostData(1, fHistMan );\r
+ PostData(2, fEventCuts);\r
+ PostData(3, fTrackCuts);\r
+ PostData(4, fPID );\r
+\r
+}\r
+//________________________________________________________________________\r
+void AliAnalysisTaskSpectraBoth::UserExec(Option_t *)\r
+{\r
+ // main event loop\r
+ Int_t ifAODEvent=AliSpectraBothTrackCuts::kotherobject;\r
+ fAOD = dynamic_cast<AliVEvent*>(InputEvent());\r
+ \r
+ // AliESDEvent* esdevent=0x0;\r
+ // AliAODEvent* aodevent=0x0;\r
+ \r
+ TString nameoftrack(fAOD->ClassName()); \r
+ if(!nameoftrack.CompareTo("AliESDEvent"))\r
+ {\r
+ ifAODEvent=AliSpectraBothTrackCuts::kESDobject;\r
+ //esdevent=dynamic_cast<AliESDEvent*>(fAOD);\r
+ }\r
+ else if(!nameoftrack.CompareTo("AliAODEvent"))\r
+ {\r
+ ifAODEvent=AliSpectraBothTrackCuts::kAODobject;\r
+ //aodevent=dynamic_cast<AliAODEvent*>(fAOD);\r
+ }\r
+ else\r
+ AliFatal("Not processing AODs or ESDS") ;\r
+ if(fdotheMCLoopAfterEventCuts)\r
+ if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
+ \r
+ \r
+ TClonesArray *arrayMC = 0;\r
+ Int_t npar=0;\r
+ AliStack* stack=0x0;\r
+ \r
+ if (fIsMC)\r
+ {\r
+ \r
+ if(ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
+ { \r
+ arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());\r
+ if (!arrayMC) \r
+ {\r
+ AliFatal("Error: MC particles branch not found!\n");\r
+ }\r
+ Int_t nMC = arrayMC->GetEntries();\r
+ for (Int_t iMC = 0; iMC < nMC; iMC++)\r
+ {\r
+ AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(iMC);\r
+ if(!partMC->Charge()) continue;//Skip neutrals\r
+ //if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax()){//charged hadron are filled inside the eta acceptance\r
+ //Printf("%f %f-%f",partMC->Eta(),fTrackCuts->GetEtaMin(),fTrackCuts->GetEtaMax());\r
+ if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())\r
+ fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
+ \r
+ //rapidity cut\r
+ if(TMath::Abs(partMC->Y()) > fTrackCuts->GetY() ) continue; \r
+ if(partMC->IsPhysicalPrimary())\r
+ npar++; \r
+ // check for true PID + and fill P_t histos \r
+ Int_t charge = partMC->Charge() > 0 ? kChPos : kChNeg ;\r
+ Int_t id = fPID->GetParticleSpecie(partMC);\r
+ if(id != kSpUndefined) \r
+ {\r
+ fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),partMC->IsPhysicalPrimary());\r
+ }\r
+ }\r
+ }\r
+ if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
+ { \r
+ AliMCEvent* mcEvent = (AliMCEvent*) MCEvent();\r
+ //Printf("MC particles: %d", mcEvent->GetNumberOfTracks()); \r
+ if (!mcEvent) \r
+ {\r
+ AliFatal("Error: MC particles branch not found!\n");\r
+ }\r
+ stack = mcEvent->Stack();\r
+ Int_t nMC = stack->GetNtrack();\r
+ for (Int_t iMC = 0; iMC < nMC; iMC++)\r
+ {\r
+ \r
+ TParticle *partMC = stack->Particle(iMC);\r
+ \r
+ if(!partMC) \r
+ continue; \r
+ \r
+ if(!partMC->GetPDG(0))\r
+ continue;\r
+ if(TMath::Abs(partMC->GetPDG(0)->Charge()/3.0)<0.01) \r
+ continue;//Skip neutrals\r
+ if(partMC->Eta() > fTrackCuts->GetEtaMin() && partMC->Eta() < fTrackCuts->GetEtaMax())\r
+ fHistMan->GetPtHistogram(kHistPtGen)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));\r
+ if(TMath::Abs(partMC->Y()) > fTrackCuts->GetY() ) \r
+ continue;\r
+ if(stack->IsPhysicalPrimary(iMC))\r
+ npar++; \r
+ // check for true PID + and fill P_t histos \r
+ Int_t charge = partMC->GetPDG(0)->Charge()/3.0 > 0 ? kChPos : kChNeg ;\r
+ Int_t id = fPID->GetParticleSpecie(partMC);\r
+ if(id != kSpUndefined) \r
+ {\r
+ fHistMan->GetHistogram2D(kHistPtGenTruePrimary,id,charge)->Fill(partMC->Pt(),stack->IsPhysicalPrimary(iMC));\r
+ }\r
+ }\r
+ }\r
+ }\r
+ \r
+ if(!fdotheMCLoopAfterEventCuts)\r
+ if(!fEventCuts->IsSelected(fAOD,fTrackCuts))return;//event selection\r
+\r
+ //main loop on tracks\r
+ Int_t ntracks=0;\r
+ //cout<<fAOD->GetNumberOfTracks()<<endl;\r
+ for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) \r
+ {\r
+ AliVTrack* track = dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));\r
+ AliAODTrack* aodtrack=0;\r
+ AliESDtrack* esdtrack=0;\r
+ Float_t dca=-999.;\r
+ if(ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
+ {\r
+ esdtrack=dynamic_cast<AliESDtrack*>(track);\r
+ Float_t dcaz=0.0;\r
+ esdtrack->GetImpactParameters(dca,dcaz);\r
+ }\r
+ else if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
+ {\r
+ aodtrack=dynamic_cast<AliAODTrack*>(track); \r
+ dca=aodtrack->DCA();\r
+ }\r
+ else\r
+ continue;\r
+ if (!fTrackCuts->IsSelected(track,kTRUE)) \r
+ continue;\r
+ ntracks++;\r
+ fPID->FillQAHistos(fHistMan, track, fTrackCuts);\r
+ \r
+ //calculate DCA for AOD track\r
+ if(dca==-999.)\r
+ {// track->DCA() does not work in old AOD production\r
+ Double_t d[2], covd[3];\r
+ AliVTrack* track_clone=(AliVTrack*)track->Clone("track_clone"); // need to clone because PropagateToDCA updates the track parameters\r
+ Bool_t isDCA = track_clone->PropagateToDCA(fAOD->GetPrimaryVertex(),fAOD->GetMagneticField(),9999.,d,covd);\r
+ delete track_clone;\r
+ if(!isDCA)\r
+ d[0]=-999.;\r
+ dca=d[0];\r
+ }\r
+ fHistMan->GetPtHistogram(kHistPtRec)->Fill(track->Pt(),dca); // PT histo\r
+ // get identity and charge\r
+ Bool_t rec[3]={false,false,false};\r
+ Int_t idRec = fPID->GetParticleSpecie(fHistMan,track, fTrackCuts,rec);\r
+ for(int irec=kSpPion;irec<kNSpecies;irec++)\r
+ {\r
+ \r
+ if(fUseMinSigma)\r
+ {\r
+ if(irec>kSpPion)\r
+ break;\r
+ }\r
+ else\r
+ { \r
+ if(!rec[irec]) \r
+ idRec = kSpUndefined;\r
+ else \r
+ idRec=irec;\r
+ } \r
+ \r
+ Int_t charge = track->Charge() > 0 ? kChPos : kChNeg;\r
+ \r
+ // Fill histograms, only if inside y and nsigma acceptance\r
+ if(idRec != kSpUndefined && fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec))\r
+ fHistMan->GetHistogram2D(kHistPtRecSigma,idRec,charge)->Fill(track->Pt(),dca);\r
+ //can't put a continue because we still have to fill allcharged primaries, done later\r
+ \r
+ /* MC Part */\r
+ if (arrayMC||stack) \r
+ {\r
+ Bool_t isPrimary = kFALSE;\r
+ Bool_t isSecondaryMaterial = kFALSE; \r
+ Bool_t isSecondaryWeak = kFALSE; \r
+ Int_t idGen =kSpUndefined;\r
+ Int_t pdgcode=0;\r
+ if (ifAODEvent==AliSpectraBothTrackCuts::kAODobject)\r
+ {\r
+ AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));\r
+ if (!partMC) \r
+ { \r
+ AliError("Cannot get MC particle");\r
+ continue; \r
+ }\r
+ // Check if it is primary, secondary from material or secondary from weak decay\r
+ isPrimary = partMC->IsPhysicalPrimary();\r
+ isSecondaryWeak = partMC->IsSecondaryFromWeakDecay();\r
+ isSecondaryMaterial = partMC->IsSecondaryFromMaterial();\r
+ //cout<<"AOD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;\r
+\r
+ if(!isPrimary&&!isSecondaryWeak&&!isSecondaryMaterial) \r
+ {\r
+ AliError("old tagging");\r
+ Int_t mfl=-999,codemoth=-999;\r
+ Int_t indexMoth=partMC->GetMother(); // FIXME ignore fakes? TO BE CHECKED, on ESD is GetFirstMother()\r
+ if(indexMoth>=0)\r
+ {//is not fake\r
+ AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);\r
+ codemoth = TMath::Abs(moth->GetPdgCode());\r
+ mfl = Int_t (codemoth/ TMath::Power(10, Int_t(TMath::Log10(codemoth))));\r
+ }\r
+ //Int_t uniqueID = partMC->GetUniqueID();\r
+ //cout<<"uniqueID: "<<partMC->GetUniqueID()<<" "<<kPDecay<<endl;\r
+ //cout<<"status: "<<partMC->GetStatus()<<" "<<kPDecay<<endl;\r
+ // if(uniqueID == kPDecay)Printf("!!!!!!!!!!!!!!!!!!!!!!!!!!!!!");\r
+ if(mfl==3) \r
+ isSecondaryWeak = kTRUE; // add if(partMC->GetStatus() & kPDecay)? FIXME\r
+ else \r
+ isSecondaryMaterial = kTRUE;\r
+ }\r
+ //cout<<"AOD 2 tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<" "<<partMC->GetMCProcessCode()<<endl;\r
+\r
+\r
+ idGen = fPID->GetParticleSpecie(partMC);\r
+ pdgcode=partMC->GetPdgCode(); \r
+ }\r
+ else if (ifAODEvent==AliSpectraBothTrackCuts::kESDobject)\r
+ {\r
+ TParticle *partMC =stack->Particle(TMath::Abs(track->GetLabel()));\r
+ if (!partMC) \r
+ { \r
+ AliError("Cannot get MC particle");\r
+ continue; \r
+ }\r
+ isPrimary = stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()));\r
+ isSecondaryWeak = stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()));\r
+ isSecondaryMaterial = stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()));\r
+ //cout<<"ESD tagging "<<isPrimary<<" "<<isSecondaryWeak<<isSecondaryMaterial<<endl;\r
+ \r
+ idGen = fPID->GetParticleSpecie(partMC);\r
+ pdgcode=partMC->GetPdgCode(); \r
+ }\r
+ else\r
+ return;\r
+ \r
+ // cout<<isPrimary<<" "<<isSecondaryWeak<<" "<<isSecondaryMaterial<<endl;\r
+ // cout<<" functions "<<partMC->IsPhysicalPrimary()<<" "<<partMC->IsSecondaryFromWeakDecay()<<" "<<partMC->IsSecondaryFromMaterial()<<endl;\r
+ \r
+ if (isPrimary&&irec==kSpPion)\r
+ fHistMan->GetPtHistogram(kHistPtRecPrimary)->Fill(track->Pt(),dca); // PT histo of primaries\r
+ \r
+ //nsigma cut (reconstructed nsigma)\r
+ if(idRec == kSpUndefined) \r
+ continue;\r
+ \r
+ // rapidity cut (reconstructed pt and identity)\r
+ if(!fTrackCuts->CheckYCut ((BothParticleSpecies_t)idRec)) continue;\r
+ \r
+ // Get true ID\r
+ \r
+ //if(TMath::Abs(partMC->Y()) > fTrackCuts->GetY() ) continue; // FIXME: do we need a rapidity cut on the generated?\r
+ // Fill histograms for primaries\r
+ \r
+ if (idRec == idGen) fHistMan->GetHistogram2D(kHistPtRecTrue, idGen, charge)->Fill(track->Pt(),dca); \r
+ \r
+ if (isPrimary) \r
+ {\r
+ fHistMan->GetHistogram2D(kHistPtRecSigmaPrimary, idRec, charge)->Fill(track->Pt(),dca); \r
+ if(idGen != kSpUndefined) \r
+ {\r
+ fHistMan->GetHistogram2D(kHistPtRecPrimary, idGen, charge)->Fill(track->Pt(),dca);\r
+ if (idRec == idGen) \r
+ fHistMan->GetHistogram2D(kHistPtRecTruePrimary, idGen, charge)->Fill(track->Pt(),dca); \r
+ }\r
+ }\r
+ //25th Apr - Muons are added to Pions -- FIXME\r
+ if ( pdgcode == 13 && idRec == kSpPion) \r
+ { \r
+ fHistMan->GetPtHistogram(kHistPtRecTrueMuonPlus)->Fill(track->Pt(),dca); \r
+ if(isPrimary)\r
+ fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonPlus)->Fill(track->Pt(),dca); \r
+ }\r
+ if ( pdgcode == -13 && idRec == kSpPion) \r
+ { \r
+ fHistMan->GetPtHistogram(kHistPtRecTrueMuonMinus)->Fill(track->Pt(),dca); \r
+ if (isPrimary) \r
+ {\r
+ fHistMan->GetPtHistogram(kHistPtRecTruePrimaryMuonMinus)->Fill(track->Pt(),dca); \r
+ }\r
+ }\r
+ \r
+ ///..... END FIXME\r
+ \r
+ // Fill secondaries\r
+ if(isSecondaryWeak ) \r
+ fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryWeakDecay, idRec, charge)->Fill(track->Pt(),dca);\r
+ if(isSecondaryMaterial) \r
+ fHistMan->GetHistogram2D(kHistPtRecSigmaSecondaryMaterial , idRec, charge)->Fill(track->Pt(),dca);\r
+ \r
+ }//end if(arrayMC)\r
+ }\r
+ \r
+ \r
+ } // end loop on tracks\r
+ // cout<< ntracks<<endl;\r
+ fHistMan->GetGenMulvsRawMulHistogram("hHistGenMulvsRawMul")->Fill(npar,ntracks);\r
+ PostData(1, fHistMan );\r
+ PostData(2, fEventCuts);\r
+ PostData(3, fTrackCuts);\r
+ PostData(4, fPID );\r
+}\r
+\r
+//_________________________________________________________________\r
+void AliAnalysisTaskSpectraBoth::Terminate(Option_t *)\r
+{\r
+ // Terminate\r
+}\r
--- /dev/null
+#ifndef ALIANALYSISTASKSPECTRABOTH_H\r
+#define ALIANALYSISTASKSPECTRABOTH_H\r
+\r
+/* See cxx source for full Copyright notice */\r
+\r
+//-------------------------------------------------------------------------\r
+// AliAnalysisTaskSpectraBoth\r
+//\r
+//\r
+//\r
+//\r
+// Author: Michele Floris, CERN\r
+//-------------------------------------------------------------------------\r
+\r
+class TH1F;\r
+class TH2F;\r
+class AliVEvent;\r
+class AliSpectraBothHistoManager;\r
+class AliSpectraBothTrackCuts;\r
+class AliSpectraBothEventCuts;\r
+class AliSpectraBothPID;\r
+class AliESDtrackCuts;\r
+#include "AliSpectraBothHistoManager.h"\r
+#include "AliAnalysisTaskSE.h"\r
+#include "AliESDtrackCuts.h"\r
+\r
+class AliAnalysisTaskSpectraBoth : public AliAnalysisTaskSE\r
+{\r
+public:\r
+\r
+ // constructors\r
+ AliAnalysisTaskSpectraBoth() : AliAnalysisTaskSE(), fAOD(0), fHistMan(0), fTrackCuts(0), fEventCuts(0), fPID(0), fIsMC(0), fNRebin(0),fUseMinSigma(0),fCuts(0),fdotheMCLoopAfterEventCuts(0)\r
+ {}\r
+ AliAnalysisTaskSpectraBoth(const char *name);\r
+ virtual ~AliAnalysisTaskSpectraBoth() {}\r
+\r
+ void SetIsMC(Bool_t isMC = kFALSE) {fIsMC = isMC; };\r
+ Bool_t GetIsMC() const { return fIsMC;};\r
+\r
+ virtual void UserCreateOutputObjects();\r
+ virtual void UserExec(Option_t *option);\r
+ virtual void Terminate(Option_t *);\r
+\r
+ AliSpectraBothHistoManager * GetHistoManager() { return fHistMan; }\r
+ AliSpectraBothTrackCuts * GetTrackCuts() { return fTrackCuts; }\r
+ AliSpectraBothEventCuts * GetEventCuts() { return fEventCuts; }\r
+ AliSpectraBothPID * GetPID() { return fPID; }\r
+ \r
+ void SetTrackCuts(AliSpectraBothTrackCuts * tc) { fTrackCuts = tc; }\r
+ void SetEventCuts(AliSpectraBothEventCuts * vc) { fEventCuts = vc; }\r
+ void SetPID (AliSpectraBothPID * pid) { fPID = pid; }\r
+ void SetNRebin(Int_t nreb){fNRebin=nreb;}\r
+ void SetUseMinSigma (Bool_t flag) {fUseMinSigma=flag;}\r
+ Int_t GetNRebin() const {return fNRebin;}\r
+ void SetAliESDtrackCuts(AliESDtrackCuts* cuts ){fCuts=cuts;}\r
+ void SetdotheMCLoopAfterEventCuts (Bool_t flag) {fdotheMCLoopAfterEventCuts=flag;}\r
+ Bool_t GetdotheMCLoopAfterEventCuts () const {return fdotheMCLoopAfterEventCuts;}\r
+private:\r
+\r
+ AliVEvent * fAOD; //! AOD object\r
+ AliSpectraBothHistoManager * fHistMan; // Histogram Manager\r
+ AliSpectraBothTrackCuts * fTrackCuts; // Track Cuts\r
+ AliSpectraBothEventCuts * fEventCuts; // Event Cuts\r
+ AliSpectraBothPID * fPID;// PID class\r
+ Bool_t fIsMC;// true if processing MC\r
+ Int_t fNRebin; //rebin of histos\r
+ Bool_t fUseMinSigma; // if true use min sigma \r
+ AliESDtrackCuts *fCuts; // ESD track cuts \r
+ Bool_t fdotheMCLoopAfterEventCuts; // if true first check the ESD event cuts than loop over MC info , if flase other approach \r
+\r
+ AliAnalysisTaskSpectraBoth(const AliAnalysisTaskSpectraBoth&);\r
+ AliAnalysisTaskSpectraBoth& operator=(const AliAnalysisTaskSpectraBoth&);\r
+\r
+ ClassDef(AliAnalysisTaskSpectraBoth, 2);\r
+};\r
+\r
+#endif\r
--- /dev/null
+
+/**************************************************************************
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+// AliSpectraBothEventCuts class
+//-----------------------------------------------------------------
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TLegend.h"
+#include "TH1F.h"
+#include "TH1I.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliAODEvent.h"
+#include "AliESDEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAnalysisTaskESDfilter.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliSpectraBothEventCuts.h"
+#include "AliSpectraBothTrackCuts.h"
+//#include "AliSpectraBothHistoManager.h"
+#include <iostream>
+
+using namespace std;
+
+ClassImp(AliSpectraBothEventCuts)
+
+AliSpectraBothEventCuts::AliSpectraBothEventCuts(const char *name) : TNamed(name, "AOD Event Cuts"), fAOD(0),fAODEvent(AliSpectraBothTrackCuts::kAODobject), fTrackBits(0),fIsMC(0),fCentFromV0(0), fUseCentPatchAOD049(0), fUseSDDPatchforLHC11a(kDoNotCheckforSDD),fTrackCuts(0),
+fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fQVectorCutMin(0), fQVectorCutMax(0), fVertexCutMin(0), fVertexCutMax(0), fMultiplicityCutMin(0), fMultiplicityCutMax(0),fMaxChi2perNDFforVertex(0),
+fHistoCuts(0),fHistoVtxBefSel(0),fHistoVtxAftSel(0),fHistoEtaBefSel(0),fHistoEtaAftSel(0),fHistoNChAftSel(0),fHistoQVector(0)
+,fHistoEP(0)
+{
+ // Constructor
+ fHistoCuts = new TH1I("fEventCuts", "Event Cuts", kNVtxCuts, -0.5, kNVtxCuts - 0.5);
+ fHistoVtxBefSel = new TH1F("fHistoVtxBefSel", "Vtx distr before event selection",500,-15,15);
+ fHistoVtxAftSel = new TH1F("fHistoVtxAftSel", "Vtx distr after event selection",500,-15,15);
+ fHistoEtaBefSel = new TH1F("fHistoEtaBefSel", "Eta distr before event selection",500,-2,2);
+ fHistoEtaAftSel = new TH1F("fHistoEtaAftSel", "Eta distr after event selection",500,-2,2);
+ fHistoNChAftSel = new TH1F("fHistoNChAftSel", "NCh distr after event selection",2000,-0.5,1999.5);
+ //fHistoQVectorPos = new TH1F("fHistoQVectorPos", "QVectorPos distribution",100,0,10);
+ //fHistoQVectorNeg = new TH1F("fHistoQVectorNeg", "QVectorNeg distribution",100,0,10);
+ fHistoQVector = new TH1F("fHistoQVector", "QVector with VZERO distribution",100,0,10);
+ fHistoEP = new TH1F("fHistoEP", "EP with VZERO distribution",100,-10,10);
+ fCentralityCutMin = 0.0; // default value of centrality cut minimum, 0 ~ no cut
+ fCentralityCutMax = 10000.0; // default value of centrality cut maximum, ~ no cut
+ // fQVectorPosCutMin=0.0;
+ // fQVectorPosCutMax=10000.0;
+ // fQVectorNegCutMin=0.0;
+ // fQVectorNegCutMax=10000.0;
+ fQVectorCutMin=0.0;
+ fQVectorCutMax=10000.0;
+ fVertexCutMin=-10.0;
+ fVertexCutMax=10.0;
+ fMultiplicityCutMin=-1.0;
+ fMultiplicityCutMax=-1.0;
+ fTrackBits=1;
+ fCentFromV0=kFALSE;
+ fMaxChi2perNDFforVertex=-1;
+}
+
+//______________________________________________________
+Bool_t AliSpectraBothEventCuts::IsSelected(AliVEvent * aod,AliSpectraBothTrackCuts *trackcuts)
+{
+ // Returns true if Event Cuts are selected and applied
+ fAOD = aod;
+ fTrackCuts = trackcuts;
+ fHistoCuts->Fill(kProcessedEvents);
+ Bool_t IsPhysSel = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);//FIXME we can add the trigger mask here
+ if(!IsPhysSel)return IsPhysSel;
+ //loop on tracks, before event selection, filling QA histos
+ AliESDEvent* esdevent=0x0;
+ AliAODEvent* aodevent=0x0;
+ Bool_t isSDD=kFALSE;
+ TString nameoftrack(fAOD->ClassName());
+ if(!nameoftrack.CompareTo("AliESDEvent"))
+ {
+ fAODEvent=AliSpectraBothTrackCuts::kESDobject;
+
+ if(fUseSDDPatchforLHC11a!=kDoNotCheckforSDD)
+ {
+ esdevent=dynamic_cast<AliESDEvent*>(fAOD);
+ if(esdevent->GetFiredTriggerClasses().Contains("ALLNOTRD"))
+ isSDD=kTRUE;
+ }
+ }
+ else if(!nameoftrack.CompareTo("AliAODEvent"))
+ {
+ fAODEvent=AliSpectraBothTrackCuts::kAODobject;
+ if(fUseSDDPatchforLHC11a!=kDoNotCheckforSDD)
+ {
+ aodevent=dynamic_cast<AliAODEvent*>(fAOD);
+ if(aodevent->GetFiredTriggerClasses().Contains("ALLNOTRD"))
+ isSDD=kTRUE;
+ }
+ }
+ else
+ return false;
+ if(fUseSDDPatchforLHC11a==kwithSDD&&isSDD==kFALSE)
+ return false;
+ if(fUseSDDPatchforLHC11a==kwithoutSDD&&isSDD==kTRUE)
+ return false;
+
+
+
+ fHistoCuts->Fill(kPhysSelEvents);
+
+
+
+
+ const AliVVertex * vertex = fAOD->GetPrimaryVertex();//FIXME vertex is recreated
+ if(vertex)fHistoVtxBefSel->Fill(vertex->GetZ());
+ fIsSelected =kFALSE;
+ if(CheckVtxRange() && CheckCentralityCut() && CheckMultiplicityCut() && CheckVtxChi2perNDF()){ //selection on vertex and Centrality
+ fIsSelected=CheckQVectorCut(); // QVector is calculated only if the centrality and vertex are correct (performance)
+ }
+ if(fIsSelected){
+ fHistoCuts->Fill(kAcceptedEvents);
+ if(vertex)
+ fHistoVtxAftSel->Fill(vertex->GetZ());
+ }
+ Int_t Nch=0;
+ for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
+ AliVTrack* track =dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));
+ /* if(fAODEvent==AliSpectraBothTrackCuts::kESDobject)
+ track=dynamic_cast<AliVTrack*>(esdevent->GetTrack(iTracks));
+ else if (fAODEvent==AliSpectraBothTrackCuts::kAODobject)
+ track=dynamic_cast<AliVTrack*>(aodevent->GetTrack(iTracks));
+ else return false;*/
+
+ if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
+ fHistoEtaBefSel->Fill(track->Eta());
+ if(fIsSelected){
+ fHistoEtaAftSel->Fill(track->Eta());
+ Nch++;
+ }
+ }
+ //Printf("NCHARGED_EvSel : %d",Nch);
+ if(fIsSelected)fHistoNChAftSel->Fill(Nch);
+ return fIsSelected;
+}
+
+//______________________________________________________
+Bool_t AliSpectraBothEventCuts::CheckVtxRange()
+{
+ // reject events outside of range
+ const AliVVertex * vertex = fAOD->GetPrimaryVertex();
+ //when moving to 2011 wìone has to add a cut using SPD vertex.
+ //The point is that for events with |z|>20 the vertexer tracks is not working (only 2011!). One has to put a safety cut using SPD vertex large e.g. 15cm
+ if (!vertex)
+ {
+ fHistoCuts->Fill(kVtxNoEvent);
+ return kFALSE;
+ }
+ if(vertex->GetNContributors()<1)
+ {
+ fHistoCuts->Fill(kVtxNoEvent);
+ return kFALSE;
+
+ }
+ TString tmp(vertex->GetTitle());
+ if(tmp.Contains("TPC"))
+ {
+ fHistoCuts->Fill(kVtxNoEvent);
+ return kFALSE;
+ }
+ if (vertex->GetZ() > fVertexCutMin && vertex->GetZ() < fVertexCutMax)
+ {
+ return kTRUE;
+ }
+ fHistoCuts->Fill(kVtxRange);
+ return kFALSE;
+}
+
+//______________________________________________________
+Bool_t AliSpectraBothEventCuts::CheckCentralityCut()
+{
+ // Check centrality cut
+ if ( fCentralityCutMax<0.0 && fCentralityCutMin<0.0 ) return kTRUE;
+ Double_t cent=0;
+ TString CentEstimator="TRK";
+ if(fCentFromV0)CentEstimator="V0M";
+ if(!fUseCentPatchAOD049)cent=fAOD->GetCentrality()->GetCentralityPercentile(CentEstimator.Data());
+ else cent=ApplyCentralityPatchAOD049();
+
+ if ( (cent <= fCentralityCutMax) && (cent >= fCentralityCutMin) ) return kTRUE;
+ fHistoCuts->Fill(kVtxCentral);
+
+ return kFALSE;
+}
+
+//______________________________________________________
+Bool_t AliSpectraBothEventCuts::CheckMultiplicityCut()
+{
+ // Check multiplicity cut
+if(fMultiplicityCutMin<0.0 && fMultiplicityCutMax<0.0)
+ return kTRUE;
+ Int_t Ncharged=0;
+ for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++){
+ AliVTrack* track = dynamic_cast<AliVTrack*>(fAOD->GetTrack(iTracks));
+
+ if (!fTrackCuts->IsSelected(track,kFALSE)) continue;
+
+ Ncharged++;
+ }
+ //Printf("NCHARGED_cut : %d",Ncharged);
+ if(Ncharged>fMultiplicityCutMin && Ncharged<fMultiplicityCutMax)return kTRUE;
+
+ return kFALSE;
+}
+
+//______________________________________________________
+Bool_t AliSpectraBothEventCuts::CheckQVectorCut()
+{
+ if(fQVectorCutMin<0.0 && fQVectorCutMax<0.0)return kTRUE;
+ // Check qvector cut
+ /// FIXME: Q vector
+ // //Selection on QVector, before ANY other selection on the event
+ // //Spectra MUST be normalized wrt events AFTER the selection on Qvector
+ // Double_t Qx2EtaPos = 0, Qy2EtaPos = 0;
+ // Double_t Qx2EtaNeg = 0, Qy2EtaNeg = 0;
+
+ // Int_t multPos = 0;
+ // Int_t multNeg = 0;
+ // for(Int_t iT = 0; iT < fAOD->GetNumberOfTracks(); iT++) {
+ // AliAODTrack* aodTrack = fAOD->GetTrack(iT);
+ // if (!aodTrack->TestFilterBit(fTrackBits)) continue;
+ // if (aodTrack->Eta() >= 0){
+ // multPos++;
+ // Qx2EtaPos += TMath::Cos(2*aodTrack->Phi());
+ // Qy2EtaPos += TMath::Sin(2*aodTrack->Phi());
+ // } else {
+ // multNeg++;
+ // Qx2EtaNeg += TMath::Cos(2*aodTrack->Phi());
+ // Qy2EtaNeg += TMath::Sin(2*aodTrack->Phi());
+ // }
+ // }
+ // Double_t qPos=-999;
+ // if(multPos!=0)qPos= TMath::Sqrt((Qx2EtaPos*Qx2EtaPos + Qy2EtaPos*Qy2EtaPos)/multPos);
+ // Double_t qNeg=-999;
+ // if(multNeg!=0)qNeg= TMath::Sqrt((Qx2EtaNeg*Qx2EtaNeg + Qy2EtaNeg*Qy2EtaNeg)/multNeg);
+ //if(qPos<fQVectorPosCutMin || qPos>fQVectorPosCutMax || qNeg<fQVectorNegCutMin || qNeg>fQVectorNegCutMax)return kFALSE;
+
+ Double_t qxEPVZERO = 0, qyEPVZERO = 0;
+ Double_t qVZERO = -999;
+ Double_t psi=fAOD->GetEventplane()->CalculateVZEROEventPlane(fAOD,10,2,qxEPVZERO,qyEPVZERO);
+
+ qVZERO= TMath::Sqrt(qxEPVZERO*qxEPVZERO + qyEPVZERO*qyEPVZERO);
+ if(qVZERO<fQVectorCutMin || qVZERO>fQVectorCutMax)return kFALSE;
+ fHistoQVector->Fill(qVZERO);
+ fHistoEP->Fill(psi);
+
+ fHistoCuts->Fill(kQVector);
+ // fHistoQVectorPos->Fill(qPos);
+ // fHistoQVectorNeg->Fill(qNeg);
+ return kTRUE;
+}
+//____________________________________________________________
+ Bool_t AliSpectraBothEventCuts::CheckVtxChi2perNDF()
+ {
+ if(fMaxChi2perNDFforVertex<0)
+ return kTRUE;
+ const AliVVertex * vertex = fAOD->GetPrimaryVertex();
+ if(TMath::Abs(vertex->GetChi2perNDF())>fMaxChi2perNDFforVertex)
+ return kFALSE;
+ return kTRUE;
+ }
+
+
+
+//______________________________________________________
+void AliSpectraBothEventCuts::PrintCuts()
+{
+ // print info about event cuts
+ cout << "Event Stats" << endl;
+ cout << " > Number of accepted events: " << fHistoCuts->GetBinContent(kAcceptedEvents + 1) << endl;
+ cout << " > Number of processed events: " << fHistoCuts->GetBinContent(kProcessedEvents + 1) << endl;
+ cout << " > Number of PhysSel events: " << fHistoCuts->GetBinContent(kPhysSelEvents + 1) << endl;
+ cout << " > Vertex out of range: " << fHistoCuts->GetBinContent(kVtxRange + 1) << endl;
+ cout << " > Events cut by centrality: " << fHistoCuts->GetBinContent(kVtxCentral + 1) << endl;
+ cout << " > Events without vertex: " << fHistoCuts->GetBinContent(kVtxNoEvent + 1) << endl;
+ cout << " > QVector cut: " << fHistoCuts->GetBinContent(kQVector + 1) << endl;
+ cout << " > Track type used for the QVector calculation: " << fTrackBits << endl;
+ // cout << " > QPosRange: [" << fQVectorPosCutMin <<"," <<fQVectorPosCutMax<<"]"<< endl;
+ // cout << " > QNegRange: [" << fQVectorNegCutMin <<"," <<fQVectorNegCutMax<<"]"<< endl;
+ cout << " > QRange: [" << fQVectorCutMin <<"," <<fQVectorCutMax<<"]"<< endl;
+ cout << " > Vertex: [" << fVertexCutMin <<"," <<fVertexCutMax<<"]"<< endl;
+ cout << " > Multiplicity: [" << fMultiplicityCutMin <<"," <<fMultiplicityCutMax<<"]"<< endl;
+ cout << " > Centrality: [" << fCentralityCutMin <<"," <<fCentralityCutMax<<"]"<< endl;
+}
+//______________________________________________________
+
+Double_t AliSpectraBothEventCuts::ApplyCentralityPatchAOD049()
+{
+ //
+ //Apply centrality patch for AOD049 outliers
+ //
+ // if (fCentralityType!="V0M") {
+ // AliWarning("Requested patch forAOD049 for wrong value (not centrality from V0).");
+ // return -999.0;
+ // }
+ AliCentrality *centrality = fAOD->GetCentrality();
+ if (!centrality) {
+ AliWarning("Cannot get centrality from AOD event.");
+ return -999.0;
+ }
+
+ Float_t cent = (Float_t)(centrality->GetCentralityPercentile("V0M"));
+ /*
+ Bool_t isSelRun = kFALSE;
+ Int_t selRun[5] = {138364, 138826, 138828, 138836, 138871};
+ if(cent<0){
+ Int_t quality = centrality->GetQuality();
+ if(quality<=1){
+ cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
+ } else {
+ Int_t runnum=aodEvent->GetRunNumber();
+ for(Int_t ir=0;ir<5;ir++){
+ if(runnum==selRun[ir]){
+ isSelRun=kTRUE;
+ break;
+ }
+ }
+ if((quality==8||quality==9)&&isSelRun) cent=(Float_t)centrality->GetCentralityPercentileUnchecked("V0M");
+ }
+ }
+ */
+ if(cent>=0.0) {
+ Float_t v0 = 0.0;
+ AliAODEvent *aodEvent = (AliAODEvent *)fAOD;
+ AliAODVZERO *aodV0 = (AliAODVZERO *) aodEvent->GetVZEROData();
+ v0+=aodV0->GetMTotV0A();
+ v0+=aodV0->GetMTotV0C();
+ if ( (cent==0) && (v0<19500) ) {
+ AliDebug(3, Form("Filtering issue in centrality -> cent = %5.2f",cent));
+ return -999.0;
+ }
+ Float_t tkl = (Float_t)(aodEvent->GetTracklets()->GetNumberOfTracklets());
+ Float_t val = 1.30552 + 0.147931 * v0;
+
+ Float_t tklSigma[101] = {176.644, 156.401, 153.789, 153.015, 142.476, 137.951, 136.127, 129.852, 127.436, 124.86,
+ 120.788, 115.611, 113.172, 110.496, 109.127, 104.421, 102.479, 99.9766, 97.5152, 94.0654,
+ 92.4602, 89.3364, 87.1342, 83.3497, 82.6216, 81.1084, 78.0793, 76.1234, 72.9434, 72.1334,
+ 68.0056, 68.2755, 66.0376, 62.9666, 62.4274, 59.65, 58.3776, 56.6361, 54.5184, 53.4224,
+ 51.932, 50.8922, 48.2848, 47.912, 46.5717, 43.4114, 43.2083, 41.3065, 40.1863, 38.5255,
+ 37.2851, 37.5396, 34.4949, 33.8366, 31.8043, 31.7412, 30.8392, 30.0274, 28.8793, 27.6398,
+ 26.6488, 25.0183, 25.1489, 24.4185, 22.9107, 21.2002, 21.6977, 20.1242, 20.4963, 19.0235,
+ 19.298, 17.4103, 16.868, 15.2939, 15.2939, 16.0295, 14.186, 14.186, 15.2173, 12.9504, 12.9504,
+ 12.9504, 15.264, 12.3674, 12.3674, 12.3674, 12.3674, 12.3674, 18.3811, 13.7544, 13.7544,
+ 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544, 13.7544
+ };
+
+ if ( TMath::Abs(tkl-val) > 6.*tklSigma[(Int_t)cent] ) {
+ AliDebug(3, Form("Outlier event in centrality -> cent = %5.2f",cent));
+ return -999.0;
+ }
+ } else {
+ //force it to be -999. whatever the negative value was
+ cent = -999.;
+ }
+ return cent;
+}
+
+//______________________________________________________
+
+
+Long64_t AliSpectraBothEventCuts::Merge(TCollection* list)
+{
+ // Merge a list of AliSpectraBothEventCuts objects with this.
+ // Returns the number of merged objects (including this).
+
+ // AliInfo("Merging");
+
+ if (!list)
+ return 0;
+
+ if (list->IsEmpty())
+ return 1;
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj;
+
+ // collections of all histograms
+ TList collections;//FIXME we should use only 1 collection
+ TList collections_histoVtxBefSel;
+ TList collections_histoVtxAftSel;
+ TList collections_histoEtaBefSel;
+ TList collections_histoEtaAftSel;
+ TList collections_histoNChAftSel;
+ // TList collections_histoQVectorPos;
+ // TList collections_histoQVectorNeg;
+ TList collections_histoQVector;
+ TList collections_histoEP;
+
+ Int_t count = 0;
+
+ while ((obj = iter->Next())) {
+ AliSpectraBothEventCuts* entry = dynamic_cast<AliSpectraBothEventCuts*> (obj);
+ if (entry == 0)
+ continue;
+
+ TH1I * histo = entry->GetHistoCuts();
+ collections.Add(histo);
+ TH1F * histo_histoVtxBefSel = entry->GetHistoVtxBefSel();
+ collections_histoVtxBefSel.Add(histo_histoVtxBefSel);
+ TH1F * histo_histoVtxAftSel = entry->GetHistoVtxAftSel();
+ collections_histoVtxAftSel.Add(histo_histoVtxAftSel);
+ TH1F * histo_histoEtaBefSel = entry->GetHistoEtaBefSel();
+ collections_histoEtaBefSel.Add(histo_histoEtaBefSel);
+ TH1F * histo_histoEtaAftSel = entry->GetHistoEtaAftSel();
+ collections_histoEtaAftSel.Add(histo_histoEtaAftSel);
+ TH1F * histo_histoNChAftSel = entry->GetHistoNChAftSel();
+ collections_histoNChAftSel.Add(histo_histoNChAftSel);
+ // TH1F * histo_histoQVectorPos = entry->GetHistoQVectorPos();
+ // collections_histoQVectorPos.Add(histo_histoQVectorPos);
+ // TH1F * histo_histoQVectorNeg = entry->GetHistoQVectorNeg();
+ // collections_histoQVectorNeg.Add(histo_histoQVectorNeg);
+ TH1F * histo_histoQVector = entry->GetHistoQVector();
+ collections_histoQVector.Add(histo_histoQVector);
+ TH1F * histo_histoEP = entry->GetHistoEP();
+ collections_histoEP.Add(histo_histoEP);
+ count++;
+ }
+
+ fHistoCuts->Merge(&collections);
+ fHistoVtxBefSel->Merge(&collections_histoVtxBefSel);
+ fHistoVtxAftSel->Merge(&collections_histoVtxAftSel);
+ fHistoEtaBefSel->Merge(&collections_histoEtaBefSel);
+ fHistoEtaAftSel->Merge(&collections_histoEtaAftSel);
+ fHistoNChAftSel->Merge(&collections_histoNChAftSel);
+ // fHistoQVectorPos->Merge(&collections_histoQVectorPos);
+ // fHistoQVectorNeg->Merge(&collections_histoQVectorNeg);
+ fHistoQVector->Merge(&collections_histoQVector);
+ fHistoEP->Merge(&collections_histoEP);
+
+
+ delete iter;
+
+ return count+1;
+}
+
--- /dev/null
+#ifndef ALISPECTRABOTHEVENTCUTS_H
+#define ALISPECTRABOTHEVENTCUTS_H
+
+/* See cxx source for full Copyright notice */
+
+//-------------------------------------------------------------------------
+// AliSpectraBothEventCuts
+//
+//
+//
+//
+// Authors: Michele Floris, CERN, Philip Versteeg, UU, Redmer Bertens, UU
+//-------------------------------------------------------------------------
+
+class AliVEvent;
+class AliSpectraBothTrackCuts;
+//class AliSpectraBothHistoManager;
+
+#include "TH1I.h"
+#include "TNamed.h"
+#include "AliSpectraBothTrackCuts.h"
+class AliSpectraBothEventCuts : public TNamed
+{
+ public:
+ enum { kProcessedEvents = 0,kPhysSelEvents,kAcceptedEvents, kVtxRange, kVtxCentral, kVtxNoEvent, kQVector, kNVtxCuts};
+enum {kDoNotCheckforSDD=0,kwithSDD,kwithoutSDD};
+
+ // Constructors
+ AliSpectraBothEventCuts() : TNamed(), fAOD(0),fAODEvent(AliSpectraBothTrackCuts::kAODobject),fTrackBits(0), fIsMC(0), fCentFromV0(0), fUseCentPatchAOD049(0),fUseSDDPatchforLHC11a(kDoNotCheckforSDD),fTrackCuts(0),
+fIsSelected(0), fCentralityCutMin(0), fCentralityCutMax(0), fQVectorCutMin(0), fQVectorCutMax(0), fVertexCutMin(0), fVertexCutMax(0), fMultiplicityCutMin(0), fMultiplicityCutMax(0), fMaxChi2perNDFforVertex(0),fHistoCuts(0),fHistoVtxBefSel(0),fHistoVtxAftSel(0),fHistoEtaBefSel(0),fHistoEtaAftSel(0),fHistoNChAftSel(0),fHistoQVector(0),fHistoEP(0) {}
+ AliSpectraBothEventCuts(const char *name);
+ virtual ~AliSpectraBothEventCuts() {}
+
+ void SetIsMC(Bool_t isMC = kFALSE) {fIsMC = isMC; };
+ Bool_t GetIsMC() const { return fIsMC;};
+ void SetCentFromV0(Bool_t centFromV0 = kFALSE) {fCentFromV0 = centFromV0; };
+ Bool_t GetCentFromV0() const { return fCentFromV0;};
+
+ void SetUseCentPatchAOD049(Bool_t useCentPatchAOD049 = kFALSE) {fUseCentPatchAOD049 = useCentPatchAOD049; };
+ Bool_t GetUseCentPatchAOD049() const { return fUseCentPatchAOD049;};
+ void SetUseSDDPatchforLHC11a(Int_t useSDDPatchforLHC11a) {fUseSDDPatchforLHC11a=useSDDPatchforLHC11a;} ;
+ Int_t GetUseSDDPatchforLHC11a() {return fUseSDDPatchforLHC11a;};
+
+ // Methods
+ Bool_t IsSelected(AliVEvent * aod,AliSpectraBothTrackCuts *trackcuts);
+ Bool_t CheckVtxRange();
+ Bool_t CheckCentralityCut();
+ Bool_t CheckMultiplicityCut();
+ Bool_t CheckQVectorCut();
+ Bool_t CheckVtxChi2perNDF();
+ void SetCentralityCutMin(Float_t cut) { fCentralityCutMin = cut; }
+ void SetCentralityCutMax(Float_t cut) { fCentralityCutMax = cut; }
+ //void SetQVectorPosCut(Float_t min,Float_t max) { fQVectorPosCutMin = min; fQVectorPosCutMax = max; }
+ //void SetQVectorNegCut(Float_t min,Float_t max) { fQVectorNegCutMin = min; fQVectorNegCutMax = max; }
+ void SetQVectorCut(Float_t min,Float_t max) { fQVectorCutMin = min; fQVectorCutMax = max; }
+ void SetVertexCut(Float_t min,Float_t max) { fVertexCutMin = min; fVertexCutMax = max; }
+ void SetMultiplicityCut(Float_t min,Float_t max) { fMultiplicityCutMin = min; fMultiplicityCutMax = max; }
+ void SetTrackBits(UInt_t TrackBits) {fTrackBits=TrackBits;}
+ void SetMaxChi2perNDFforVertex(Float_t cut) { fMaxChi2perNDFforVertex=cut;}
+
+ UInt_t GetTrackType() const { return fTrackBits;}
+ TH1I * GetHistoCuts() { return fHistoCuts; }
+ TH1F * GetHistoVtxBefSel() { return fHistoVtxBefSel; }
+ TH1F * GetHistoVtxAftSel() { return fHistoVtxAftSel; }
+ TH1F * GetHistoEtaBefSel() { return fHistoEtaBefSel; }
+ TH1F * GetHistoEtaAftSel() { return fHistoEtaAftSel; }
+ TH1F * GetHistoNChAftSel() { return fHistoNChAftSel; }
+ /* TH1F * GetHistoQVectorPos() { return fHistoQVectorPos; } */
+ /* TH1F * GetHistoQVectorNeg() { return fHistoQVectorNeg; } */
+ TH1F * GetHistoQVector() { return fHistoQVector; }
+ TH1F * GetHistoEP() { return fHistoEP; }
+ Float_t GetCentralityMin() const { return fCentralityCutMin; }
+ Float_t GetCentralityMax() const { return fCentralityCutMax; }
+ //Float_t GetQVectorPosCutMin() const { return fQVectorPosCutMin; }
+ //Float_t GetQVectorPosCutMax() const { return fQVectorPosCutMax; }
+ //Float_t GetQVectorNegCutMin() const { return fQVectorNegCutMin; }
+ //Float_t GetQVectorNegCutMax() const { return fQVectorNegCutMax; }
+ Float_t GetQVectorCutMin() const { return fQVectorCutMin; }
+ Float_t GetQVectorCutMax() const { return fQVectorCutMax; }
+ Float_t GetVertexCutMin() const { return fVertexCutMin; }
+ Float_t GetVertexCutMax() const { return fVertexCutMax; }
+ Float_t GetMultiplicityCutMin() const { return fMultiplicityCutMin; }
+ Float_t GetMultiplicityCutMax() const { return fMultiplicityCutMax; }
+ Float_t GetMaxChi2perNDFforVertex() const {return fMaxChi2perNDFforVertex;}
+
+ void PrintCuts();
+ Double_t ApplyCentralityPatchAOD049();
+
+ Float_t NumberOfEvents() { return fHistoCuts->GetBinContent(kAcceptedEvents+1); }
+ Float_t NumberOfProcessedEvents() { return fHistoCuts->GetBinContent(kProcessedEvents+1); }
+ Float_t NumberOfPhysSelEvents() { return fHistoCuts->GetBinContent(kPhysSelEvents+1); }
+
+ Long64_t Merge(TCollection* list);
+
+
+ private:
+
+ AliVEvent *fAOD; //! AOD event
+ Int_t fAODEvent; // flag what type event is conected use values form AliSpeatraAODTrackCuts
+ UInt_t fTrackBits; // Type of track to be used in the Qvector calculation
+ Bool_t fIsMC;// true if processing MC
+ Bool_t fCentFromV0;// default centrality with tracks
+ Bool_t fUseCentPatchAOD049;// Patch for centrality selection on AOD049
+ Int_t fUseSDDPatchforLHC11a; // if true will check for ALLNOTRD in fired trigger class
+
+ AliSpectraBothTrackCuts *fTrackCuts; //! track cuts
+ Bool_t fIsSelected; // True if cuts are selected
+ Float_t fCentralityCutMin; // minimum centrality percentile
+ Float_t fCentralityCutMax; // maximum centrality percentile
+ /* Float_t fQVectorPosCutMin; // minimum qvecPos */
+ /* Float_t fQVectorPosCutMax; // maximum qvecPos */
+ /* Float_t fQVectorNegCutMin; // minimum qvecNeg */
+ /* Float_t fQVectorNegCutMax; // maximum qvecNeg */
+ Float_t fQVectorCutMin; // minimum qvecPos
+ Float_t fQVectorCutMax; // maximum qvecPos
+ Float_t fVertexCutMin; // minimum vertex position
+ Float_t fVertexCutMax; // maximum vertex position
+ Float_t fMultiplicityCutMin; // minimum multiplicity position
+ Float_t fMultiplicityCutMax; // maximum multiplicity position
+ Float_t fMaxChi2perNDFforVertex; // maximum value of Chi2perNDF of vertex
+ TH1I *fHistoCuts; // Cuts statistics
+ TH1F *fHistoVtxBefSel; // Vtx distr before event selection
+ TH1F *fHistoVtxAftSel; // Vtx distr after event selection
+ TH1F *fHistoEtaBefSel; // Eta distr before event selection
+ TH1F *fHistoEtaAftSel; // Eta distr after event selection
+ TH1F *fHistoNChAftSel; // NCh distr after event selection
+ //TH1F *fHistoQVectorPos; // QVectorPos
+ //TH1F *fHistoQVectorNeg; // QVectorNeg
+ TH1F *fHistoQVector; // QVector
+ TH1F *fHistoEP; // EP
+ AliSpectraBothEventCuts(const AliSpectraBothEventCuts&);
+ AliSpectraBothEventCuts& operator=(const AliSpectraBothEventCuts&);
+
+ ClassDef(AliSpectraBothEventCuts, 3);
+
+};
+#endif
+
--- /dev/null
+
+/**************************************************************************
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+// AliSpectraBothHistoManager class
+//-----------------------------------------------------------------
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TLegend.h"
+#include "TH1.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAnalysisTaskESDfilter.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliSpectraBothHistoManager.h"
+#include <iostream>
+
+using namespace std;
+
+ClassImp(AliSpectraBothHistoManager)
+
+
+using namespace AliSpectraNameSpaceBoth;
+#include "HistogramNamesBoth.cxx" // generate this automatically running createNames.py
+
+
+//const char* kParticleSpecies[] =
+ // {
+ // "PionPlus",
+ // "KaonPlus",
+ // "ProtonPlus",
+ // "PionMinus",
+ // "KaonMinus",
+ // "ProtonMinus",
+ // };
+
+
+AliSpectraBothHistoManager::AliSpectraBothHistoManager(const char *name,Int_t nrebin): TNamed(name, "AOD Spectra Histo Manager"), fOutputList(0), fNRebin(0)
+{
+ // ctor
+ fNRebin=nrebin;
+ fOutputList = new TList;
+ fOutputList->SetOwner();
+ Bool_t oldStatus = TH1::AddDirectoryStatus();
+ TH1::AddDirectory(kFALSE);
+ for (Int_t ihist = 0; ihist < kNHist ; ihist++)
+ {
+ if (ihist <= kNPtGenHist) BookPtGenHistogram(kHistNameBoth[ihist]); // PT histos
+ if (ihist > kNPtGenHist && ihist <= kNPtGenAllChHist) BookPtGenAllChHistogram(kHistNameBoth[ihist]); // PT histos
+ if (ihist > kNPtGenAllChHist && ihist <= kNPtRecHist) BookPtRecHistogram(kHistNameBoth[ihist]); // PT histos
+ if (ihist > kNPtRecHist && ihist <= kNPtRecAllChHist) BookPtRecAllChHistogram(kHistNameBoth[ihist]); // PT histos
+ if (ihist > kNPtRecAllChHist && ihist <= kNHistPID) BookPIDHistogram(kHistNameBoth[ihist]); // PID histos
+ if (ihist > kNHistPID && ihist <= kNHistNSig) BookNSigHistogram(kHistNameBoth[ihist]); // NSigmaSep histos
+ if (ihist > kNHistNSig && ihist<kHistGenMulvsRawMul) BookqVecHistogram(kHistNameBoth[ihist]); // qDistr histos
+ if(ihist==kHistGenMulvsRawMul) BookGenMulvsRawMulHistogram(kHistNameBoth[ihist]);
+ }
+
+ TH1::AddDirectory(oldStatus);
+}
+
+//_______________________________________________________
+
+TH2F* AliSpectraBothHistoManager::BookPtGenHistogram(const char * name)
+{
+ // Return a pt histogram with predefined binning, set the ID and add it to the output list
+ AliInfo(Form("Booking pt gen histogram %s", name));
+
+ //standard histo
+ const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
+ Int_t nbinsTempl=52;
+
+ TH2F * hist = new TH2F(name,Form("P_{T} distribution (%s)", name),nbinsTempl,templBins,2,-0.5,1.5);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
+ hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
+ hist->GetYaxis()->SetTitle("IsPhysicalPrimary()");
+ hist->SetMarkerStyle(kFullCircle);
+ hist->Sumw2();
+ fOutputList->Add(hist);
+
+ return hist;
+}
+
+//_______________________________________________________
+
+TH2F* AliSpectraBothHistoManager::BookPtGenAllChHistogram(const char * name)
+{
+ // Return a pt histogram with predefined binning, set the ID and add it to the output list
+ AliInfo(Form("Booking pt gen histogram - no PID %s", name));
+
+ //standard histo
+ const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0,5.5,6.0,6.5,7,7.5,8,8.5,9,9.5,10};
+ Int_t nbinsTempl=62;
+
+ TH2F * hist = new TH2F(name,Form("P_{T} distribution (All Ch) (%s)", name),nbinsTempl,templBins,2,-0.5,1.5);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
+ hist->GetXaxis()->SetTitle("generated P_{T} (GeV / c)");
+ hist->GetYaxis()->SetTitle("IsPhysicalPrimary()");
+ hist->SetMarkerStyle(kFullCircle);
+ hist->Sumw2();
+ fOutputList->Add(hist);
+
+ return hist;
+}
+
+
+//_______________________________________________________
+TH2F* AliSpectraBothHistoManager::BookPtRecHistogram(const char * name)
+{
+ // Return a pt histogram with predefined binning, set the ID and add it to the output list
+ AliInfo(Form("Booking pt rec histogram %s, rebin:%d", name, fNRebin));
+
+ //standard histo
+ const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
+ Int_t nbinsTempl=52;
+
+ TH2F * hist = new TH2F(name,Form("reconstructed P_{T} distribution (%s)", name),nbinsTempl,templBins,3000,-3,3);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
+ if(fNRebin!=0)hist->RebinY(fNRebin);
+ hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
+ hist->GetYaxis()->SetTitle("DCA xy");
+ hist->SetMarkerStyle(kFullCircle);
+ hist->Sumw2();
+ fOutputList->Add(hist);
+
+ return hist;
+}
+
+//_______________________________________________________
+TH2F* AliSpectraBothHistoManager::BookPtRecAllChHistogram(const char * name)
+{
+ // Return a pt histogram with predefined binning, set the ID and add it to the output list
+ AliInfo(Form("Booking pt rec histogram %s, rebin:%d", name, fNRebin));
+
+ //standard histo
+ const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0,5.5,6.0,6.5,7,7.5,8,8.5,9,9.5,10};
+ Int_t nbinsTempl=62;
+
+ TH2F * hist = new TH2F(name,Form("reconstructed P_{T} distribution (All Ch) (%s)", name),nbinsTempl,templBins,3000,-3,3);//need to be at least 1 becuase the generated are filled with (pt,IsPhysPrim)
+ if(fNRebin!=0)hist->RebinY(fNRebin);
+ hist->GetXaxis()->SetTitle("P_{T} (GeV / c)");
+ hist->GetYaxis()->SetTitle("DCA xy");
+ hist->SetMarkerStyle(kFullCircle);
+ hist->Sumw2();
+ fOutputList->Add(hist);
+
+ return hist;
+}
+
+//_____________________________________________________________________________
+
+TH2F* AliSpectraBothHistoManager::BookPIDHistogram(const char * name)
+{
+ // Return a pt histogram with predefined binning, set the ID and add it to the output list
+ AliInfo(Form("Booking PID histogram %s, rebin:%d", name, fNRebin));
+TString tmp(name);
+ TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 50, 0, 2.5, 200, -1000, 1000);
+ if(fNRebin!=0){
+ hist->RebinY(fNRebin);
+ // hist->RebinX(fNRebin);
+ }
+ if(tmp.Contains("Pt"))
+ hist->GetXaxis()->SetTitle("P_{t} (GeV / c)");
+ else
+ hist->GetXaxis()->SetTitle("P (GeV / c)");
+
+ hist->GetYaxis()->SetTitle("PID signal");
+ // hist->Sumw2();
+ fOutputList->Add(hist);
+
+ return hist;
+}
+
+//_____________________________________________________________________________
+
+TH2F* AliSpectraBothHistoManager::BookNSigHistogram(const char * name)
+{
+ // Return a pt histogram with predefined binning, set the ID and add it to the output list
+ AliInfo(Form("Booking NSigma histogram %s, rebin:%d", name, fNRebin));
+ Int_t nbins=200;
+ Float_t miny=-20;
+ TString tmp(name);
+ if(tmp.Contains("TPCTOF"))
+ {
+ nbins=100;
+ miny=0.0;
+ }
+ TH2F * hist = new TH2F(name, Form("Particle Identification (%s)", name), 50, 0, 2.5,nbins,miny, 20);
+ if(fNRebin!=0){
+ hist->RebinY(fNRebin);
+ //hist->RebinX(fNRebin);
+ }
+ if(tmp.Contains("Pt"))
+ hist->GetXaxis()->SetTitle("P_{t} (GeV / c)");
+ else
+ hist->GetXaxis()->SetTitle("P (GeV / c)");
+ //hist->GetYaxis()->SetTitle("TPC");
+ //hist->Sumw2();
+ fOutputList->Add(hist);
+
+ return hist;
+}
+
+//_____________________________________________________________________________
+
+TH2F* AliSpectraBothHistoManager::BookqVecHistogram(const char * name)
+{
+ // Return a pt histogram with predefined binning, set the ID and add it to the output list
+ AliInfo(Form("Booking q Vector histogram %s", name));
+
+ TH2F * hist = new TH2F(name, Form("q Vector distribution vs Centrality (%s)", name), 100, 0, 10, 100, 0, 100);
+ hist->GetXaxis()->SetTitle("q vector");
+ hist->GetYaxis()->SetTitle("Centrality (V0)");
+ // hist->Sumw2();
+ fOutputList->Add(hist);
+
+ return hist;
+}
+//____________________________________________________________________________________
+ TH2F* AliSpectraBothHistoManager::BookGenMulvsRawMulHistogram(const char * name)
+{
+ // Return a pt histogram with predefined binning, set the ID and add it to the output list
+ AliInfo(Form("Booking Gen vs Raw multilicity histogram %s", name));
+
+ TH2F * hist = new TH2F(name, Form("Gen vs Raw multilicity (%s)", name), 100, -0.5, 99.5, 100, -0.5, 99.5);
+ hist->GetXaxis()->SetTitle("gen mul ");
+ hist->GetYaxis()->SetTitle("raw mul ");
+ // hist->Sumw2();
+ fOutputList->Add(hist);
+
+ return hist;
+}
+//_____________________________________________________________________________
+
+TH1F* AliSpectraBothHistoManager::GetPtHistogram1D(const char * name,Double_t minDCA,Double_t maxDCA)
+{
+ // //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
+ // //if minDCA=-1 && maxDCA=-1 the projection is done using the full DCA range
+ TH2F *hist=(TH2F*)fOutputList->FindObject(name);
+ TH1F *outhist=0x0;
+ Printf("--- Projecting %s on Xaxis[%f,%f]:",name,minDCA,maxDCA);
+ if(minDCA==-1 && maxDCA==-1){
+ outhist=(TH1F*)hist->ProjectionX("_px",0,-1,"e");
+ Printf("Full Range");
+ }else {
+ Int_t firstbin=hist->GetYaxis()->FindBin(minDCA);
+ Int_t lastbin=hist->GetYaxis()->FindBin(maxDCA);
+ Printf("firstbin: %d lastbin: %d",firstbin,lastbin);
+ outhist=(TH1F*)hist->ProjectionX("_px",firstbin,lastbin,"e");
+ }
+ Printf("Entries outhist: %.0f Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries());
+ return outhist;
+}
+
+//_____________________________________________________________________________
+
+TH1F* AliSpectraBothHistoManager::GetDCAHistogram1D(const char * name,Double_t minPt,Double_t maxPt)
+{
+ // //return the projection of the TH2 (pt,DCA) in the DCA bin range [firstDCAbin,lastDCAbin]
+ // //if minPt=-1 && maxPt=-1 the projection is done using the full DCA range
+ TH2F *hist=(TH2F*)fOutputList->FindObject(name);
+ TH1F *outhist=0x0;
+ Printf("--- Projecting %s on Yaxis[%f,%f]:",name,minPt,maxPt);
+ if(minPt==-1 && maxPt==-1){
+ outhist=(TH1F*)hist->ProjectionY("_py",0,-1,"e");
+ Printf("Full Range");
+ }else {
+ Int_t firstbin=hist->GetXaxis()->FindBin(minPt);
+ Int_t lastbin=hist->GetXaxis()->FindBin(maxPt);
+ Printf("firstbin: %d lastbin: %d",firstbin,lastbin);
+ outhist=(TH1F*)hist->ProjectionY("_py",firstbin,lastbin,"e");
+ Printf("GetDCAHistogram1D(%s) BinRange:%d %d Pt Range: %f %f",hist->GetName(),firstbin,lastbin,hist->GetXaxis()->GetBinLowEdge(firstbin),hist->GetXaxis()->GetBinLowEdge(firstbin)+hist->GetXaxis()->GetBinWidth(lastbin));
+ }
+ Printf("Entries outhist: %.0f Entries hist: %.0f",outhist->GetEntries(),hist->GetEntries());
+ return outhist;
+}
+
+//_____________________________________________________________________________
+
+Long64_t AliSpectraBothHistoManager::Merge(TCollection* list)
+{
+ // Merge a list of AliSpectraBothHistoManager objects with this.
+ // Returns the number of merged objects (including this).
+
+ // AliInfo("Merging");
+
+ if (!list)
+ return 0;
+
+ if (list->IsEmpty())
+ return 1;
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj;
+
+ // collections of all histograms
+ TList collections;
+
+ Int_t count = 0;
+
+ while ((obj = iter->Next())) {
+ AliSpectraBothHistoManager* entry = dynamic_cast<AliSpectraBothHistoManager*> (obj);
+ if (entry == 0)
+ continue;
+
+ TList * hlist = entry->GetOutputList();
+ collections.Add(hlist);
+ count++;
+ }
+
+ fOutputList->Merge(&collections);
+
+ delete iter;
+
+ return count+1;
+}
+
+
+TH1* AliSpectraBothHistoManager::GetHistogram1D(UInt_t histoType, UInt_t particleType, UInt_t charge) {
+ // GetHistogram using particle ID and histogram type
+ Int_t baseId = -1;
+
+ if (particleType == kSpUndefined) {
+ AliError ("Trying to get histo for undefined particle");
+ return 0;
+ }
+
+ switch(histoType) {
+ case kHistPtGenTruePrimary:
+ baseId = kHistPtGenTruePrimaryPionPlus;
+ break;
+ case kHistPtRecSigma:
+ baseId = kHistPtRecSigmaPionPlus;
+ break;
+ case kHistPtRecTrue:
+ baseId = kHistPtRecTruePionPlus;
+ break;
+ case kHistPtRecTruePrimary:
+ baseId = kHistPtRecTruePrimaryPionPlus;
+ break;
+ case kHistPtRecPrimary:
+ baseId = kHistPtRecPrimaryPionPlus;
+ break;
+ case kHistPtRecSigmaPrimary:
+ baseId = kHistPtRecSigmaPrimaryPionPlus;
+ break;
+ case kHistPtRecSigmaSecondaryMaterial:
+ baseId = kHistPtRecSigmaSecondaryMaterialPionPlus;
+ break;
+ case kHistPtRecSigmaSecondaryWeakDecay:
+ baseId = kHistPtRecSigmaSecondaryWeakDecayPionPlus;
+ break;
+ case kHistNSigTPC:
+ baseId = kHistNSigPionTPC;
+ break;
+ case kHistNSigTOF:
+ baseId = kHistNSigPionTOF;
+ break;
+ case kHistNSigTPCTOF:
+ baseId = kHistNSigPionTPCTOF;
+ break;
+ default:
+ baseId = -1;
+ }
+
+ if (baseId < 0)
+ AliFatal(Form("Wrong histogram type %d", histoType));
+
+ //cout << "T[" << histoType << "] ID["<< baseId <<"] P["<<particleType<<"] C[" << charge
+ // << " --> ["<< baseId + particleType + 3*(charge) <<"] = " ;
+
+ baseId = baseId + particleType + 3*(charge);
+
+ //cout << GetHistogram(baseId)->GetName() << endl;
+
+ return GetHistogram(baseId);
+}
+
+TH2* AliSpectraBothHistoManager::GetHistogram2D(UInt_t histoType, UInt_t particleType, UInt_t charge){
+ // returns histo based on ids, casting it to TH2*
+ return (TH2*) GetHistogram1D(histoType,particleType,charge);
+
+
+}
--- /dev/null
+#ifndef ALISPECTRABOTHHISTOMANAGER_H
+#define ALISPECTRABOTHHISTOMANAGER_H
+
+/* See cxx source for full Copyright notice */
+
+//-------------------------------------------------------------------------
+// AliSpectraBothHistoManager
+//
+//
+//
+//
+// Authors: Michele Floris, CERN, Philip Versteeg, UU, Redmer Bertens, UU
+//-------------------------------------------------------------------------
+
+class AliAODEvent;
+class TH1F;
+class TH2F;
+class TH1;
+class TH2;
+
+//class TList;
+#include "TNamed.h"
+#include "TList.h"
+#include "HistogramsBoth.h" // Change this file if you want to add an histogram
+#include "HistogramNamesBoth.h" // generate this automatically running createNames.py
+
+namespace AliSpectraNameSpaceBoth
+{
+
+
+ enum BothParticleSpecies_t
+ {
+ kSpPion,
+ kSpKaon,
+ kSpProton,
+ kNSpecies,
+ kSpUndefined,
+ }; // Particle species used in plotting
+
+ //extern const char * kParticleSpecies[];
+
+ enum BothHistoType_t
+ {
+ kHistPtGenTruePrimary,
+ kHistPtRecSigma,
+ kHistPtRecTrue,
+ kHistPtRecTruePrimary,
+ kHistPtRecSigmaPrimary,
+ kHistPtRecSigmaSecondaryMaterial,
+ kHistPtRecSigmaSecondaryWeakDecay,
+ kHistNSigTPC,
+ kHistNSigTOF,
+ kHistNSigTPCTOF,
+ kNHistoTypes
+ }; // Types of histos
+
+ enum BothCharge_t
+ {
+ kChPos = 0,
+ kChNeg,
+ kNCharge
+ };
+
+}
+
+using namespace AliSpectraNameSpaceBoth;
+
+class AliSpectraBothHistoManager : public TNamed
+{
+public:
+ AliSpectraBothHistoManager() : TNamed(), fOutputList(0), fNRebin(0) {}
+ AliSpectraBothHistoManager(const char *name,Int_t nrebin);
+ virtual ~AliSpectraBothHistoManager() {}
+
+
+ TH2F* BookPtGenHistogram(const char * name);
+ TH2F* BookPtGenAllChHistogram(const char * name);
+ TH2F* BookPtRecHistogram(const char * name);
+ TH2F* BookPtRecAllChHistogram(const char * name);
+ TH2F* BookPIDHistogram(const char * name);
+ TH2F* BookNSigHistogram(const char * name);
+ TH2F* BookqVecHistogram(const char * name);
+ TH2F* BookGenMulvsRawMulHistogram(const char * name);
+
+ TH1F* GetPtHistogram1D(const char * name,Double_t minDCA,Double_t maxDCA);
+ TH1F* GetDCAHistogram1D(const char * name,Double_t minPt,Double_t maxPt);
+ TH2* GetHistogram(UInt_t id) { return (TH2*) fOutputList->At(id); }
+ // TH1* GetHistogram(BothHistoType_t histoType, BothParticleSpecies_t particleType, UInt_t charge);
+ TH1* GetHistogram1D(UInt_t histoType, UInt_t particleType, UInt_t charge);
+ TH2* GetHistogram2D(UInt_t histoType, UInt_t particleType, UInt_t charge);
+ TH2* GetPtHistogram(UInt_t id) { return (TH2*) fOutputList->At(id); }
+ TH2* GetPtHistogram(const char * name) { return (TH2*) fOutputList->FindObject(name); }
+ TH2* GetPtHistogramByName(UInt_t id) { return (TH2*) fOutputList->FindObject(kHistNameBoth[id]); } // Use this if you want to read a file saved with a different histo list
+ TH2* GetPIDHistogram(UInt_t id) { return (TH2*) fOutputList->At(id); }
+ TH2* GetPIDHistogram(const char * name) { return (TH2*) fOutputList->FindObject(name); }
+ TH2* GetPIDHistogramByName(UInt_t id) { return (TH2*) fOutputList->FindObject(kHistNameBoth[id]); }// Use this if you want to read a file saved with a different histo list
+ TH2* GetNSigHistogram(UInt_t id) { return (TH2*) fOutputList->At(id); }
+ TH2* GetNSigHistogram(const char * name) { return (TH2*) fOutputList->FindObject(name); }
+ TH2* GetNSigHistogramByName(UInt_t id) { return (TH2*) fOutputList->FindObject(kHistNameBoth[id]); }// Use this if you want to read a file saved with a different histo list
+ TH2* GetqVecHistogram(UInt_t id) { return (TH2*) fOutputList->At(id); }
+ TH2* GetqVecHistogram(const char * name) { return (TH2*) fOutputList->FindObject(name); }
+ TH2* GetqVecHistogramByName(UInt_t id) { return (TH2*) fOutputList->FindObject(kHistNameBoth[id]); }// Use this if you want to read a file saved with a different histo list
+ TH2* GetGenMulvsRawMulHistogram(const char * name) { return (TH2*) fOutputList->FindObject(name); }
+ TH2* GetGenMulvsRawMulHistogramByName(UInt_t id) { return (TH2*) fOutputList->FindObject(kHistNameBoth[id]); }// Use this if you want to read a file saved with a different histo list
+ //TH1F* GetTH1F(UInt_t id) { return (TH1F*) GetPtHistogram(id); }
+ //TH2F* GetTH2F(UInt_t id) { return (TH2F*) GetPIDHistogram(id); }
+
+ TList * GetOutputList() {return fOutputList;}
+ void SetNRebin(Int_t nreb){fNRebin=nreb;}
+ Int_t GetNRebin() {return fNRebin;}
+
+ Long64_t Merge(TCollection* list);
+
+
+private:
+ TList *fOutputList; // List of Pt Histo's
+ Int_t fNRebin; //rebin of histos
+ AliSpectraBothHistoManager(const AliSpectraBothHistoManager&);
+ AliSpectraBothHistoManager& operator=(const AliSpectraBothHistoManager&);
+
+ ClassDef(AliSpectraBothHistoManager, 1);
+
+};
+#endif
+
--- /dev/null
+#include "AliSpectraBothPID.h"
+#include "AliAODEvent.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TList.h"
+#include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
+#include "AliPIDResponse.h"
+#include "AliAnalysisManager.h"
+#include "AliInputEventHandler.h"
+#include "AliSpectraBothTrackCuts.h"
+
+ClassImp(AliSpectraBothPID)
+
+AliSpectraBothPID::AliSpectraBothPID() : TNamed("PID", "PID object"), fPIDType(kNSigmaTPCTOF), fNSigmaPID(3), fPIDResponse(0),fshiftTPC(0),fshiftTOF(0)
+{
+
+}
+
+AliSpectraBothPID::AliSpectraBothPID(BothPIDType_t pidType) : TNamed("PID", "PID object"), fPIDType(pidType), fNSigmaPID(3), fPIDResponse(0), fshiftTPC(0),fshiftTOF(0)
+{
+
+
+
+}
+
+
+
+void AliSpectraBothPID::FillQAHistos(AliSpectraBothHistoManager * hman, AliVTrack * track, AliSpectraBothTrackCuts * trackCuts)
+{
+
+ // fill a bunch of QA histos
+
+
+ // Get PID response object, if needed
+ if(!fPIDResponse)
+ {
+ AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
+ fPIDResponse = inputHandler->GetPIDResponse();
+ }
+
+ //Response
+ AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(track);
+
+ hman->GetPIDHistogram(kHistPIDTPC)->Fill(track->GetTPCmomentum(), track->GetTPCsignal()*track->Charge()); // PID histo
+ hman->GetPIDHistogram(kHistPIDTPCPion)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kPion)*track->Charge()); // Expected PIDPion histo
+ hman->GetPIDHistogram(kHistPIDTPCKaon)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kKaon)*track->Charge()); // Expected PIDKaon histo
+ hman->GetPIDHistogram(kHistPIDTPCProton)->Fill(track->GetTPCmomentum(),fPIDResponse->GetTPCResponse().GetExpectedSignal(track->GetTPCmomentum(),AliPID::kProton)*track->Charge()); // Expected PIDProton histo
+
+ Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
+ Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
+ Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
+ Double_t nsigmaTOFkProton=100.0,nsigmaTOFkKaon=100.0,nsigmaTOFkPion=100.0;
+
+ Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
+ Double_t nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
+ Double_t nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
+ if(fPIDType==kNSigmaTOF)
+ {
+ nsigmaTPCTOFkProton = 100.0;
+ nsigmaTPCTOFkKaon = 100.0;
+ nsigmaTPCTOFkPion =100.0;
+
+ }
+
+ if(track->Pt()>trackCuts->GetPtTOFMatching())
+ {
+
+ hman->GetPIDHistogram(kHistPIDTOF)->Fill(track->P(),(track->GetTOFsignal()/100)*track->Charge()); // PID histo
+
+ nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
+ nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF);
+ nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF);
+
+ //TOF
+ hman->GetPtHistogram(kHistNSigProtonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
+ hman->GetPtHistogram(kHistNSigKaonTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));
+ hman->GetPtHistogram(kHistNSigPionTOF)->Fill(track->P(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));
+ hman->GetPtHistogram(kHistNSigProtonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton));
+ hman->GetPtHistogram(kHistNSigKaonPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon));
+ hman->GetPtHistogram(kHistNSigPionPtTOF)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion));
+
+ nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2);
+ nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2);
+ nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2);
+
+ }
+ else
+ {
+ if (fPIDType==kNSigmaTOF)
+ return;
+ }
+
+ //TPC
+ hman->GetPtHistogram(kHistNSigProtonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
+ hman->GetPtHistogram(kHistNSigKaonTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
+ hman->GetPtHistogram(kHistNSigPionTPC)->Fill(track->P(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
+ hman->GetPtHistogram(kHistNSigProtonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
+ hman->GetPtHistogram(kHistNSigKaonPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
+ hman->GetPtHistogram(kHistNSigPionPtTPC)->Fill(track->Pt(),fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
+ //TPCTOF
+ hman->GetPtHistogram(kHistNSigProtonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkProton);
+ hman->GetPtHistogram(kHistNSigKaonTPCTOF)->Fill(track->P(),nsigmaTPCTOFkKaon);
+ hman->GetPtHistogram(kHistNSigPionTPCTOF)->Fill(track->P(),nsigmaTPCTOFkPion);
+ hman->GetPtHistogram(kHistNSigProtonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkProton);
+ hman->GetPtHistogram(kHistNSigKaonPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkKaon);
+ hman->GetPtHistogram(kHistNSigPionPtTPCTOF)->Fill(track->Pt(),nsigmaTPCTOFkPion);
+
+}
+
+Int_t AliSpectraBothPID::GetParticleSpecie(AliAODMCParticle * part)
+{
+ // return PID according to MC truth
+ switch(TMath::Abs(part->PdgCode()))
+ {
+ case 2212:
+ return kSpProton;
+ break;
+ case 321:
+ return kSpKaon;
+ break;
+ case 211:
+ return kSpPion;
+ break;
+ default:
+ return kSpUndefined;
+ }
+}
+
+Int_t AliSpectraBothPID::GetParticleSpecie(TParticle* part)
+{
+ // return PID according to MC truth
+ switch(TMath::Abs(part->GetPdgCode()))
+ {
+ case 2212:
+ return kSpProton;
+ break;
+ case 321:
+ return kSpKaon;
+ break;
+ case 211:
+ return kSpPion;
+ break;
+ default:
+ return kSpUndefined;
+ }
+}
+
+Int_t AliSpectraBothPID::GetParticleSpecie(AliSpectraBothHistoManager * hman,AliVTrack * trk, AliSpectraBothTrackCuts * trackCuts, Bool_t* rec)
+{
+ // return PID according to detectors
+ // Get PID response object, if needed
+
+ // guess the particle based on the smaller nsigma
+
+ rec[kSpPion]=false;
+ rec[kSpKaon]=false;
+ rec[kSpProton]=false;
+
+ if(!fPIDResponse)
+ {
+ AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
+ AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
+ fPIDResponse = inputHandler->GetPIDResponse();
+ }
+
+ if(!fPIDResponse)
+ {
+ AliFatal("Cannot get pid response");
+ return 0;
+ }
+
+
+ // Compute nsigma for each hypthesis
+ AliVParticle *inEvHMain = dynamic_cast<AliVParticle *>(trk);
+
+ // --- TPC
+ Double_t nsigmaTPCkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kProton)+fshiftTPC);
+ Double_t nsigmaTPCkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kKaon)+fshiftTPC);
+ Double_t nsigmaTPCkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(inEvHMain, AliPID::kPion)+fshiftTPC);
+
+
+
+ Double_t nsigmaTOFkProton=100.0,nsigmaTOFkKaon=100.0,nsigmaTOFkPion=100.0;
+
+ Double_t nsigmaTPCTOFkProton = TMath::Abs(nsigmaTPCkProton);
+ Double_t nsigmaTPCTOFkKaon = TMath::Abs(nsigmaTPCkKaon);
+ Double_t nsigmaTPCTOFkPion = TMath::Abs(nsigmaTPCkPion);
+ if(fPIDType==kNSigmaTOF)
+ {
+ nsigmaTPCTOFkProton = 100.0;
+ nsigmaTPCTOFkKaon = 100.0;
+ nsigmaTPCTOFkPion =100.0;
+
+ }
+
+
+
+ if(trk->Pt()>trackCuts->GetPtTOFMatching())
+ {
+ nsigmaTOFkProton = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kProton)+fshiftTOF);
+ nsigmaTOFkKaon = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kKaon)+fshiftTOF);
+ nsigmaTOFkPion = TMath::Abs(fPIDResponse->NumberOfSigmasTOF(inEvHMain, AliPID::kPion)+fshiftTOF);
+ // the TOF info is used in combined
+ nsigmaTPCTOFkProton = TMath::Sqrt((nsigmaTPCkProton*nsigmaTPCkProton+nsigmaTOFkProton*nsigmaTOFkProton)/2.);
+ nsigmaTPCTOFkKaon = TMath::Sqrt((nsigmaTPCkKaon*nsigmaTPCkKaon+nsigmaTOFkKaon*nsigmaTOFkKaon)/2.);
+ nsigmaTPCTOFkPion = TMath::Sqrt((nsigmaTPCkPion*nsigmaTPCkPion+nsigmaTOFkPion*nsigmaTOFkPion)/2.);
+ }
+ else
+ {
+ if (fPIDType==kNSigmaTOF)
+ return kSpUndefined;
+
+ }
+
+ // select the nsigma to be used for the actual PID
+ Double_t nsigmaPion=100;
+ Double_t nsigmaKaon=100;
+ Double_t nsigmaProton=100;
+
+ switch (fPIDType)
+ {
+ case kNSigmaTPC:
+ nsigmaProton = nsigmaTPCkProton;
+ nsigmaKaon = nsigmaTPCkKaon ;
+ nsigmaPion = nsigmaTPCkPion ;
+ break;
+ case kNSigmaTOF:
+ nsigmaProton = nsigmaTOFkProton;
+ nsigmaKaon = nsigmaTOFkKaon ;
+ nsigmaPion = nsigmaTOFkPion ;
+ break;
+ case kNSigmaTPCTOF:
+ nsigmaProton = nsigmaTPCTOFkProton;
+ nsigmaKaon = nsigmaTPCTOFkKaon ;
+ nsigmaPion = nsigmaTPCTOFkPion ;
+ break;
+ }
+
+ if(nsigmaPion < fNSigmaPID)
+ rec[kSpPion]=true;
+ if(nsigmaKaon < fNSigmaPID)
+ rec[kSpKaon]=true;
+ if(nsigmaProton < fNSigmaPID)
+ rec[kSpProton]=true;
+
+ if( ( nsigmaKaon==nsigmaPion ) && ( nsigmaKaon==nsigmaProton ))
+ return kSpUndefined;//if is the default value for the three
+ if( ( nsigmaKaon < nsigmaPion ) && ( nsigmaKaon < nsigmaProton ) && (nsigmaKaon < fNSigmaPID))
+ {
+ hman->GetPIDHistogram(kHistPIDTPCKaonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDKaon histo
+ return kSpKaon;
+ }
+ if( ( nsigmaPion < nsigmaKaon ) && ( nsigmaPion < nsigmaProton ) && (nsigmaPion < fNSigmaPID))
+ {
+ hman->GetPIDHistogram(kHistPIDTPCPionRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDPion histo
+ return kSpPion;
+ }
+ if( ( nsigmaProton < nsigmaKaon ) && ( nsigmaProton < nsigmaPion ) && (nsigmaProton < fNSigmaPID))
+ {
+ hman->GetPIDHistogram(kHistPIDTPCProtonRec)->Fill(trk->GetTPCmomentum(), trk->GetTPCsignal()*trk->Charge()); // Reconstructed PIDProton histo
+ return kSpProton;
+ }
+ // else, return undefined
+ return kSpUndefined;
+
+}
+
+Long64_t AliSpectraBothPID::Merge(TCollection* list)
+{
+ // Merging interface.
+ // Returns the number of merged objects (including this).
+
+ Printf("Merging");
+
+ if (!list)
+ return 0;
+
+ if (list->IsEmpty())
+ return 1;
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj;
+
+ // Actually, we don't do anything here...
+ // collections of all histograms
+ // TList collections;
+
+ Int_t count = 0;
+
+ while ((obj = iter->Next())) {
+ AliSpectraBothPID* entry = dynamic_cast<AliSpectraBothPID*> (obj);
+ if (entry == 0)
+ continue;
+
+ // TH1I * histo = entry->GetHistoCuts();
+ // collections.Add(histo);
+ count++;
+ }
+
+ // fHistoCuts->Merge(&collections);
+
+ delete iter;
+ Printf("OK");
+ return count+1;
+}
+
--- /dev/null
+
+#ifndef ALISPECTRABOTHPID_H
+#define ALISPECTRABOTHPID_H
+
+/* See cxx source for full Copyright notice */
+
+//-------------------------------------------------------------------------
+// AliSpectraBothPID
+//
+//
+//
+//
+// Authors: Michele Floris, CERN, Leonardo Milano, Torino
+//-------------------------------------------------------------------------
+
+class AliAODEvent;
+class TH1F;
+class TH2F;
+class TList;
+class AliVTrack;
+class AliAODMCParticle;
+class TParticle;
+class AliPIDResponse;
+class AliSpectraBothTrackCuts;
+
+#include "TNamed.h"
+#include "TParticle.h"
+#include "AliSpectraBothHistoManager.h"
+
+namespace AliSpectraNameSpaceBoth {
+
+ enum BothPIDType_t
+ {
+ kNSigmaTPC,
+ kNSigmaTOF,
+ kNSigmaTPCTOF, // squared sum
+ };
+
+
+
+}
+
+using namespace AliSpectraNameSpaceBoth;
+
+class AliSpectraBothPID : public TNamed
+{
+public:
+ AliSpectraBothPID() ;
+ AliSpectraBothPID(BothPIDType_t pidType);
+ virtual ~AliSpectraBothPID() {}
+
+ void FillQAHistos(AliSpectraBothHistoManager * hman, AliVTrack * track, AliSpectraBothTrackCuts * trackCuts) ;
+ void SetNSigmaCut(Float_t nsigma) { fNSigmaPID = nsigma; }
+ void SetPIDtype(BothPIDType_t pidType){fPIDType=pidType;}
+ void SetShiftTPC(Float_t shift){fshiftTPC=shift;}
+ void SetShiftTOF(Float_t shift){fshiftTOF=shift;}
+ Float_t GetNSigmaCut() {return fNSigmaPID; }
+
+ Int_t GetParticleSpecie(AliSpectraBothHistoManager * hman,AliVTrack * trk, AliSpectraBothTrackCuts * trackCuts, Bool_t* rec);
+ Int_t GetParticleSpecie(AliAODMCParticle * trk);
+ Int_t GetParticleSpecie(TParticle * trk);
+
+
+
+ Long64_t Merge(TCollection* list);
+
+
+private:
+
+ BothPIDType_t fPIDType; // PID type
+ Float_t fNSigmaPID; // number of sigma for PID cut
+ AliPIDResponse *fPIDResponse; // ! PID response object
+ Float_t fshiftTPC; // shift of the nsigma TPC
+ Float_t fshiftTOF; // shift of the nsigma TPC
+
+ AliSpectraBothPID(const AliSpectraBothPID&);
+ AliSpectraBothPID& operator=(const AliSpectraBothPID&);
+
+ ClassDef(AliSpectraBothPID, 1);
+
+};
+#endif
+
--- /dev/null
+
+/**************************************************************************
+ * Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+//-----------------------------------------------------------------
+// AliSpectraBothTrackCuts class
+//-----------------------------------------------------------------
+
+#include "TChain.h"
+#include "TTree.h"
+#include "TLegend.h"
+#include "TH1F.h"
+#include "TH1I.h"
+#include "TH2F.h"
+#include "TCanvas.h"
+#include "AliAnalysisTask.h"
+#include "AliAnalysisManager.h"
+#include "AliAODTrack.h"
+#include "AliVTrack.h"
+#include "AliExternalTrackParam.h"
+#include "AliAODMCParticle.h"
+#include "AliAODEvent.h"
+#include "AliAODInputHandler.h"
+#include "AliAnalysisTaskESDfilter.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliSpectraBothTrackCuts.h"
+//#include "AliSpectraBothHistoManager.h"
+#include <iostream>
+
+using namespace std;
+
+const char * AliSpectraBothTrackCuts::kBinLabel[] ={"TrkBit",
+ "TrkCuts",
+ "TrkEta",
+ "TrkDCA",
+ "TrkP",
+ "TrkPt",
+ "TrkPtTOF",
+ "TOFMatching",
+ "kTOFout",
+ "kTIME",
+ "kTOFpid",
+ "Accepted"};
+
+
+ClassImp(AliSpectraBothTrackCuts)
+
+
+AliSpectraBothTrackCuts::AliSpectraBothTrackCuts(const char *name) : TNamed(name, "AOD Track Cuts"), fIsSelected(0), fTrackBits(0), fMinTPCcls(0), fEtaCutMin(0), fEtaCutMax(0), fDCACut(0), fPCut(0), fPtCut(0), fYCut(0),
+ fPtCutTOFMatching(0),fAODtrack(kotherobject), fHashitinSPD1(0),
+fHistoCuts(0), fHistoNSelectedPos(0), fHistoNSelectedNeg(0), fHistoNMatchedPos(0), fHistoNMatchedNeg(0), fHistoEtaPhiHighPt(0), fHistoNclustersITS(0),fTrack(0),fCuts(0)
+
+{
+ // Constructor
+ fHistoCuts = new TH1I("fTrkCuts", "Track Cuts", kNTrkCuts, -0.5, kNTrkCuts - 0.5);
+ for(Int_t ibin=1;ibin<=kNTrkCuts;ibin++)fHistoCuts->GetXaxis()->SetBinLabel(ibin,kBinLabel[ibin-1]);
+ //standard histo
+ const Double_t templBins[] = {0.05,0.1,0.12,0.14,0.16,0.18,0.20,0.25,0.30,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.2,3.4,3.6,3.8,4.0,4.2,4.4,4.6,4.8,5.0};
+ Int_t nbinsTempl=52;
+
+ fHistoNSelectedPos=new TH1F("fHistoNSelectedPos","fHistoNSelectedPos",nbinsTempl,templBins);
+ fHistoNSelectedPos->GetXaxis()->SetTitle("P_{T} (GeV / c)");
+ fHistoNSelectedNeg=new TH1F("fHistoNSelectedNeg","fHistoNSelectedNeg",nbinsTempl,templBins);
+ fHistoNSelectedNeg->GetXaxis()->SetTitle("P_{T} (GeV / c)");
+ fHistoNMatchedPos=new TH1F("fHistoNMatchedPos","fHistoNMatchedPos",nbinsTempl,templBins);
+ fHistoNMatchedPos->GetXaxis()->SetTitle("P_{T} (GeV / c)");
+ fHistoNMatchedNeg=new TH1F("fHistoNMatchedNeg","fHistoNMatchedNeg",nbinsTempl,templBins);
+ fHistoNMatchedNeg->GetXaxis()->SetTitle("P_{T} (GeV / c)");
+ fHistoEtaPhiHighPt=new TH2F("fHistoEtaPhiHighPt","fHistoEtaPhiHighPt",200,-1,1,400,0,7);
+ fHistoEtaPhiHighPt->SetXTitle("eta");
+ fHistoEtaPhiHighPt->SetYTitle("phi");
+ fHistoNclustersITS=new TH1F("fHistoNclustersITS","fHistoNclustersITS;N;ITSLayer",6,-0.5,5.5);
+
+ fEtaCutMin = -100000.0; // default value of eta cut ~ no cut
+ fEtaCutMax = 100000.0; // default value of eta cut ~ no cut
+ fDCACut = 100000.0; // default value of dca cut ~ no cut
+ fPCut = 100000.0; // default value of p cut ~ no cut
+ fPtCut = 100000.0; // default value of pt cut ~ no cut
+ fPtCutTOFMatching=0.6; //default value fot matching with TOF
+ fYCut = 100000.0; // default value of y cut ~ no cut
+ fMinTPCcls=70; // ncls in TPC
+}
+
+//_______________________________________________________
+Bool_t AliSpectraBothTrackCuts::IsSelected(AliVTrack * track,Bool_t FillHistStat)
+{
+// Returns true if Track Cuts are selected and applied
+ if (!track)
+ {
+ printf("ERROR: Could not receive track");
+ return kFALSE;
+ }
+ fTrack = track;
+ TString nameoftrack(track->ClassName());
+ if(!nameoftrack.CompareTo("AliESDtrack"))
+ fAODtrack=kESDobject;
+ else if(!nameoftrack.CompareTo("AliAODTrack"))
+ fAODtrack=kAODobject;
+ else
+ fAODtrack=kotherobject;
+ if(!CheckTrackType()){
+ return kFALSE;
+ }
+ if(FillHistStat)fHistoCuts->Fill(kTrkBit);
+ if(!CheckTrackCuts()){
+ return kFALSE;
+ }
+ if(FillHistStat)fHistoCuts->Fill(kTrkCuts);
+ if(!CheckEtaCut()){
+ return kFALSE;
+ }
+ if(FillHistStat)fHistoCuts->Fill(kTrkEta);
+ if(!CheckDCACut()){
+ return kFALSE;
+ }
+ if(FillHistStat)fHistoCuts->Fill(kTrkDCA);
+ if(!CheckPCut()){
+ return kFALSE;
+ }
+ if(FillHistStat)fHistoCuts->Fill(kTrkP);
+ if(!CheckPtCut()){
+ return kFALSE;
+ }
+ if(FillHistStat)fHistoCuts->Fill(kTrkPt);
+ if(!CheckTOFMatching(FillHistStat)){
+ return kFALSE;
+ }
+ if(FillHistStat)fHistoCuts->Fill(kAccepted);
+ //Printf("-------- %d,%d",kTOFMatching,kAccepted);
+
+ return kTRUE;
+}
+//_________________________________________________________
+
+Bool_t AliSpectraBothTrackCuts::CheckTrackType()
+{
+ // Check track Type
+ if(fAODtrack==kESDobject)
+ {
+ AliESDtrack* esdtrack=dynamic_cast<AliESDtrack*>(fTrack);
+
+ if(fCuts->AcceptTrack(esdtrack)) return kTRUE;
+ return kFALSE;
+ }
+ else if(fAODtrack==kAODobject)
+ {
+ AliAODTrack* aodtrack=dynamic_cast<AliAODTrack*>(fTrack);
+ if (aodtrack->TestFilterBit(fTrackBits)) return kTRUE;
+ return kFALSE;
+ }
+
+ else
+ return kFALSE;
+
+}
+//_________________________________________________________
+
+Bool_t AliSpectraBothTrackCuts::CheckTrackCuts()
+{
+ // Check additional track Cuts
+ Bool_t PassTrackCuts=kTRUE;
+ AliAODTrack* aodtrack=0;
+ AliESDtrack* esdtrack=0;
+ if(fAODtrack==kESDobject)
+ {
+ esdtrack=dynamic_cast<AliESDtrack*>(fTrack);
+ if (!esdtrack->HasPointOnITSLayer(0) && !esdtrack->HasPointOnITSLayer(1))PassTrackCuts=kFALSE; //FIXME 1 SPD for the moment
+ if (fHashitinSPD1&&!esdtrack->HasPointOnITSLayer(0)) PassTrackCuts=kFALSE;
+ if (esdtrack->GetTPCNcls()<fMinTPCcls)PassTrackCuts=kFALSE;
+ if(!esdtrack->IsOn(AliESDtrack::kTPCrefit))PassTrackCuts=kFALSE;
+ if(!esdtrack->IsOn(AliESDtrack::kITSrefit))PassTrackCuts=kFALSE;
+ if(PassTrackCuts)
+ {
+ for(int i=0;i<6;i++)
+ if(esdtrack->HasPointOnITSLayer(i))
+ fHistoNclustersITS->Fill(i);
+ }
+ }
+ else if (fAODtrack==kAODobject)
+ {
+ aodtrack=dynamic_cast<AliAODTrack*>(fTrack);
+ if (!aodtrack->HasPointOnITSLayer(0) && !aodtrack->HasPointOnITSLayer(1))PassTrackCuts=kFALSE; //FIXME 1 SPD for the moment
+ if (fHashitinSPD1&&!aodtrack->HasPointOnITSLayer(0)) PassTrackCuts=kFALSE;
+ if (aodtrack->GetTPCNcls()<fMinTPCcls)PassTrackCuts=kFALSE;
+ if(!aodtrack->IsOn(AliAODTrack::kTPCrefit))PassTrackCuts=kFALSE;
+ if(!aodtrack->IsOn(AliAODTrack::kITSrefit))PassTrackCuts=kFALSE;
+ if(PassTrackCuts)
+ {
+ for(int i=0;i<6;i++)
+ if(aodtrack->HasPointOnITSLayer(i))
+ fHistoNclustersITS->Fill(i);
+ }
+ }
+ else
+ return kFALSE;
+
+
+
+ return PassTrackCuts;
+}
+//________________________________________________________
+Bool_t AliSpectraBothTrackCuts::CheckEtaCut()
+{
+ // Check eta cut
+ if (fTrack->Eta() < fEtaCutMax && fTrack->Eta() > fEtaCutMin) return kTRUE;
+ return kFALSE;
+}
+
+Bool_t AliSpectraBothTrackCuts::CheckYCut(BothParticleSpecies_t species)
+{
+ // check if the rapidity is within the set range
+ Double_t y;
+
+ Double_t pz=fTrack->Pz();
+ Double_t p=fTrack->P();
+ Double_t mass=-1.0;
+ /*
+ if (species == kSpProton) { y = fTrack->Y(9.38271999999999995e-01); }
+ if ( species == kSpKaon ) { y = fTrack->Y(4.93676999999999977e-01); }
+ if ( species == kSpPion) { y = fTrack->Y(1.39570000000000000e-01); }
+
+ */
+ if (species == kSpProton) { mass=9.38271999999999995e-01; }
+ if ( species == kSpKaon ) { mass=4.93676999999999977e-01; }
+ if ( species == kSpPion) { mass=1.39570000000000000e-01; }
+ if(mass<0.0)
+ y =-999.0 ;
+ else
+ y=0.5*TMath::Log((TMath::Sqrt(mass*mass+p*p)+pz)/(TMath::Sqrt(mass*mass+p*p)-pz));
+ if (TMath::Abs(y) > fYCut || y < -998.) return kFALSE;
+ return kTRUE;
+}
+//_______________________________________________________
+Bool_t AliSpectraBothTrackCuts::CheckDCACut()
+{
+ // Check DCA cut
+ // if (TMath::Abs(fTrack->DCA()) < fDCACut) return kTRUE; //FIXME for newest AOD fTrack->DCA() always gives -999
+
+ AliAODTrack* aodtrack=0;
+ AliESDtrack* esdtrack=0;
+ if(fAODtrack==kESDobject)
+ {
+ esdtrack=dynamic_cast<AliESDtrack*>(fTrack);
+ Float_t dcaxy=0.0;
+ Float_t dcaz=0.0;
+ esdtrack->GetImpactParameters(dcaxy,dcaz);
+ if (TMath::Abs(dcaxy) < fDCACut)
+ return kTRUE;
+ else
+ return kFALSE;
+ }
+ else if (fAODtrack==kAODobject)
+ {
+ aodtrack=dynamic_cast<AliAODTrack*>(fTrack);
+ if (TMath::Abs(aodtrack->DCA()) < fDCACut) return kTRUE;
+ else
+ return kFALSE;
+
+ }
+ else
+ return kFALSE;
+}
+//________________________________________________________
+Bool_t AliSpectraBothTrackCuts::CheckPCut()
+{
+ // Check P cut
+ if (fTrack->P() < fPCut) return kTRUE;
+ return kFALSE;
+}
+//_______________________________________________________
+Bool_t AliSpectraBothTrackCuts::CheckPtCut()
+{
+ // check Pt cut
+// if ((fTrack->Pt() < fPtCut) && (fTrack->Pt() > 0.3 )) return kTRUE;
+ if (fTrack->Pt() < fPtCut) return kTRUE;
+ return kFALSE;
+}
+
+//_______________________________________________________
+Bool_t AliSpectraBothTrackCuts::CheckTOFMatching(Bool_t FillHistStat)
+{
+ // check Pt cut
+ // if ((fTrack->Pt() < fPtCut) && (fTrack->Pt() > 0.3 )) return kTRUE;
+
+ if (fTrack->Pt() < fPtCutTOFMatching) return kTRUE;
+ else{
+ if(FillHistStat)fHistoCuts->Fill(kTrkPtTOF);
+ if(fTrack->Charge()>0)fHistoNSelectedPos->Fill(fTrack->Pt());
+ else fHistoNSelectedNeg->Fill(fTrack->Pt());
+ UInt_t status;
+ status=fTrack->GetStatus();
+ if((status&AliAODTrack::kTOFout)&&FillHistStat)fHistoCuts->Fill(kTrTOFout);
+ if((status&AliAODTrack::kTIME)&&FillHistStat)fHistoCuts->Fill(kTrTIME);
+ if((status&AliAODTrack::kTOFpid)&&FillHistStat)fHistoCuts->Fill(kTrTOFpid);
+
+ if((status&AliAODTrack::kTOFout)==0 || (status&AliAODTrack::kTIME)==0){//kTOFout and kTIME
+ return kFALSE;
+ }
+ if(FillHistStat)fHistoCuts->Fill(kTOFMatching);
+ if(fTrack->Charge()>0)fHistoNMatchedPos->Fill(fTrack->Pt());
+ else fHistoNMatchedNeg->Fill(fTrack->Pt());
+ if(fTrack->Pt()>1.5){
+ //fHistoEtaPhiHighPt->Fill(fTrack->GetOuterParam()->Eta(),fTrack->GetOuterParam()->Phi());
+ //Printf("AliExternalTrackParam * extpar=(AliExternalTrackParam*)fTrack->GetOuterParam();");
+ //AliExternalTrackParam * extpar=(AliExternalTrackParam*)fTrack->GetOuterParam();
+ fHistoEtaPhiHighPt->Fill(fTrack->Eta(),fTrack->Phi());
+ //Printf("fHistoEtaPhiHighPt->Fill(extpar->Eta(),extpar->Phi());");
+ //fHistoEtaPhiHighPt->Fill(extpar->Eta(),extpar->Phi());
+ //delete extpar;
+ }
+ return kTRUE;
+ }
+}
+//_______________________________________________________
+void AliSpectraBothTrackCuts::PrintCuts() const
+{
+ // Print cuts
+ cout << "Track Cuts" << endl;
+ cout << " > TrackBit\t" << fTrackBits << endl;
+ cout << " > Eta cut\t" << fEtaCutMin <<","<< fEtaCutMax << endl;
+ cout << " > DCA cut\t" << fDCACut << endl;
+ cout << " > P cut\t" << fPCut << endl;
+ cout << " > Pt cut \t" << fPtCut << endl;
+ cout << " > TPC cls \t" << fMinTPCcls << endl;
+}
+//_______________________________________________________
+void AliSpectraBothTrackCuts::SetTrackType(UInt_t bit)
+{
+ // Set the type of track to be used. The argument should be the bit number. The mask is produced automatically.
+ fTrackBits = (0x1 << (bit - 1));
+}
+//_______________________________________________________
+
+Long64_t AliSpectraBothTrackCuts::Merge(TCollection* list)
+{
+ // Merge a list of AliSpectraBothTrackCuts objects with this.
+ // Returns the number of merged objects (including this).
+
+ // AliInfo("Merging");
+
+
+ if (!list)
+ return 0;
+
+ if (list->IsEmpty())
+ return 1;
+
+ TIterator* iter = list->MakeIterator();
+ TObject* obj;
+
+ // collections of all histograms
+ TList collections;//FIXME we should only 1 collection
+ TList collections_histoNSelectedPos;
+ TList collections_histoNSelectedNeg;
+ TList collections_histoNMatchedPos;
+ TList collections_histoNMatchedNeg;
+ TList collections_histoEtaPhiHighPt;
+
+ Int_t count = 0;
+
+ while ((obj = iter->Next())) {
+ AliSpectraBothTrackCuts* entry = dynamic_cast<AliSpectraBothTrackCuts*> (obj);
+ if (entry == 0)
+ continue;
+
+ TH1I * histo = entry->GetHistoCuts();
+ collections.Add(histo);
+ TH1F * histoNSelectedPos = entry->GetHistoNSelectedPos();
+ collections_histoNSelectedPos.Add(histoNSelectedPos);
+ TH1F * histoNSelectedNeg = entry->GetHistoNSelectedNeg();
+ collections_histoNSelectedNeg.Add(histoNSelectedNeg);
+ TH1F * histoNMatchedPos = entry->GetHistoNMatchedPos();
+ collections_histoNMatchedPos.Add(histoNMatchedPos);
+ TH1F * histoNMatchedNeg = entry->GetHistoNMatchedNeg();
+ collections_histoNMatchedNeg.Add(histoNMatchedNeg);
+ TH2F * histoEtaPhiHighPt = entry->GetHistoEtaPhiHighPt();
+ collections_histoEtaPhiHighPt.Add(histoEtaPhiHighPt);
+ count++;
+ }
+
+ fHistoCuts->Merge(&collections);
+ fHistoNSelectedPos->Merge(&collections_histoNSelectedPos);
+ fHistoNSelectedNeg->Merge(&collections_histoNSelectedNeg);
+ fHistoNMatchedPos->Merge(&collections_histoNMatchedPos);
+ fHistoNMatchedNeg->Merge(&collections_histoNMatchedNeg);
+ fHistoEtaPhiHighPt->Merge(&collections_histoEtaPhiHighPt);
+
+ delete iter;
+
+ return count+1;
+}
+
--- /dev/null
+#ifndef ALISPECTRABOTHTRACKCUTS_H
+#define ALISPECTRABOTHTRACKCUTS_H
+
+/* See cxx source for full Copyright notice */
+
+//-------------------------------------------------------------------------
+// AliSpectraBothTrackCuts
+//
+//
+//
+//
+// Authors: Michele Floris, CERN, Philip Versteeg, UU, Redmer Bertens, UU
+//-------------------------------------------------------------------------
+
+class AliAODEvent;
+//class AliSpectraBothHistoManager;
+class TH1I;
+class AliAODMCParticle;
+class AliAODTrack;
+class AliESDtrackCuts;
+#include "AliSpectraBothHistoManager.h"
+#include "TNamed.h"
+#include "AliESDtrackCuts.h"
+
+using namespace AliSpectraNameSpaceBoth;
+
+class AliSpectraBothTrackCuts : public TNamed
+{
+ public:
+
+ enum { kTrkBit = 0, kTrkCuts, kTrkEta, kTrkDCA, kTrkP, kTrkPt,kTrkPtTOF,kTOFMatching,kTrTOFout,kTrTIME,kTrTOFpid,kAccepted,kNTrkCuts};
+ enum {kAODobject=0,kESDobject,kotherobject};
+
+ AliSpectraBothTrackCuts() : TNamed(), fIsSelected(0), fTrackBits(0), fMinTPCcls(0), fEtaCutMin(0), fEtaCutMax(0),fDCACut(0) ,fPCut(0), fPtCut(0),fYCut(0), fPtCutTOFMatching(0),fAODtrack(0), fHashitinSPD1(0),fHistoCuts(0), fHistoNSelectedPos(0), fHistoNSelectedNeg(0), fHistoNMatchedPos(0), fHistoNMatchedNeg(0), fHistoEtaPhiHighPt(0), fHistoNclustersITS(0),fTrack(0),fCuts(0) {}
+
+ AliSpectraBothTrackCuts(const char *name);
+ virtual ~AliSpectraBothTrackCuts() {} // To be implemented
+
+ Bool_t IsSelected(AliVTrack * track,Bool_t FillHistStat);
+
+ void SetTrackType(UInt_t bit);
+ Bool_t CheckTrackType();
+ Bool_t CheckTrackCuts();
+ Bool_t CheckEtaCut();
+ Bool_t CheckYCut(BothParticleSpecies_t specie); // not included in standard cuts
+ Bool_t CheckDCACut();
+ Bool_t CheckPCut();
+ Bool_t CheckPtCut();
+ Bool_t CheckTOFMatching(Bool_t FillHistStat);
+ void PrintCuts() const;
+
+ UInt_t GetTrackType() const { return fTrackBits;}
+ TH1I * GetHistoCuts() { return fHistoCuts; }
+ TH1F * GetHistoNSelectedPos() { return fHistoNSelectedPos; }
+ TH1F * GetHistoNSelectedNeg() { return fHistoNSelectedNeg; }
+ TH1F * GetHistoNMatchedPos() { return fHistoNMatchedPos; }
+ TH1F * GetHistoNMatchedNeg() { return fHistoNMatchedNeg; }
+ TH2F * GetHistoEtaPhiHighPt() { return fHistoEtaPhiHighPt; }
+ TH1F * GetHistoNclustersITS() {return fHistoNclustersITS;}
+ void SetEta(Float_t etamin,Float_t etamax) { fEtaCutMin = etamin;fEtaCutMax = etamax; }
+ void SetDCA(Float_t dca) { fDCACut = dca; }
+ void SetP(Float_t p) { fPCut = p; }
+ void SetPt(Float_t pt) { fPtCut = pt; }
+ void SetY(Float_t y) { fYCut = y;}
+ void SetPtTOFMatching(Float_t pt) { fPtCutTOFMatching = pt; }
+ void SetTrackBits(UInt_t TrackBits) {fTrackBits=TrackBits;}
+ void SetMinTPCcls(UInt_t MinTPCcls) {fMinTPCcls=MinTPCcls;}
+ void SetHashitinSPD1 (Bool_t value) {fHashitinSPD1=value;}
+ Float_t GetEtaMin() const { return fEtaCutMin; }
+ Float_t GetEtaMax() const { return fEtaCutMax; }
+ Float_t GetY() const { return fYCut; }
+ Float_t GetDCA() const { return fDCACut; }
+ Float_t GetP() const { return fPCut; }
+ Float_t GetPt() const { return fPtCut; }
+ Float_t GetPtTOFMatching() const { return fPtCutTOFMatching; }
+ Long64_t Merge(TCollection* list);
+ void SetAliESDtrackCuts(AliESDtrackCuts* cuts ){fCuts=cuts;}
+
+ private:
+
+ Bool_t fIsSelected; // True if cuts are selected
+ UInt_t fTrackBits; // Type of track to be used
+ UInt_t fMinTPCcls; // min number of clusters in the TPC
+ Float_t fEtaCutMin; // Allowed absolute maximum value of Eta
+ Float_t fEtaCutMax; // Allowed absolute maximum value of Eta
+ Float_t fDCACut; // Maximum value of DCA
+ Float_t fPCut; // Maximum value of P
+ Float_t fPtCut; // Maximum value of Pt
+ Float_t fYCut; // Maximum value of Y
+ Float_t fPtCutTOFMatching; // TOF Matching
+ Int_t fAODtrack; // 0 ESD track connected , 1 AOD track conected , else nothing
+ Bool_t fHashitinSPD1; // Check if SPD1 has a hit
+ TH1I *fHistoCuts; // Cuts statistics
+ TH1F *fHistoNSelectedPos; // Selected positive tracks
+ TH1F *fHistoNSelectedNeg; // Selected negative tracks
+ TH1F *fHistoNMatchedPos; // Matched positive tracks
+ TH1F *fHistoNMatchedNeg; // Matched negative tracks
+ TH2F *fHistoEtaPhiHighPt; // EtaPhi distr at high pt (>1.5 GeV/c)
+ TH1F *fHistoNclustersITS; // Number of clusters in ITS
+ AliVTrack *fTrack; //! Track pointer
+ AliESDtrackCuts *fCuts; //! cuts
+ static const char * kBinLabel[]; // labels of stat histo
+
+
+ AliSpectraBothTrackCuts(const AliSpectraBothTrackCuts&);
+ AliSpectraBothTrackCuts& operator=(const AliSpectraBothTrackCuts&);
+
+ ClassDef(AliSpectraBothTrackCuts, 3);
+};
+#endif
+
--- /dev/null
+//#include "Histograms.h"
+
+const char * AliSpectraNameSpaceBoth::kHistNameBoth[kNHist] =
+ {
+ "hHistPtGenTruePrimaryPionPlus",
+ "hHistPtGenTruePrimaryKaonPlus",
+ "hHistPtGenTruePrimaryProtonPlus",
+ "hHistPtGenTruePrimaryPionMinus",
+ "hHistPtGenTruePrimaryKaonMinus",
+ "hHistPtGenTruePrimaryProtonMinus",
+ "hHistPtGen",
+ "hHistPtRecSigmaPionPlus",
+ "hHistPtRecSigmaKaonPlus",
+ "hHistPtRecSigmaProtonPlus",
+ "hHistPtRecSigmaPionMinus",
+ "hHistPtRecSigmaKaonMinus",
+ "hHistPtRecSigmaProtonMinus",
+ "hHistPtRecTruePionPlus",
+ "hHistPtRecTrueKaonPlus",
+ "hHistPtRecTrueProtonPlus",
+ "hHistPtRecTruePionMinus",
+ "hHistPtRecTrueKaonMinus",
+ "hHistPtRecTrueProtonMinus",
+ "hHistPtRecTrueMuonPlus",
+ "hHistPtRecTrueMuonMinus",
+ "hHistPtRecPrimaryPionPlus",
+ "hHistPtRecPrimaryKaonPlus",
+ "hHistPtRecPrimaryProtonPlus",
+ "hHistPtRecPrimaryPionMinus",
+ "hHistPtRecPrimaryKaonMinus",
+ "hHistPtRecPrimaryProtonMinus",
+ "hHistPtRecPrimaryMuonPlus",
+ "hHistPtRecPrimaryMuonMinus",
+ "hHistPtRecSigmaPrimaryPionPlus",
+ "hHistPtRecSigmaPrimaryKaonPlus",
+ "hHistPtRecSigmaPrimaryProtonPlus",
+ "hHistPtRecSigmaPrimaryPionMinus",
+ "hHistPtRecSigmaPrimaryKaonMinus",
+ "hHistPtRecSigmaPrimaryProtonMinus",
+ "hHistPtRecSigmaSecondaryMaterialPionPlus",
+ "hHistPtRecSigmaSecondaryMaterialKaonPlus",
+ "hHistPtRecSigmaSecondaryMaterialProtonPlus",
+ "hHistPtRecSigmaSecondaryMaterialPionMinus",
+ "hHistPtRecSigmaSecondaryMaterialKaonMinus",
+ "hHistPtRecSigmaSecondaryMaterialProtonMinus",
+ "hHistPtRecSigmaSecondaryWeakDecayPionPlus",
+ "hHistPtRecSigmaSecondaryWeakDecayKaonPlus",
+ "hHistPtRecSigmaSecondaryWeakDecayProtonPlus",
+ "hHistPtRecSigmaSecondaryWeakDecayPionMinus",
+ "hHistPtRecSigmaSecondaryWeakDecayKaonMinus",
+ "hHistPtRecSigmaSecondaryWeakDecayProtonMinus",
+ "hHistPtRecTruePrimaryPionPlus",
+ "hHistPtRecTruePrimaryKaonPlus",
+ "hHistPtRecTruePrimaryProtonPlus",
+ "hHistPtRecTruePrimaryPionMinus",
+ "hHistPtRecTruePrimaryKaonMinus",
+ "hHistPtRecTruePrimaryProtonMinus",
+ "hHistPtRecTruePrimaryMuonPlus",
+ "hHistPtRecTruePrimaryMuonMinus",
+ "hHistPtRec",
+ "hHistPtRecPrimary",
+ "hHistPIDTPC",
+ "hHistPIDTOF",
+ "hHistPIDTPCPion",
+ "hHistPIDTPCKaon",
+ "hHistPIDTPCProton",
+ "hHistPIDTPCPionRec",
+ "hHistPIDTPCKaonRec",
+ "hHistPIDTPCProtonRec",
+ "hHistNSigPionTPC",
+ "hHistNSigKaonTPC",
+ "hHistNSigProtonTPC",
+ "hHistNSigPionPtTPC",
+ "hHistNSigKaonPtTPC",
+ "hHistNSigProtonPtTPC",
+ "hHistNSigPionTOF",
+ "hHistNSigKaonTOF",
+ "hHistNSigProtonTOF",
+ "hHistNSigPionPtTOF",
+ "hHistNSigKaonPtTOF",
+ "hHistNSigProtonPtTOF",
+ "hHistNSigPionTPCTOF",
+ "hHistNSigKaonTPCTOF",
+ "hHistNSigProtonTPCTOF",
+ "hHistNSigPionPtTPCTOF",
+ "hHistNSigKaonPtTPCTOF",
+ "hHistNSigProtonPtTPCTOF",
+ "hHistGenMulvsRawMul",
+ };
--- /dev/null
+#ifndef HISTOGRAMNAMESBOTH_H
+#define HISTOGRAMNAMESBOTH_H
+//#include "Histograms.h"
+//This file was generated automatically, please do not edit!!
+
+namespace AliSpectraNameSpaceBoth
+{
+ extern const char * kHistNameBoth[kNHist] ;
+}
+
+#endif
--- /dev/null
+
+// This file is used to give a list of histograms to be created by the manager.
+// the histogram names are automatically generated by the createNames.py script
+// the type/binning of the histograms depends on the range.
+// DON'T FORGET TO RUN createNames.py AFTER EDITING THIS FILE
+// IMPORTANT CONVENTIONS:
+// - don't assign numerical value explicitly to the entries (they would be skipped in the authomatic name generation)
+// - If you add an histogram set, please respect the order:
+// PionPlus, KaonPlus, ProtonPlus, PionMinus, KaonMinus, ProtonMinus (needed for getters)
+
+namespace AliSpectraNameSpaceBoth
+{
+ enum AODPtHist_t
+ {
+
+ // 6 Pt Generated True Primary
+ kHistPtGenTruePrimaryPionPlus, // Pt histo for pions +, generated tracks, true ID, primary Event
+ kHistPtGenTruePrimaryKaonPlus, // Pt histo for kaons +, generated tracks, true ID, primary Event
+ kHistPtGenTruePrimaryProtonPlus, // Pt histo for protons +, generated tracks, true ID, primary Event
+ kHistPtGenTruePrimaryPionMinus, // Pt histo for pions -, generated tracks, true ID, primary Event
+ kHistPtGenTruePrimaryKaonMinus, // Pt histo for kaons -, generated tracks, true ID, primary Event
+ kHistPtGenTruePrimaryProtonMinus, // Pt histo for protons -, generated tracks, true ID, primary Event
+ kNPtGenHist = kHistPtGenTruePrimaryProtonMinus, // Number of ptGen-likehistos histos - PID
+ kHistPtGen, // Pt histo for all particles, generated tracks
+ kNPtGenAllChHist = kHistPtGen, // Number of ptGen-likehistos histos - AllCh
+
+ // 6 Pt Reconstructed Sigma
+ kHistPtRecSigmaPionPlus, // Pt histo for pions +, reconstructed tracks, sigma ID
+ kHistPtRecSigmaKaonPlus, // Pt histo for kaons +, reconsructed tracks, sigma ID
+ kHistPtRecSigmaProtonPlus, // Pt histo for protons +, reconstructed tracks, sigma ID
+ kHistPtRecSigmaPionMinus, // Pt histo for pions -, reconstructed tracks, sigma ID
+ kHistPtRecSigmaKaonMinus, // Pt histo for kaons -, reconstructed tracks, sigma ID
+ kHistPtRecSigmaProtonMinus, // Pt histo for protons -, reconstructed tracks, sigma ID
+
+ // 6 Pt Reconstructed True & identified with nsigma
+ kHistPtRecTruePionPlus, // Pt histo for pions +, reconstructed tracks, true ID
+ kHistPtRecTrueKaonPlus, // Pt histo for kaons +, reconsructed tracks, true ID
+ kHistPtRecTrueProtonPlus, // Pt histo for protons +, reconstructed tracks, true ID
+ kHistPtRecTruePionMinus, // Pt histo for pions -, reconstructed tracks, true ID
+ kHistPtRecTrueKaonMinus, // Pt histo for kaons -, reconstructed tracks, true ID
+ kHistPtRecTrueProtonMinus, // Pt histo for protons -, reconstructed tracks, true ID
+ kHistPtRecTrueMuonPlus, // Pt histo for muons +, reconstructed tracks, true ID,
+ kHistPtRecTrueMuonMinus, // Pt histo for muons +, reconstructed tracks, true ID,
+
+ // 6 Pt Reconstructed True & (regardless of the offline nsigma identification)
+ kHistPtRecPrimaryPionPlus, // Pt histo for pions +, reconstructed tracks, true ID
+ kHistPtRecPrimaryKaonPlus, // Pt histo for kaons +, reconsructed tracks, true ID
+ kHistPtRecPrimaryProtonPlus, // Pt histo for protons +, reconstructed tracks, true ID
+ kHistPtRecPrimaryPionMinus, // Pt histo for pions -, reconstructed tracks, true ID
+ kHistPtRecPrimaryKaonMinus, // Pt histo for kaons -, reconstructed tracks, true ID
+ kHistPtRecPrimaryProtonMinus, // Pt histo for protons -, reconstructed tracks, true ID
+ kHistPtRecPrimaryMuonPlus, // Pt histo for muons +, reconstructed tracks, true ID,
+ kHistPtRecPrimaryMuonMinus, // Pt histo for muons +, reconstructed tracks, true ID,
+
+ // 6 Pt Reconstructed Sigma Primary
+ kHistPtRecSigmaPrimaryPionPlus, // Pt histo for pions +, reconstructed tracks, sigma ID, primary Event
+ kHistPtRecSigmaPrimaryKaonPlus, // Pt histo for kaons +, reconsructed tracks, sigma ID, primary Event
+ kHistPtRecSigmaPrimaryProtonPlus, // Pt histo for protons +, reconstructed tracks, sigma ID, primary Event
+ kHistPtRecSigmaPrimaryPionMinus, // Pt histo for pions -, reconstructed tracks, sigma ID, primary Event
+ kHistPtRecSigmaPrimaryKaonMinus, // Pt histo for kaons -, reconstructed tracks, sigma ID, primary Event
+ kHistPtRecSigmaPrimaryProtonMinus, // Pt histo for protons -, reconstructed tracks, sigma ID, primary Event
+
+ // 6 Pt Reconstructed Sigma Secondary Material
+ kHistPtRecSigmaSecondaryMaterialPionPlus, // Pt histo for pions +, reconstructed tracks, sigma ID, secondary Event
+ kHistPtRecSigmaSecondaryMaterialKaonPlus, // Pt histo for kaons +, reconsructed tracks, sigma ID, secondary Event
+ kHistPtRecSigmaSecondaryMaterialProtonPlus, // Pt histo for protons +, reconstructed tracks, sigma ID, secondary Event
+ kHistPtRecSigmaSecondaryMaterialPionMinus, // Pt histo for pions -, reconstructed tracks, sigma ID, secondary Event
+ kHistPtRecSigmaSecondaryMaterialKaonMinus, // Pt histo for kaons -, reconstructed tracks, sigma ID, secondary Event
+ kHistPtRecSigmaSecondaryMaterialProtonMinus, // Pt histo for protons -, reconstructed tracks, sigma ID, secondary Event
+
+ // 6 Pt Reconstructed Sigma Secondary WeakDecay
+ kHistPtRecSigmaSecondaryWeakDecayPionPlus, // Pt histo for pions +, reconstructed tracks, sigma ID, secondary Event
+ kHistPtRecSigmaSecondaryWeakDecayKaonPlus, // Pt histo for kaons +, reconsructed tracks, sigma ID, secondary Event
+ kHistPtRecSigmaSecondaryWeakDecayProtonPlus, // Pt histo for protons +, reconstructed tracks, sigma ID, secondary Event
+ kHistPtRecSigmaSecondaryWeakDecayPionMinus, // Pt histo for pions -, reconstructed tracks, sigma ID, secondary Event
+ kHistPtRecSigmaSecondaryWeakDecayKaonMinus, // Pt histo for kaons -, reconstructed tracks, sigma ID, secondary Event
+ kHistPtRecSigmaSecondaryWeakDecayProtonMinus, // Pt histo for protons -, reconstructed tracks, sigma ID, secondary Event
+
+ // 6 Pt Reconstructed True Primary
+ kHistPtRecTruePrimaryPionPlus, // Pt histo for pions +, reconstructed tracks, true ID, primary event
+ kHistPtRecTruePrimaryKaonPlus, // Pt histo for kaons +, reconsructed tracks, true ID, primary event
+ kHistPtRecTruePrimaryProtonPlus, // Pt histo for protons +, reconstructed tracks, true ID, primary event
+ kHistPtRecTruePrimaryPionMinus, // Pt histo for pions -, reconstructed tracks, true ID, primary event
+ kHistPtRecTruePrimaryKaonMinus, // Pt histo for kaons -, reconstructed tracks, true ID, primary event
+ kHistPtRecTruePrimaryProtonMinus, // Pt histo for protons -, reconstructed tracks, true ID, primary event
+ kHistPtRecTruePrimaryMuonPlus, // Pt histo for muons +, reconstructed tracks, true ID, primary event
+ kHistPtRecTruePrimaryMuonMinus, // Pt histo for muons +, reconstructed tracks, true ID, primary event
+ kNPtRecHist = kHistPtRecTruePrimaryMuonMinus, // Number of ptRec-likehistos histos
+
+ // Rest
+ kHistPtRec, // Pt histo for all particles, reconstructed tracks
+ kHistPtRecPrimary, // Pt histo for all particles, reconstructed tracks
+ kNPtRecAllChHist = kHistPtRecPrimary, // Number of ptRec-likehistos histos - no PID
+
+ kHistPIDTPC, // Particle Identification histo
+ kHistPIDTOF,
+ kHistPIDTPCPion,
+ kHistPIDTPCKaon,
+ kHistPIDTPCProton,
+ kHistPIDTPCPionRec,
+ kHistPIDTPCKaonRec,
+ kHistPIDTPCProtonRec,
+ kNHistPID =kHistPIDTPCProtonRec,
+
+ kHistNSigPionTPC,
+ kHistNSigKaonTPC,
+ kHistNSigProtonTPC, // NSigma separation plot
+ kHistNSigPionPtTPC,
+ kHistNSigKaonPtTPC,
+ kHistNSigProtonPtTPC,
+
+ kHistNSigPionTOF,
+ kHistNSigKaonTOF,
+ kHistNSigProtonTOF,
+ kHistNSigPionPtTOF,
+ kHistNSigKaonPtTOF,
+ kHistNSigProtonPtTOF,
+
+ kHistNSigPionTPCTOF,
+ kHistNSigKaonTPCTOF,
+ kHistNSigProtonTPCTOF,
+ kHistNSigPionPtTPCTOF,
+ kHistNSigKaonPtTPCTOF,
+ kHistNSigProtonPtTPCTOF,
+ kNHistNSig=kHistNSigProtonPtTPCTOF,
+
+ kHistGenMulvsRawMul,
+
+ kNHist, // Total number of histos
+ }; // Type of events plotted in Pt Histogram
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff