From 8df8e3822d6c956934c2d28eee71d8fa66cc4f7b Mon Sep 17 00:00:00 2001 From: andronic Date: Thu, 8 Jul 2010 16:53:05 +0000 Subject: [PATCH] Code updates with bugfixes; new PID class; filtered AOD config macro --- .../AliAnalysisTaskDielectronFilter.cxx | 72 +- .../AliAnalysisTaskDielectronFilter.h | 4 + PWG3/dielectron/AliDielectron.cxx | 12 +- PWG3/dielectron/AliDielectron.h | 10 +- PWG3/dielectron/AliDielectronCF.cxx | 20 +- PWG3/dielectron/AliDielectronCFdraw.cxx | 5 +- PWG3/dielectron/AliDielectronDebugTree.cxx | 20 +- PWG3/dielectron/AliDielectronDebugTree.h | 6 +- PWG3/dielectron/AliDielectronHistos.cxx | 8 + PWG3/dielectron/AliDielectronMC.cxx | 129 +++- PWG3/dielectron/AliDielectronMC.h | 33 +- PWG3/dielectron/AliDielectronPID.cxx | 377 ++++++++++ PWG3/dielectron/AliDielectronPID.h | 133 ++++ PWG3/dielectron/AliDielectronPair.cxx | 102 ++- PWG3/dielectron/AliDielectronPair.h | 12 +- PWG3/dielectron/AliDielectronPairLegCuts.cxx | 8 +- PWG3/dielectron/AliDielectronPairLegCuts.h | 2 +- PWG3/dielectron/AliDielectronSignalBase.cxx | 24 +- PWG3/dielectron/AliDielectronSignalBase.h | 10 +- PWG3/dielectron/AliDielectronSignalExt.h | 5 + PWG3/dielectron/AliDielectronSignalFunc.cxx | 28 +- PWG3/dielectron/AliDielectronSignalFunc.h | 2 + PWG3/dielectron/AliDielectronSpectrum.cxx | 40 + PWG3/dielectron/AliDielectronSpectrum.h | 16 +- PWG3/dielectron/AliDielectronVarManager.cxx | 18 + PWG3/dielectron/AliDielectronVarManager.h | 140 +++- PWG3/dielectron/{ => macros}/AddTaskJPSI.C | 3 - PWG3/dielectron/macros/AddTaskJPSIFilter.C | 54 ++ PWG3/dielectron/macros/ConfigJpsi2eeData.C | 58 +- PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C | 2 - PWG3/dielectron/macros/ConfigJpsi2eeEff.C | 254 +++++++ PWG3/dielectron/macros/ConfigJpsi2eeFilter.C | 156 ++++ PWG3/dielectron/macros/MakeDataReport.C | 708 ++++++++++++++++++ PWG3/dielectron/macros/MakeEfficiencyReport.C | 135 ++++ 34 files changed, 2482 insertions(+), 124 deletions(-) create mode 100644 PWG3/dielectron/AliDielectronPID.cxx create mode 100644 PWG3/dielectron/AliDielectronPID.h rename PWG3/dielectron/{ => macros}/AddTaskJPSI.C (95%) create mode 100644 PWG3/dielectron/macros/AddTaskJPSIFilter.C create mode 100644 PWG3/dielectron/macros/ConfigJpsi2eeEff.C create mode 100644 PWG3/dielectron/macros/ConfigJpsi2eeFilter.C create mode 100644 PWG3/dielectron/macros/MakeDataReport.C create mode 100644 PWG3/dielectron/macros/MakeEfficiencyReport.C diff --git a/PWG3/dielectron/AliAnalysisTaskDielectronFilter.cxx b/PWG3/dielectron/AliAnalysisTaskDielectronFilter.cxx index 89826cf79ec..e5d4dab8500 100644 --- a/PWG3/dielectron/AliAnalysisTaskDielectronFilter.cxx +++ b/PWG3/dielectron/AliAnalysisTaskDielectronFilter.cxx @@ -25,17 +25,22 @@ #include #include #include +#include +#include #include "AliDielectron.h" +#include "AliDielectronMC.h" #include "AliDielectronHistos.h" +#include "AliDielectronVarManager.h" #include "AliAnalysisTaskDielectronFilter.h" ClassImp(AliAnalysisTaskDielectronFilter) //_________________________________________________________________________________ AliAnalysisTaskDielectronFilter::AliAnalysisTaskDielectronFilter() : - AliAnalysisTaskSE(), - fDielectron(0) + AliAnalysisTaskSE(), + fDielectron(0), + fSelectPhysics(kTRUE) { // // Constructor @@ -44,8 +49,9 @@ AliAnalysisTaskDielectronFilter::AliAnalysisTaskDielectronFilter() : //_________________________________________________________________________________ AliAnalysisTaskDielectronFilter::AliAnalysisTaskDielectronFilter(const char *name) : - AliAnalysisTaskSE(name), - fDielectron(0) + AliAnalysisTaskSE(name), + fDielectron(0), + fSelectPhysics(kTRUE) { // // Constructor @@ -69,6 +75,7 @@ void AliAnalysisTaskDielectronFilter::Init() Error("Init","Dielectron framework class required. Please create and instance with proper cuts and set it via 'SetDielectron' before executing this task!!!"); return; } + fDielectron->Init(); aodH->AddFilteredAOD("AliAOD.Dielectron.root", "DielectronEvents"); // AddAODBranch("AliDielectronCandidates",fDielectron->GetPairArraysPointer(),"deltaAOD.Dielectron.root"); @@ -80,18 +87,69 @@ void AliAnalysisTaskDielectronFilter::UserExec(Option_t *) // // Main loop. Called for every event // - + if (!fDielectron) return; + AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager(); + AliESDInputHandler *esdHandler=0x0; + if ( (esdHandler=dynamic_cast(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){ + AliDielectronVarManager::SetESDpid(esdHandler->GetESDpid()); + } else { + //load esd pid bethe bloch parameters depending on the existance of the MC handler + // yes: MC parameters + // no: data parameters + if (!AliDielectronVarManager::GetESDpid()){ + if (AliDielectronMC::Instance()->HasMC()) { + AliDielectronVarManager::InitESDpid(); + } else { + AliDielectronVarManager::InitESDpid(1); + } + } + } + // Was event selected ? + AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler()); + Bool_t isSelected = kTRUE; + if( fSelectPhysics && inputHandler && inputHandler->GetEventSelection() ) { + isSelected = inputHandler->IsEventSelected(); + } + + if (!isSelected) return; + //bz for AliKF Double_t bz = InputEvent()->GetMagneticField(); AliKFParticle::SetField( bz ); fDielectron->Process(InputEvent()); - + if(fDielectron->HasCandidates()){ + //If input event is an AliESDevent + // replace the references of the legs with the AOD references + if(man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class()){ + AliAODEvent *aod = ((AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler()))->GetAOD(); + TObjArray *obj = 0x0; + AliAODTrack *leg1 = 0x0; + AliAODTrack *leg2 = 0x0; + for(Int_t i=0; i < 10; i++ ){ + obj = (TObjArray*)(*(fDielectron->GetPairArraysPointer()))->UncheckedAt(i); + if(!obj) continue; + for(int j=0;jGetEntriesFast();j++) + { + AliDielectronPair *pairObj = (AliDielectronPair*)obj->UncheckedAt(j); + Int_t id1 = ((AliESDtrack*)pairObj->GetFirstDaughter())->GetID(); + Int_t id2 = ((AliESDtrack*)pairObj->GetSecondDaughter())->GetID(); + for(Int_t it=0;itGetNumberOfTracks();it++){ + if(aod->GetTrack(it)->GetID() == id1) leg1 = aod->GetTrack(it); + if(aod->GetTrack(it)->GetID() == id2) leg2 = aod->GetTrack(it); + } + if(!leg1 || !leg2) continue; + pairObj->SetRefFirstDaughter(leg1); + pairObj->SetRefSecondDaughter(leg2); + } + } + } + AliAODExtension *extDielectron = dynamic_cast - ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler())->GetFilteredAOD("AliAOD.Dielectron.root"); + ((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler())->GetFilteredAOD("AliAOD.Dielectron.root"); extDielectron->SelectEvent(); //see if dielectron candidate branch exists, if not create is TTree *t=extDielectron->GetTree(); diff --git a/PWG3/dielectron/AliAnalysisTaskDielectronFilter.h b/PWG3/dielectron/AliAnalysisTaskDielectronFilter.h index 167d5957b7f..ff41d04e067 100644 --- a/PWG3/dielectron/AliAnalysisTaskDielectronFilter.h +++ b/PWG3/dielectron/AliAnalysisTaskDielectronFilter.h @@ -37,6 +37,8 @@ public: virtual void UserExec(Option_t *option); virtual void Init(); virtual void LocalInit() {Init();} + + void UsePhysicsSelection(Bool_t phy=kTRUE) {fSelectPhysics=phy;} void SetDielectron(AliDielectron * const die) { fDielectron = die; } @@ -44,6 +46,8 @@ private: AliDielectron *fDielectron; // J/psi framework object + Bool_t fSelectPhysics; // Whether to use physics selection + AliAnalysisTaskDielectronFilter(const AliAnalysisTaskDielectronFilter &c); AliAnalysisTaskDielectronFilter& operator= (const AliAnalysisTaskDielectronFilter &c); diff --git a/PWG3/dielectron/AliDielectron.cxx b/PWG3/dielectron/AliDielectron.cxx index a76cb2e92e7..56ee8c268aa 100644 --- a/PWG3/dielectron/AliDielectron.cxx +++ b/PWG3/dielectron/AliDielectron.cxx @@ -90,6 +90,8 @@ AliDielectron::AliDielectron() : fTrackFilter("TrackFilter"), fPairFilter("PairFilter"), fPdgMother(443), + fPdgLeg1(11), + fPdgLeg2(11), fHistos(0x0), fPairCandidates(new TObjArray(10)), fCfManagerPair(0x0), @@ -108,6 +110,8 @@ AliDielectron::AliDielectron(const char* name, const char* title) : fTrackFilter("TrackFilter"), fPairFilter("PairFilter"), fPdgMother(443), + fPdgLeg1(11), + fPdgLeg2(11), fHistos(0x0), fPairCandidates(new TObjArray(10)), fCfManagerPair(0x0), @@ -137,7 +141,7 @@ void AliDielectron::Init() // Initialise objects // if (fCfManagerPair) fCfManagerPair->InitialiseContainer(fPairFilter); - + if (fDebugTree) fDebugTree->SetDielectron(this); } //________________________________________________________________ @@ -295,9 +299,9 @@ void AliDielectron::FillPairArrays(Int_t arr1, Int_t arr2) { Int_t end=ntrack2; if (arr1==arr2) end=itrack1; for (Int_t itrack2=0; itrack2SetTracks(static_cast(fTracks[arr1].UncheckedAt(itrack1)), 11, - static_cast(fTracks[arr2].UncheckedAt(itrack2)), 11); + //create the pair + candidate->SetTracks(static_cast(fTracks[arr1].UncheckedAt(itrack1)), fPdgLeg1, + static_cast(fTracks[arr2].UncheckedAt(itrack2)), fPdgLeg2); candidate->SetType(pairIndex); candidate->SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother)); diff --git a/PWG3/dielectron/AliDielectron.h b/PWG3/dielectron/AliDielectron.h index 6f6f8267283..0f202bb7fdf 100644 --- a/PWG3/dielectron/AliDielectron.h +++ b/PWG3/dielectron/AliDielectron.h @@ -54,7 +54,11 @@ public: AliAnalysisFilter& GetTrackFilter() { return fTrackFilter; } AliAnalysisFilter& GetPairFilter() { return fPairFilter; } - void SetMotherPdg( Int_t pdgMother ) { fPdgMother=pdgMother; } + void SetMotherPdg( Int_t pdgMother ) { fPdgMother=pdgMother; } + void SetLegPdg(Int_t pdgLeg1, Int_t pdgLeg2) { fPdgLeg1=pdgLeg1; fPdgLeg2=pdgLeg2; } + Int_t GetMotherPdg() const { return fPdgMother; } + Int_t GetLeg1Pdg() const { return fPdgLeg1; } + Int_t GetLeg2Pdg() const { return fPdgLeg2; } const TObjArray* GetTrackArray(Int_t i) const {return (i>=0&&i<4)?&fTracks[i]:0;} const TObjArray* GetPairArray(Int_t i) const {return (i>=0&&i<10)? @@ -85,6 +89,8 @@ private: AliAnalysisFilter fPairFilter; // pair cuts Int_t fPdgMother; // pdg code of mother tracks + Int_t fPdgLeg1; // pdg code leg1 + Int_t fPdgLeg2; // pdg code leg2 AliDielectronHistos *fHistos; // Histogram manager // Streaming and merging should be handled @@ -124,7 +130,7 @@ private: AliDielectron(const AliDielectron &c); AliDielectron &operator=(const AliDielectron &c); - ClassDef(AliDielectron,2); + ClassDef(AliDielectron,3); }; inline void AliDielectron::InitPairCandidateArrays() diff --git a/PWG3/dielectron/AliDielectronCF.cxx b/PWG3/dielectron/AliDielectronCF.cxx index 02e689ae929..3d6c68eb17a 100644 --- a/PWG3/dielectron/AliDielectronCF.cxx +++ b/PWG3/dielectron/AliDielectronCF.cxx @@ -240,7 +240,7 @@ void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter) //Pure MC truth if (fStepForMCtruth){ - fCfContainer->SetStepTitle(step++,"MC truth"); + fCfContainer->SetStepTitle(step++,"MC truth"); } //before cuts (MC truth) @@ -317,8 +317,8 @@ void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter) cutName+="&"; cutName+=filter.GetCuts()->At(iCut)->GetName(); } - fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut } + fCfContainer->SetStepTitle(step++, cutName.Data()); //Step for the cut if (fHasMC){ if (fStepsForSignal) fCfContainer->SetStepTitle(step++, (cutName+" (Signal)").Data()); //Step for the cut with MC truth @@ -328,7 +328,7 @@ void AliDielectronCF::InitialiseContainer(const AliAnalysisFilter& filter) } if (step!=fNSteps) { - AliError("Something went wrong in the naming of the steps!!!"); + AliError(Form("Something went wrong in the naming of the steps!!! (%d != %d)",step,fNSteps)); } } @@ -371,7 +371,7 @@ void AliDielectronCF::Fill(UInt_t mask, const AliDielectronPair *particle) // step 0 would be full MC truth and is handled in FillMC Int_t step=0; if (fStepForMCtruth) ++step; - + //No cuts (MC truth) if (fStepForNoCutsMCmotherPid){ if (isMCTruth) fCfContainer->Fill(fValues,step); @@ -483,6 +483,15 @@ void AliDielectronCF::FillMC(const TObject *particle) Double_t valuesPair[AliDielectronVarManager::kNMaxValues]; AliDielectronVarManager::Fill(particle,valuesPair); + + AliVParticle *d1=0x0; + AliVParticle *d2=0x0; + AliDielectronMC::Instance()->GetDaughters(particle,d1,d2); + + valuesPair[AliDielectronVarManager::kThetaHE]=AliDielectronPair::ThetaPhiCM(d1,d2,kTRUE,kTRUE); + valuesPair[AliDielectronVarManager::kPhiHE]=AliDielectronPair::ThetaPhiCM(d1,d2,kTRUE,kFALSE); + valuesPair[AliDielectronVarManager::kThetaCS]=AliDielectronPair::ThetaPhiCM(d1,d2,kFALSE,kTRUE); + valuesPair[AliDielectronVarManager::kPhiCS]=AliDielectronPair::ThetaPhiCM(d1,d2,kFALSE,kFALSE); //TODO: temporary solution, set manually the pair type to 1: unlikesign SE valuesPair[AliDielectronVarManager::kPairType]=1; @@ -493,9 +502,6 @@ void AliDielectronCF::FillMC(const TObject *particle) } if (fNVarsLeg>0){ - AliVParticle *d1=0x0; - AliVParticle *d2=0x0; - AliDielectronMC::Instance()->GetDaughters(particle,d1,d2); Double_t valuesLeg1[AliDielectronVarManager::kNMaxValues]; Double_t valuesLeg2[AliDielectronVarManager::kNMaxValues]; if (d1->Pt()>d2->Pt()){ diff --git a/PWG3/dielectron/AliDielectronCFdraw.cxx b/PWG3/dielectron/AliDielectronCFdraw.cxx index 5656f9cb5b0..ca1088f1cae 100644 --- a/PWG3/dielectron/AliDielectronCFdraw.cxx +++ b/PWG3/dielectron/AliDielectronCFdraw.cxx @@ -605,10 +605,9 @@ void AliDielectronCFdraw::Draw(const TObjArray *arr, const char* opt) // } TCanvas *c=gPad->GetCanvas(); - if (!gVirtualPS&&!drawSamePlus) c->Clear(); + if (!gVirtualPS&&!drawSamePlus&&!drawSame&&nPads>1) c->Clear(); - - if (!drawSame){ + if (!drawSame&&nPads>1){ //optimised division Int_t nCols = (Int_t)TMath::Ceil( TMath::Sqrt(nPads) ); Int_t nRows = (Int_t)TMath::Ceil( (Double_t)nPads/(Double_t)nCols ); diff --git a/PWG3/dielectron/AliDielectronDebugTree.cxx b/PWG3/dielectron/AliDielectronDebugTree.cxx index 373d6473b9c..82b9d09849b 100644 --- a/PWG3/dielectron/AliDielectronDebugTree.cxx +++ b/PWG3/dielectron/AliDielectronDebugTree.cxx @@ -38,6 +38,8 @@ NOTE: Please use with extream care! Only for debugging and test purposes!!! #include #include +#include "AliDielectron.h" +#include "AliDielectronMC.h" #include "AliDielectronPair.h" #include "AliDielectronDebugTree.h" @@ -49,7 +51,8 @@ AliDielectronDebugTree::AliDielectronDebugTree() : fFileName("jpsi_debug.root"), fNVars(0), fNVarsLeg(0), - fStreamer(0x0) + fStreamer(0x0), + fDielectron(0x0) { // // Default Constructor @@ -66,7 +69,8 @@ AliDielectronDebugTree::AliDielectronDebugTree(const char* name, const char* tit fFileName("jpsi_debug.root"), fNVars(0), fNVarsLeg(0), - fStreamer(0x0) + fStreamer(0x0), + fDielectron(0x0) { // // Named Constructor @@ -106,6 +110,7 @@ void AliDielectronDebugTree::Fill(AliDielectronPair *pair) //Get File and event information TObjString fileName; Int_t eventInFile=-1; + Int_t runNumber=-1; TTree *t=man->GetTree(); if (t) { @@ -117,6 +122,7 @@ void AliDielectronDebugTree::Fill(AliDielectronPair *pair) if (han){ AliESDEvent *ev=dynamic_cast(han->GetEvent()); eventInFile=ev->GetEventNumberInFile(); + runNumber=ev->GetRunNumber(); } if (!fStreamer) fStreamer=new TTreeSRedirector(fFileName.Data()); @@ -126,9 +132,19 @@ void AliDielectronDebugTree::Fill(AliDielectronPair *pair) (*fStreamer) << "Pair" << "File.=" << &fileName << "EventInFile=" << eventInFile + << "Run=" << runNumber << "Leg1_ID=" << id1 << "Leg2_ID=" << id2; + //Fill MC information + Bool_t hasMC=AliDielectronMC::Instance()->HasMC(); + if (hasMC){ + Int_t pdg=443; + if (fDielectron) pdg=fDielectron->GetMotherPdg(); + Bool_t isMotherMC = AliDielectronMC::Instance()->IsMotherPdg(pair,pdg); + (*fStreamer) << "Pair" + << "mcTruth=" << isMotherMC; + } Int_t var=0; Double_t values[AliDielectronVarManager::kNMaxValues]; diff --git a/PWG3/dielectron/AliDielectronDebugTree.h b/PWG3/dielectron/AliDielectronDebugTree.h index bf46792f695..a4c003ffe95 100644 --- a/PWG3/dielectron/AliDielectronDebugTree.h +++ b/PWG3/dielectron/AliDielectronDebugTree.h @@ -26,6 +26,7 @@ class TTreeSRedirector; class AliDielectronPair; +class AliDielectron; class AliDielectronDebugTree : public TNamed { public: @@ -41,6 +42,8 @@ public: void Fill(AliDielectronPair *pair); + void SetDielectron(AliDielectron * const dielectron) { fDielectron=dielectron; } + void DeleteStreamer(); void WriteTree(); private: @@ -52,7 +55,8 @@ private: Int_t fVariablesLeg[AliDielectronVarManager::kNMaxValues]; //configured variables for the legs TTreeSRedirector *fStreamer; //! Tree Redirector - + AliDielectron *fDielectron; //! pointer to mother dielectron manager + AliDielectronDebugTree(const AliDielectronDebugTree &c); AliDielectronDebugTree &operator=(const AliDielectronDebugTree &c); diff --git a/PWG3/dielectron/AliDielectronHistos.cxx b/PWG3/dielectron/AliDielectronHistos.cxx index 160d673ac5f..d5e3cd8bc2c 100644 --- a/PWG3/dielectron/AliDielectronHistos.cxx +++ b/PWG3/dielectron/AliDielectronHistos.cxx @@ -445,6 +445,14 @@ void AliDielectronHistos::Draw(const Option_t* option) if ( TMath::Abs(h->GetXaxis()->GetBinWidth(1)-h->GetXaxis()->GetBinWidth(2))>1e-10 ) gPad->SetLogx(); if ( TMath::Abs(h->GetYaxis()->GetBinWidth(1)-h->GetYaxis()->GetBinWidth(2))>1e-10 ) gPad->SetLogy(); if ( TMath::Abs(h->GetZaxis()->GetBinWidth(1)-h->GetZaxis()->GetBinWidth(2))>1e-10 ) gPad->SetLogz(); + TString histOpt=h->GetOption(); + histOpt.ToLower(); + if (histOpt.Contains("logx")) gPad->SetLogx(); + if (histOpt.Contains("logy")) gPad->SetLogy(); + if (histOpt.Contains("logz")) gPad->SetLogz(); + histOpt.ReplaceAll("logx",""); + histOpt.ReplaceAll("logy",""); + histOpt.ReplaceAll("logz",""); h->Draw(drawOpt.Data()); } if (gVirtualPS) c->Update(); diff --git a/PWG3/dielectron/AliDielectronMC.cxx b/PWG3/dielectron/AliDielectronMC.cxx index e2b55fefaf5..a612dae8aff 100644 --- a/PWG3/dielectron/AliDielectronMC.cxx +++ b/PWG3/dielectron/AliDielectronMC.cxx @@ -51,7 +51,7 @@ AliDielectronMC* AliDielectronMC::Instance() else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODHandler::Class()) type=kAOD; AliMCEventHandler* mcHandler = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); - + fgInstance=new AliDielectronMC(type); fgInstance->SetHasMC(mcHandler!=0x0); @@ -115,6 +115,7 @@ AliVParticle* AliDielectronMC::GetMCTrackFromMCEvent(Int_t _itrk) // // return MC track directly from MC event // + if (_itrk<0) return NULL; if (!fMCEvent){ AliError("No fMCEvent"); return NULL;} AliVParticle * track = fMCEvent->GetTrack(_itrk); // tracks from MC event return track; @@ -126,8 +127,12 @@ Bool_t AliDielectronMC::ConnectMCEvent() // // connect stack object from the mc handler // + fMCEvent=0x0; AliMCEventHandler* mcHandler = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); if (!mcHandler){ AliError("Could not retrive MC event handler!"); return kFALSE; } + if (!mcHandler->InitOk() ) return kFALSE; + if (!mcHandler->TreeK() ) return kFALSE; + if (!mcHandler->TreeTR() ) return kFALSE; AliMCEvent* mcEvent = mcHandler->MCEvent(); if (!mcEvent){ AliError("Could not retrieve MC event!"); return kFALSE; } @@ -151,7 +156,7 @@ Bool_t AliDielectronMC::UpdateStack() } //____________________________________________________________ -AliMCParticle* AliDielectronMC::GetMCTrack(AliESDtrack* _track) +AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track) { // // return MC track @@ -162,7 +167,7 @@ AliMCParticle* AliDielectronMC::GetMCTrack(AliESDtrack* _track) } //____________________________________________________________ -TParticle* AliDielectronMC::GetMCTrackFromStack(AliESDtrack* _track) +TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track) { // // return MC track from stack @@ -175,21 +180,38 @@ TParticle* AliDielectronMC::GetMCTrackFromStack(AliESDtrack* _track) } //____________________________________________________________ -AliMCParticle* AliDielectronMC::GetMCTrackMother(AliESDtrack* _track) +AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliESDtrack* _track) { // // return MC track mother // AliMCParticle* mcpart = GetMCTrack(_track); if (!mcpart) return NULL; -printf("mcpart->GetMother() : %d\n",mcpart->GetMother()); AliMCParticle* mcmother = dynamic_cast(fMCEvent->GetTrack(mcpart->GetMother())); if (!mcmother) return NULL; return mcmother; } //____________________________________________________________ -TParticle* AliDielectronMC::GetMCTrackMotherFromStack(AliESDtrack* _track) +AliMCParticle* AliDielectronMC::GetMCTrackMother(const AliMCParticle* _particle){ + // + // return MC track mother + // + AliMCParticle* mcmother = dynamic_cast(fMCEvent->GetTrack(_particle->GetMother())); + return mcmother; +} + +//____________________________________________________________ +AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){ + // + // return MC track mother + // + AliAODMCParticle* mcmother = dynamic_cast(fMCEvent->GetTrack(_particle->GetMother())); + return mcmother; +} + +//____________________________________________________________ +TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track) { // // return MC track mother from stack @@ -202,7 +224,7 @@ TParticle* AliDielectronMC::GetMCTrackMotherFromStack(AliESDtrack* _track) } //____________________________________________________________ -Int_t AliDielectronMC::GetMCPID(AliESDtrack* _track) +Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track) { // // return PDG code of the track from the MC truth info @@ -213,7 +235,7 @@ Int_t AliDielectronMC::GetMCPID(AliESDtrack* _track) } //____________________________________________________________ -Int_t AliDielectronMC::GetMCPIDFromStack(AliESDtrack* _track) +Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track) { // // return MC PDG code from stack @@ -224,7 +246,7 @@ Int_t AliDielectronMC::GetMCPIDFromStack(AliESDtrack* _track) } //____________________________________________________________ -Int_t AliDielectronMC::GetMotherPDG(AliESDtrack* _track) +Int_t AliDielectronMC::GetMotherPDG( const AliESDtrack* _track) { // // return PDG code of the mother track from the MC truth info @@ -235,7 +257,7 @@ Int_t AliDielectronMC::GetMotherPDG(AliESDtrack* _track) } //____________________________________________________________ -Int_t AliDielectronMC::GetMotherPDGFromStack(AliESDtrack* _track) +Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track) { // // return PDG code of the mother track from stack @@ -246,7 +268,7 @@ Int_t AliDielectronMC::GetMotherPDGFromStack(AliESDtrack* _track) } //____________________________________________________________ -Int_t AliDielectronMC::GetMCProcess(AliESDtrack* _track) +Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track) { // // return process number of the track @@ -257,7 +279,7 @@ Int_t AliDielectronMC::GetMCProcess(AliESDtrack* _track) } //____________________________________________________________ -Int_t AliDielectronMC::GetMCProcessFromStack(AliESDtrack* _track) +Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track) { // // return process number of the track @@ -268,7 +290,42 @@ Int_t AliDielectronMC::GetMCProcessFromStack(AliESDtrack* _track) } //____________________________________________________________ -Int_t AliDielectronMC::GetMCProcessMother(AliESDtrack* _track) +Int_t AliDielectronMC::NumberOfDaughters(const AliESDtrack* track) +{ + // + // returns the number of daughters + // + AliMCParticle *mcmother=GetMCTrackMother(track); + if(!mcmother||!mcmother->Particle()) return -999; +// return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0; + return mcmother->Particle()->GetNDaughters(); +} + +//____________________________________________________________ +Int_t AliDielectronMC::NumberOfDaughters(const AliMCParticle* particle) +{ + // + // returns the number of daughters + // + AliMCParticle *mcmother=GetMCTrackMother(particle); + if(!mcmother||!mcmother->Particle()) return -999; +// return mcmother->GetFirstDaughter()>0?mcmother->GetLastDaughter()-mcmother->GetFirstDaughter()+1:0; + return mcmother->Particle()->GetNDaughters(); +} + +//____________________________________________________________ +Int_t AliDielectronMC::NumberOfDaughters(const AliAODMCParticle* particle) +{ + // + // returns the number of daughters + // + AliAODMCParticle *mcmother=GetMCTrackMother(particle); + if(!mcmother) return -999; + return mcmother->GetNDaughters(); +} + +//____________________________________________________________ +Int_t AliDielectronMC::GetMCProcessMother(const AliESDtrack* _track) { // // return process number of the mother of the track @@ -279,7 +336,7 @@ Int_t AliDielectronMC::GetMCProcessMother(AliESDtrack* _track) } //____________________________________________________________ -Int_t AliDielectronMC::GetMCProcessMotherFromStack(AliESDtrack* _track) +Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track) { // // return process number of the mother of the track @@ -297,6 +354,7 @@ Bool_t AliDielectronMC::IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMo // if (!fMCEvent) return kFALSE; + if (!particle) return kFALSE; if (particle->IsA()==AliMCParticle::Class()){ return IsMCMotherToEEesd(static_cast(particle),pdgMother); @@ -326,7 +384,8 @@ Bool_t AliDielectronMC::IsMCMotherToEEesd(const AliMCParticle *particle, Int_t p if ((ilast-ifirst)!=1) return kFALSE; AliMCParticle *firstD=static_cast(GetMCTrackFromMCEvent(ifirst)); AliMCParticle *secondD=static_cast(GetMCTrackFromMCEvent(ilast)); - + + //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!! if (firstD->Charge()>0){ if (firstD->PdgCode()!=-11) return kFALSE; if (secondD->PdgCode()!=11) return kFALSE; @@ -356,6 +415,7 @@ Bool_t AliDielectronMC::IsMCMotherToEEaod(const AliAODMCParticle *particle, Int_ AliAODMCParticle *firstD=static_cast(GetMCTrackFromMCEvent(ifirst)); AliAODMCParticle *secondD=static_cast(GetMCTrackFromMCEvent(ilast)); + //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!! if (firstD->Charge()>0){ if (firstD->GetPdgCode()!=11) return kFALSE; if (secondD->GetPdgCode()!=-11) return kFALSE; @@ -387,6 +447,7 @@ Int_t AliDielectronMC::GetLabelMotherWithPdgESD(const AliVParticle *particle1, c // test if mother of particle 1 and 2 has pdgCode +-11 (electron), // have the same mother and the mother had pdg code pdgMother // ESD case + //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!! // AliMCParticle *mcPart1=static_cast(GetMCTrackFromMCEvent(particle1->GetLabel())); @@ -414,6 +475,7 @@ Int_t AliDielectronMC::GetLabelMotherWithPdgAOD(const AliVParticle *particle1, c // test if mother of particle 1 and 2 has pdgCode +-11 (electron), // have the same mother and the mother had pdg code pdgMother // AOD case + //TODO: check how you can get rid of the hardcoded numbers. One should make use of the PdgCodes set in AliDielectron!!! // AliAODMCParticle *mcPart1=static_cast(GetMCTrackFromMCEvent(particle1->GetLabel())); AliAODMCParticle *mcPart2=static_cast(GetMCTrackFromMCEvent(particle2->GetLabel())); @@ -458,3 +520,40 @@ void AliDielectronMC::GetDaughters(const TObject *mother, AliVParticle* &d1, Ali d1=fMCEvent->GetTrack(lblD1); d2=fMCEvent->GetTrack(lblD2); } + +//____________________________________________________________ +Bool_t AliDielectronMC::HaveSameMother(const AliDielectronPair * pair) +{ + // + // Check whether two particles have the same mother + // + + const AliVParticle * daughter1 = pair->GetFirstDaughter(); + const AliVParticle * daughter2 = pair->GetSecondDaughter(); + + AliMCParticle *mcDaughter1=static_cast(GetMCTrackFromMCEvent(daughter1->GetLabel())); + AliMCParticle *mcDaughter2=static_cast(GetMCTrackFromMCEvent(daughter2->GetLabel())); + if (!mcDaughter1 || !mcDaughter2) return 0; + + TParticle *tDaughter1 = mcDaughter1->Particle(); + TParticle *tDaughter2 = mcDaughter2->Particle(); + if(!tDaughter1 || !tDaughter2)return 0; + + TParticle *mother1 = fStack->Particle(tDaughter1->GetMother(0)); + TParticle *mother2 = fStack->Particle(tDaughter2->GetMother(0)); + if(!mother1 || !mother2) return 0; + + + if(mother1->IsEqual(mother2)) return 1 ; + + return 0; + +} + + + + + + + + diff --git a/PWG3/dielectron/AliDielectronMC.h b/PWG3/dielectron/AliDielectronMC.h index 0079dcaa9ad..ab774156dc0 100644 --- a/PWG3/dielectron/AliDielectronMC.h +++ b/PWG3/dielectron/AliDielectronMC.h @@ -42,14 +42,14 @@ public: void Initialize(); // initialization Int_t GetNMCTracks(); // return number of generated tracks Int_t GetNMCTracksFromStack(); // return number of generated tracks from stack - Int_t GetMCPID(AliESDtrack* _track); // return MC PID - Int_t GetMCPIDFromStack(AliESDtrack* _track); // return MC PID - Int_t GetMotherPDG(AliESDtrack* _track); // return mother PID from the MC stack - Int_t GetMotherPDGFromStack(AliESDtrack* _track); // return mother PID from the MC stack - Int_t GetMCProcess(AliESDtrack* _track); // return process number - Int_t GetMCProcessFromStack(AliESDtrack* _track); // return process number - Int_t GetMCProcessMother(AliESDtrack* _track); // return process number of the mother track - Int_t GetMCProcessMotherFromStack(AliESDtrack* _track); // return process number of the mother track + Int_t GetMCPID(const AliESDtrack* _track); // return MC PID + Int_t GetMCPIDFromStack(const AliESDtrack* _track); // return MC PID + Int_t GetMotherPDG(const AliESDtrack* _track); // return mother PID from the MC stack + Int_t GetMotherPDGFromStack(const AliESDtrack* _track); // return mother PID from the MC stack + Int_t GetMCProcess(const AliESDtrack* _track); // return process number + Int_t GetMCProcessFromStack(const AliESDtrack* _track); // return process number + Int_t GetMCProcessMother(const AliESDtrack* _track); // return process number of the mother track + Int_t GetMCProcessMotherFromStack(const AliESDtrack* _track); // return process number of the mother track Bool_t ConnectMCEvent(); Bool_t UpdateStack(); @@ -57,16 +57,25 @@ public: Bool_t IsMotherPdg(const AliDielectronPair* pair, Int_t pdgMother); Bool_t IsMotherPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother); Bool_t IsMCMotherToEE(const AliVParticle *particle, Int_t pdgMother); + + Bool_t HaveSameMother(const AliDielectronPair *pair); Int_t GetLabelMotherWithPdg(const AliDielectronPair* pair, Int_t pdgMother); Int_t GetLabelMotherWithPdg(const AliVParticle *particle1, const AliVParticle *particle2, Int_t pdgMother); AliVParticle* GetMCTrackFromMCEvent(AliVParticle *track); // return MC track directly from MC event AliVParticle* GetMCTrackFromMCEvent(Int_t _itrk); // return MC track directly from MC event - TParticle* GetMCTrackFromStack(AliESDtrack* _track); // return MC track from stack - AliMCParticle* GetMCTrack(AliESDtrack* _track); // return MC track associated with reco track - TParticle* GetMCTrackMotherFromStack(AliESDtrack* _track); // return MC mother track from stack - AliMCParticle* GetMCTrackMother(AliESDtrack* _track); // return MC mother track from stack + TParticle* GetMCTrackFromStack(const AliESDtrack* _track); // return MC track from stack + AliMCParticle* GetMCTrack(const AliESDtrack* _track); // return MC track associated with reco track + + TParticle* GetMCTrackMotherFromStack(const AliESDtrack* _track); // return MC mother track from stack + AliMCParticle* GetMCTrackMother(const AliESDtrack* _track); // return MC mother track from stack + AliMCParticle* GetMCTrackMother(const AliMCParticle* _particle); // return MC mother track from stack + AliAODMCParticle* GetMCTrackMother(const AliAODMCParticle* _particle); // return MC mother track from stack + + Int_t NumberOfDaughters(const AliESDtrack* track); // return number of daughters + Int_t NumberOfDaughters(const AliMCParticle* particle); // return number of daughters + Int_t NumberOfDaughters(const AliAODMCParticle* particle); // return number of daughters void GetDaughters(const TObject *mother, AliVParticle* &d1, AliVParticle* &d2); diff --git a/PWG3/dielectron/AliDielectronPID.cxx b/PWG3/dielectron/AliDielectronPID.cxx new file mode 100644 index 00000000000..308cbdcf962 --- /dev/null +++ b/PWG3/dielectron/AliDielectronPID.cxx @@ -0,0 +1,377 @@ +/************************************************************************* +* 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. * +**************************************************************************/ + +/////////////////////////////////////////////////////////////////////////// +// Dielectron PID // +// // +// // +/* +Detailed description + + +*/ +// // +/////////////////////////////////////////////////////////////////////////// + +#include +#include + +#include +#include +#include + +#include "AliDielectronVarManager.h" + +#include "AliDielectronPID.h" + +ClassImp(AliDielectronPID) + +AliDielectronPID::AliDielectronPID() : + AliAnalysisCuts(), + fNcuts(0), + fRequirePIDbit(kTRUE), + fESDpid(0x0) +{ + // + // Default Constructor + // + for (Int_t icut=0; icut(track); + //TODO: Which momentum to use? + // Different momenta for different detectors? + Double_t mom=part->P(); + + Bool_t selected=kFALSE; + fESDpid=AliDielectronVarManager::GetESDpid(); + + for (UChar_t icut=0; icut1e-20) && (mom<=pMin || mom>pMax) ) continue; + + // test if we are supposed to use a function for the cut + if (fFunUpperCut[icut]) fNsigmaUp[icut] =fFunUpperCut[icut]->Eval(mom); + if (fFunLowerCut[icut]) fNsigmaLow[icut]=fFunLowerCut[icut]->Eval(mom); + + switch (fDetType[icut]){ + case kITS: + selected = IsSelectedITS(part,icut); + break; + case kTPC: + selected = IsSelectedTPC(part,icut); + break; + case kTRD: + selected = IsSelectedTRD(part,icut); + break; + case kTOF: + selected = IsSelectedTOF(part,icut); + break; + } + if (!selected) return kFALSE; + } + + return selected; +} + +//______________________________________________ +Bool_t AliDielectronPID::IsSelectedITS(AliVParticle * const part, Int_t icut) const +{ + // + // ITS part of the PID check + // Don't accept the track if there was no pid bit set + // + Float_t numberOfSigmas=-1000.; + + if (part->IsA()==AliESDtrack::Class()){ + // ESD case in case the PID bit is not set, don't use this track! + AliESDtrack *track=static_cast(part); + if (fRequirePIDbit&&!(track->GetStatus()&AliESDtrack::kITSpid)) return kFALSE; + + numberOfSigmas=fESDpid->NumberOfSigmasITS(track, fPartType[icut]); + }else{ + // AOD case + // FIXME: Is there a place to check whether the PID is was set in ESD??? + AliAODTrack *track=static_cast(part); + numberOfSigmas=NumberOfSigmasITS(track, fPartType[icut]); + } + Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut]; + return selected; +} + +//______________________________________________ +Bool_t AliDielectronPID::IsSelectedTPC(AliVParticle * const part, Int_t icut) const +{ + // + // TPC part of the PID check + // Don't accept the track if there was no pid bit set + // + Float_t numberOfSigmas=-1000.; + + if (part->IsA()==AliESDtrack::Class()){ + // ESD case in case the PID bit is not set, don't use this track! + AliESDtrack *track=static_cast(part); + if (fRequirePIDbit&&!(track->GetStatus()&AliESDtrack::kTPCpid)) return kFALSE; + + numberOfSigmas=fESDpid->NumberOfSigmasTPC(track, fPartType[icut]); + }else{ + // AOD case + // FIXME: Is there a place to check whether the PID is was set in ESD??? + AliAODTrack *track=static_cast(part); + numberOfSigmas=NumberOfSigmasTPC(track, fPartType[icut]); + } + + Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut]; + return selected; +} + +//______________________________________________ +Bool_t AliDielectronPID::IsSelectedTRD(AliVParticle * const /*part*/, Int_t /*icut*/) const +{ + // + // TRD part of the pid check + // + return kFALSE; +} + +//______________________________________________ +Bool_t AliDielectronPID::IsSelectedTOF(AliVParticle * const part, Int_t icut) const +{ + // + // TOF part of the PID check + // Don't accept the track if there was no pid bit set + // + Float_t numberOfSigmas=-1000.; + + if (part->IsA()==AliESDtrack::Class()){ + // ESD case in case the PID bit is not set, don't use this track! + AliESDtrack *track=static_cast(part); + if (fRequirePIDbit&&!(track->GetStatus()&AliESDtrack::kTOFpid)) return kFALSE; + + numberOfSigmas=fESDpid->NumberOfSigmasTOF(track, fPartType[icut], fESDpid->GetTOFResponse().GetTimeZero()); + }else{ + // AOD case + // FIXME: Is there a place to check whether the PID is was set in ESD??? + AliAODTrack *track=static_cast(part); + numberOfSigmas=NumberOfSigmasTOF(track, fPartType[icut]); + } + + Bool_t selected=((numberOfSigmas>=fNsigmaLow[icut])&&(numberOfSigmas<=fNsigmaUp[icut]))^fExclude[icut]; + return selected; +} + +//______________________________________________ +void AliDielectronPID::SetDefaults(Int_t def){ + // + // initialise default pid strategies + // + + if (def==0){ + // 2sigma bands TPC: + // - include e + // - exclude mu,K,pi,p + // -complete p range + AddCut(kTPC,AliPID::kElectron,2); + AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,0.,kTRUE); + AddCut(kTPC,AliPID::kPion,-2.,2.,0.,0.,kTRUE); + AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,0.,kTRUE); + AddCut(kTPC,AliPID::kProton,-2.,2.,0.,0.,kTRUE); + } else if (def==1) { + // 2sigma bands TPC: + // - include e 0SetParameters(-2.7,-0.4357); + AddCut(kTPC,AliPID::kElectron,lowerCut,3.); + AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5); + } else if (def==7) { + // lower cut TPC: parametrisation by HFE + // upper cut TPC: 3 sigma + // TOF ele band 3sigma 0 +#include +#include +#include + +#include + +class TF1; + +class TList; + +class AliDielectronPID : public AliAnalysisCuts { +public: + enum DetType {kITS, kTPC, kTRD, kTOF}; + + AliDielectronPID(); + AliDielectronPID(const char*name, const char* title); + + virtual ~AliDielectronPID(); + + void AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, Double_t nSigmaUp=-99999., + Double_t pMin=0, Double_t pMax=0, Bool_t exclude=kFALSE); + + void AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, TF1 * const funUp, + Double_t pMin=0, Double_t pMax=0, Bool_t exclude=kFALSE); + + void AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, Double_t nSigmaUp, + Double_t pMin=0, Double_t pMax=0, Bool_t exclude=kFALSE); + + void AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, TF1 * const funUp, + Double_t pMin=0, Double_t pMax=0, Bool_t exclude=kFALSE); + + void SetDefaults(Int_t def); + + void SetRequirePIDbit(Bool_t require) {fRequirePIDbit=require;} + // + //Analysis cuts interface + //const + virtual Bool_t IsSelected(TObject* track); + virtual Bool_t IsSelected(TList* /* list */ ) {return kFALSE;} + +private: + enum {kNmaxPID=10}; + + DetType fDetType[kNmaxPID]; //detector type of nsigma cut + AliPID::EParticleType fPartType[kNmaxPID]; //particle type + Float_t fNsigmaLow[kNmaxPID]; //lower nsigma bound + Float_t fNsigmaUp[kNmaxPID]; //upper nsigma bound + Double_t fPmin[kNmaxPID]; //lower momentum + Double_t fPmax[kNmaxPID]; //upper momentum + Bool_t fExclude[kNmaxPID]; //use as exclusion band + TF1 *fFunUpperCut[kNmaxPID];//use function as upper cut + TF1 *fFunLowerCut[kNmaxPID];//use function as lower cut + UChar_t fNcuts; //number of cuts + + Bool_t fRequirePIDbit; //If to require the PID bits to be set. + + AliESDpid *fESDpid; //! esd pid object + + + Bool_t IsSelectedITS(AliVParticle * const part, Int_t icut) const; + Bool_t IsSelectedTPC(AliVParticle * const part, Int_t icut) const; + Bool_t IsSelectedTRD(AliVParticle * const part, Int_t icut) const; + Bool_t IsSelectedTOF(AliVParticle * const part, Int_t icut) const; + + Float_t NumberOfSigmasITS(const AliAODTrack *track, AliPID::EParticleType type) const; + Float_t NumberOfSigmasTPC(const AliAODTrack *track, AliPID::EParticleType type) const; + Float_t NumberOfSigmasTOF(const AliAODTrack *track, AliPID::EParticleType type) const; + + AliDielectronPID(const AliDielectronPID &c); + AliDielectronPID &operator=(const AliDielectronPID &c); + + ClassDef(AliDielectronPID,1) // Dielectron PID +}; + + +// +// Inline functions for AOD as long as ther is no AOD pid object we have to fake it +// + +inline Float_t AliDielectronPID::NumberOfSigmasITS(const AliAODTrack *track, AliPID::EParticleType type) const { + AliAODPid *pid=track->GetDetPid(); + if (!pid) return -1000.; + + return fESDpid->GetITSResponse().GetNumberOfSigmas(track->P(),pid->GetITSsignal(),type); +} + +inline Float_t AliDielectronPID::NumberOfSigmasTPC(const AliAODTrack *track, AliPID::EParticleType type) const { + AliAODPid *pid=track->GetDetPid(); + if (!pid) return -1000.; + + Double_t mom = pid->GetTPCmomentum(); + if (mom<0) mom=track->P(); + + //FIXME: rough estimate of the number of clusters used for PID. Needs to be fixed!!! + Int_t ncl=(Int_t)track->GetTPCClusterMap().CountBits(); + return fESDpid->GetTPCResponse().GetNumberOfSigmas(mom,pid->GetTPCsignal(),ncl,type); +} + +inline Float_t AliDielectronPID::NumberOfSigmasTOF(const AliAODTrack *track, AliPID::EParticleType type) const { + AliAODPid *pid=track->GetDetPid(); + if (!pid) return -1000.; + + Double_t times[AliPID::kSPECIES]; + pid->GetIntegratedTimes(times); + Double_t tofRes = fESDpid->GetTOFResponse().GetExpectedSigma(track->P(),times[type],AliPID::ParticleMass(type)); + return (pid->GetTOFsignal() - times[type])/ tofRes; +} + + + +#endif diff --git a/PWG3/dielectron/AliDielectronPair.cxx b/PWG3/dielectron/AliDielectronPair.cxx index ef30d2795b2..ad6994c8e41 100644 --- a/PWG3/dielectron/AliDielectronPair.cxx +++ b/PWG3/dielectron/AliDielectronPair.cxx @@ -22,6 +22,7 @@ #include "AliDielectronPair.h" #include "AliVTrack.h" +#include "AliPID.h" ClassImp(AliDielectronPair) @@ -77,13 +78,13 @@ void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1, fPair.Initialize(); fD1.Initialize(); fD2.Initialize(); - + AliKFParticle kf1(*particle1,pid1); AliKFParticle kf2(*particle2,pid2); fPair.AddDaughter(kf1); fPair.AddDaughter(kf2); - + if (particle1->Pt()>particle2->Pt()){ fRefD1 = particle1; fRefD2 = particle2; @@ -97,3 +98,100 @@ void AliDielectronPair::SetTracks(AliVTrack * const particle1, Int_t pid1, } } +//______________________________________________ +Double_t AliDielectronPair::ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2, + const Bool_t isHE, const Bool_t isTheta) { + // The function calculates theta and phi in the mother rest frame with + // respect to the helicity coordinate system and Collins-Soper coordinate system + // TO DO: generalize for different decays (only J/Psi->e+e- now) + + // Laboratory frame 4-vectors: + // projectile beam & target beam 4-mom + const Double_t kBeamEnergy = 3500.; //TODO: need to retrieve the beam energy from somewhere + TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))); + TLorentzVector targMom(0.,0.,kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))); + + // first & second daughter 4-mom + TLorentzVector p1Mom(d1->Px(),d1->Py(),d1->Pz(),TMath::Sqrt(d1->Px()*d1->Px()+d1->Py()*d1->Py()+d1->Pz()*d1->Pz()+AliPID::ParticleMass(AliPID::kElectron)*AliPID::ParticleMass(AliPID::kElectron))); + TLorentzVector p2Mom(d2->Px(),d2->Py(),d2->Pz(),TMath::Sqrt(d2->Px()*d2->Px()+d2->Py()*d2->Py()+d2->Pz()*d2->Pz()+AliPID::ParticleMass(AliPID::kElectron)*AliPID::ParticleMass(AliPID::kElectron))); + // J/Psi 4-momentum vector + TLorentzVector motherMom=p1Mom+p2Mom; + + // boost all the 4-mom vectors to the mother rest frame + TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect(); + p1Mom.Boost(beta); + p2Mom.Boost(beta); + projMom.Boost(beta); + targMom.Boost(beta); + + // x,y,z axes + TVector3 zAxis; + if(isHE) zAxis = (motherMom.Vect()).Unit(); + else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit(); + TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit(); + TVector3 xAxis = (yAxis.Cross(zAxis)).Unit(); + + // return either theta or phi + if(isTheta) { + if(d1->Charge()>0) + return zAxis.Dot((p1Mom.Vect()).Unit()); + else + return zAxis.Dot((p2Mom.Vect()).Unit()); + + } + else { + if(d1->Charge()>0) + return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis)); + else + return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis)); + } +} + +//______________________________________________ +Double_t AliDielectronPair::ThetaPhiCM(const Bool_t isHE, const Bool_t isTheta) const { + // The function calculates theta and phi in the mother rest frame with + // respect to the helicity coordinate system and Collins-Soper coordinate system + // TO DO: generalize for different decays (only J/Psi->e+e- now) + + // Laboratory frame 4-vectors: + // projectile beam & target beam 4-mom + const Double_t kBeamEnergy = 3500.; //TODO: need to retrieve the beam energy from somewhere + TLorentzVector projMom(0.,0.,-kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))); + TLorentzVector targMom(0.,0.,kBeamEnergy,TMath::Sqrt(kBeamEnergy*kBeamEnergy+AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))); + + // first & second daughter 4-mom + AliVParticle *d1 = dynamic_cast(fRefD1.GetObject()); + AliVParticle *d2 = dynamic_cast(fRefD2.GetObject()); + TLorentzVector p1Mom(d1->Px(),d1->Py(),d1->Pz(),TMath::Sqrt(d1->Px()*d1->Px()+d1->Py()*d1->Py()+d1->Pz()*d1->Pz()+AliPID::ParticleMass(AliPID::kElectron)*AliPID::ParticleMass(AliPID::kElectron))); + TLorentzVector p2Mom(d2->Px(),d2->Py(),d2->Pz(),TMath::Sqrt(d2->Px()*d2->Px()+d2->Py()*d2->Py()+d2->Pz()*d2->Pz()+AliPID::ParticleMass(AliPID::kElectron)*AliPID::ParticleMass(AliPID::kElectron))); + // J/Psi 4-momentum vector + TLorentzVector motherMom=p1Mom+p2Mom; + + // boost all the 4-mom vectors to the mother rest frame + TVector3 beta = (-1.0/motherMom.E())*motherMom.Vect(); + p1Mom.Boost(beta); + p2Mom.Boost(beta); + projMom.Boost(beta); + targMom.Boost(beta); + + // x,y,z axes + TVector3 zAxis; + if(isHE) zAxis = (motherMom.Vect()).Unit(); + else zAxis = ((projMom.Vect()).Unit()-(targMom.Vect()).Unit()).Unit(); + TVector3 yAxis = ((projMom.Vect()).Cross(targMom.Vect())).Unit(); + TVector3 xAxis = (yAxis.Cross(zAxis)).Unit(); + + // return either theta or phi + if(isTheta) { + if(fD1.GetQ()>0) + return zAxis.Dot((p1Mom.Vect()).Unit()); + else + return zAxis.Dot((p2Mom.Vect()).Unit()); + } + else { + if(fD1.GetQ()>0) + return TMath::ATan2((p1Mom.Vect()).Dot(yAxis), (p1Mom.Vect()).Dot(xAxis)); + else + return TMath::ATan2((p2Mom.Vect()).Dot(yAxis), (p2Mom.Vect()).Dot(xAxis)); + } +} diff --git a/PWG3/dielectron/AliDielectronPair.h b/PWG3/dielectron/AliDielectronPair.h index fa5addd15eb..335c4511500 100644 --- a/PWG3/dielectron/AliDielectronPair.h +++ b/PWG3/dielectron/AliDielectronPair.h @@ -68,7 +68,10 @@ public: virtual Double_t M() const { return fPair.GetMass(); } virtual Double_t Eta() const { return fPair.GetEta();} - virtual Double_t Y() const { return TLorentzVector(Px(),Py(),Pz(),E()).Rapidity();} + virtual Double_t Y() const { + if((E()*E()-Px()*Px()-Py()*Py()-Pz()*Pz())>0.) return TLorentzVector(Px(),Py(),Pz(),E()).Rapidity(); + else return -1111.; + } virtual Short_t Charge() const { return fPair.GetQ();} virtual Int_t GetLabel() const { return fLabel; } @@ -89,6 +92,10 @@ public: Double_t DistanceDaughtersXY() const { return fD1.GetDistanceFromParticleXY(fD2); } Double_t DeviationDaughters() const { return fD1.GetDeviationFromParticle(fD2); } Double_t DeviationDaughtersXY() const { return fD1.GetDeviationFromParticleXY(fD2); } + // calculate cos(theta*) and phi* in HE and CS pictures + Double_t ThetaPhiCM(Bool_t isHE, Bool_t isTheta) const; + static Double_t ThetaPhiCM(const AliVParticle* d1, const AliVParticle* d2, + const Bool_t isHE, const Bool_t isTheta); // internal KF particle const AliKFParticle& GetKFParticle() const { return fPair; } @@ -96,6 +103,9 @@ public: const AliKFParticle& GetKFSecondDaughter() const { return fD2; } // daughter references + void SetRefFirstDaughter(AliVParticle * const track) {fRefD1 = track;} + void SetRefSecondDaughter(AliVParticle * const track) {fRefD2 = track;} + AliVParticle* GetFirstDaughter() const { return dynamic_cast(fRefD1.GetObject()); } AliVParticle* GetSecondDaughter() const { return dynamic_cast(fRefD2.GetObject()); } diff --git a/PWG3/dielectron/AliDielectronPairLegCuts.cxx b/PWG3/dielectron/AliDielectronPairLegCuts.cxx index 0ed9120e033..0c62499ca57 100644 --- a/PWG3/dielectron/AliDielectronPairLegCuts.cxx +++ b/PWG3/dielectron/AliDielectronPairLegCuts.cxx @@ -83,10 +83,16 @@ Bool_t AliDielectronPairLegCuts::IsSelected(TObject* track) Bool_t isLeg1selected=(fFilterLeg1.IsSelected(leg1)==selectedMaskLeg1); Bool_t isLeg2selected=(fFilterLeg2.IsSelected(leg2)==selectedMaskLeg2); + Bool_t isLeg1selectedMirror=(fFilterLeg1.IsSelected(leg2)==selectedMaskLeg1); + Bool_t isLeg2selectedMirror=(fFilterLeg2.IsSelected(leg1)==selectedMaskLeg2); + Bool_t isSelected=isLeg1selected&&isLeg2selected; if (fCutType==kAnyLeg) isSelected=isLeg1selected||isLeg2selected; - + + if (fCutType==kMixLegs) + isSelected=(isLeg1selected&&isLeg2selected)||(isLeg1selectedMirror&&isLeg2selectedMirror); + SetSelected(isSelected); return isSelected; } diff --git a/PWG3/dielectron/AliDielectronPairLegCuts.h b/PWG3/dielectron/AliDielectronPairLegCuts.h index c242a382dc2..2e10a8d3a3c 100644 --- a/PWG3/dielectron/AliDielectronPairLegCuts.h +++ b/PWG3/dielectron/AliDielectronPairLegCuts.h @@ -26,7 +26,7 @@ class AliDielectronPairLegCuts : public AliAnalysisCuts { public: - enum CutType { kBothLegs=0, kAnyLeg }; + enum CutType { kBothLegs=0, kAnyLeg, kMixLegs }; AliDielectronPairLegCuts(); AliDielectronPairLegCuts(const char* name, const char* title); diff --git a/PWG3/dielectron/AliDielectronSignalBase.cxx b/PWG3/dielectron/AliDielectronSignalBase.cxx index 44adbd31f18..8efaf1aa7fa 100644 --- a/PWG3/dielectron/AliDielectronSignalBase.cxx +++ b/PWG3/dielectron/AliDielectronSignalBase.cxx @@ -35,9 +35,9 @@ ClassImp(AliDielectronSignalBase) AliDielectronSignalBase::AliDielectronSignalBase() : TNamed(), - fValues(4), - fErrors(4), - fIntMin(2.99), + fValues(6), + fErrors(6), + fIntMin(3.01), fIntMax(3.19) { // @@ -49,9 +49,9 @@ AliDielectronSignalBase::AliDielectronSignalBase() : //______________________________________________ AliDielectronSignalBase::AliDielectronSignalBase(const char* name, const char* title) : TNamed(name, title), - fValues(4), - fErrors(4), - fIntMin(2.99), + fValues(6), + fErrors(6), + fIntMin(3.01), fIntMax(3.19) { // @@ -85,10 +85,14 @@ TPaveText* AliDielectronSignalBase::DrawStats(Double_t x1/*=0.*/, Double_t y1/*= t->SetFillColor(kWhite); t->SetBorderSize(1); t->SetTextAlign(12); - t->AddText(Form("Singal : %.2f #pm %.2f",GetSignal(),GetSignalError())); - t->AddText(Form("Backgnd: %.2f #pm %.2f",GetBackground(),GetBackgroundError())); - t->AddText(Form("Signif.: %.2f #pm %.2f",GetSignificance(),GetSignificanceError())); - t->AddText(Form("SoB : %.2f #pm %.2f",GetSignalOverBackground(),GetSignalOverBackgroundError())); + t->AddText(Form("Signal : %.5g #pm %.5g",GetSignal(),GetSignalError())); + t->AddText(Form("Backgnd: %.5g #pm %.5g",GetBackground(),GetBackgroundError())); + t->AddText(Form("Signif.: %.5g #pm %.5g",GetSignificance(),GetSignificanceError())); + t->AddText(Form("SoB : %.5g #pm %.5g",GetSignalOverBackground(),GetSignalOverBackgroundError())); + if (GetMass()>0){ + t->AddText(Form("Mass: %.5g #pm %.5g", GetMass(), GetMassError())); + t->AddText(Form("Mass res.: %.5g #pm %.5g", GetMassWidth(), GetMassWidthError())); + } t->Draw(); return t; diff --git a/PWG3/dielectron/AliDielectronSignalBase.h b/PWG3/dielectron/AliDielectronSignalBase.h index 88260809810..36dda37e6f3 100644 --- a/PWG3/dielectron/AliDielectronSignalBase.h +++ b/PWG3/dielectron/AliDielectronSignalBase.h @@ -42,7 +42,9 @@ public: void SetBackground(Double_t val, Double_t valErr) {fValues(1)=val; fErrors(1)=valErr;} void SetSignificance(Double_t val, Double_t valErr) {fValues(2)=val; fErrors(2)=valErr;} void SetSignalOverBackground(Double_t val, Double_t valErr) {fValues(3)=val; fErrors(3)=valErr;} - + void SetMass(Double_t val, Double_t valErr) {fValues(4)=val; fErrors(4)=valErr;} + void SetMassWidth(Double_t val, Double_t valErr) {fValues(5)=val; fErrors(5)=valErr;} + const TVectorD& GetValues() const {return fValues;} const TVectorD& GetErrors() const {return fErrors;} @@ -53,11 +55,15 @@ public: Double_t GetBackground() {return fValues(1);} Double_t GetSignificance() {return fValues(2);} Double_t GetSignalOverBackground() {return fValues(3);} + Double_t GetMass() {return fValues(4);} + Double_t GetMassWidth() {return fValues(5);} Double_t GetSignalError() {return fErrors(0);} Double_t GetBackgroundError() {return fErrors(1);} Double_t GetSignificanceError() {return fErrors(2);} Double_t GetSignalOverBackgroundError() {return fErrors(3);} + Double_t GetMassError() {return fErrors(4);} + Double_t GetMassWidthError() {return fValues(5);} void GetSignal(Double_t &val, Double_t &valErr) {val=fValues(0); valErr=fErrors(0);} void GetBackground(Double_t &val, Double_t &valErr) {val=fValues(1); valErr=fErrors(1);} @@ -75,6 +81,8 @@ public: 1: Background entries 2: Significance ( S/sqr(S+B) ) 3: Signal/Background + 4: Mass + 5: Mass width It is enough to calculate the signal and background and then call SetSignificanceAndSOB(TVectorD &values, TVectorD &errors) diff --git a/PWG3/dielectron/AliDielectronSignalExt.h b/PWG3/dielectron/AliDielectronSignalExt.h index 028bb10668e..ea36afa0c4c 100644 --- a/PWG3/dielectron/AliDielectronSignalExt.h +++ b/PWG3/dielectron/AliDielectronSignalExt.h @@ -37,6 +37,11 @@ public: void ProcessLS(TObjArray* const arrhist); // like-sign method void ProcessEM(TObjArray* const arrhist); // event mixing method + // getters + TH1F* GetHistogramSignal() const { return fSignal; } + TH1F* GetHistogramBackground() const { return fBackground; } + + // setters void SetMethod(Int_t method){ fMethod=method; } void SetHistograms(TH1F* const unlike, TH1F* const backg, TH1F* const signal); void SetRebin(Int_t rebin) {fRebin = rebin;} diff --git a/PWG3/dielectron/AliDielectronSignalFunc.cxx b/PWG3/dielectron/AliDielectronSignalFunc.cxx index a990a533646..b1c4e5b4a06 100644 --- a/PWG3/dielectron/AliDielectronSignalFunc.cxx +++ b/PWG3/dielectron/AliDielectronSignalFunc.cxx @@ -52,7 +52,8 @@ AliDielectronSignalFunc::AliDielectronSignalFunc() : fFitMin(2.5), fFitMax(4), fParM(-1), - fParMres(-1) + fParMres(-1), + fSignalHist(0x0) { // // Default Constructor @@ -72,7 +73,8 @@ AliDielectronSignalFunc::AliDielectronSignalFunc(const char* name, const char* t fFitMin(2.5), fFitMax(4), fParM(-1), - fParMres(-1) + fParMres(-1), + fSignalHist(0x0) { // // Named Constructor @@ -105,6 +107,9 @@ void AliDielectronSignalFunc::Process(TH1 * const hist) return; } + //set current signal hist pointer + fSignalHist=hist; + //initial the fit function to its default parameters if (fVInitParam.GetMatrixArray()) fSigBack->SetParameters(fVInitParam.GetMatrixArray()); @@ -153,7 +158,10 @@ void AliDielectronSignalFunc::Process(TH1 * const hist) SetSignal(signal,signal_err); SetBackground(background,background_err); SetSignificanceAndSOB(); - + if (fParM>-1){ + SetMass(fSigBack->GetParameter(fParM), fSigBack->GetParError(fParM)); + SetMassWidth(fSigBack->GetParameter(fParMres), fSigBack->GetParError(fParMres)); + } //cleanup delete hSignal; delete hBackground; @@ -295,7 +303,12 @@ void AliDielectronSignalFunc::Draw(const Option_t* option) fBackground->SetRange(fFitMin,fFitMax); if (!drawOpt.Contains("same")){ - grSig->Draw("af"); + if (fSignalHist){ + fSignalHist->Draw(); + grSig->Draw("f"); + } else { + grSig->Draw("af"); + } } else { grSig->Draw("f"); } @@ -303,10 +316,7 @@ void AliDielectronSignalFunc::Draw(const Option_t* option) fSigBack->Draw("same"); fBackground->Draw("same"); - if (optStat) { - TPaveText* pave=DrawStats(); - if (fParM>=0) pave->AddText(Form("Mass: %.3f #pm %.3f", fSigBack->GetParameter(fParM), fSigBack->GetParError(fParM))); - if (fParMres>=0) pave->AddText(Form("Mass res.: %.3f #pm %.3f", fSigBack->GetParameter(fParMres), fSigBack->GetParError(fParMres))); - } + if (optStat) DrawStats(); + } diff --git a/PWG3/dielectron/AliDielectronSignalFunc.h b/PWG3/dielectron/AliDielectronSignalFunc.h index 11f9da85c17..b3da5537c61 100644 --- a/PWG3/dielectron/AliDielectronSignalFunc.h +++ b/PWG3/dielectron/AliDielectronSignalFunc.h @@ -66,6 +66,8 @@ private: Int_t fParM; // Paramter which defines the mass Int_t fParMres; // Paramter which defines the resolution of the mass + + TH1 *fSignalHist; // Current signal histogram AliDielectronSignalFunc(const AliDielectronSignalFunc &c); AliDielectronSignalFunc &operator=(const AliDielectronSignalFunc &c); diff --git a/PWG3/dielectron/AliDielectronSpectrum.cxx b/PWG3/dielectron/AliDielectronSpectrum.cxx index 6d1221a2a36..31a73b412f8 100644 --- a/PWG3/dielectron/AliDielectronSpectrum.cxx +++ b/PWG3/dielectron/AliDielectronSpectrum.cxx @@ -29,6 +29,8 @@ Detailed description #include #include #include +#include +#include #include #include @@ -49,12 +51,15 @@ AliDielectronSpectrum::AliDielectronSpectrum() : fStepSignal(kTRUE), fStepSignificance(kFALSE), fStepSOB(kFALSE), + fStepMass(kFALSE), + fStepMassWidth(kFALSE), fSignalStep(-1), fCorrNom(-1), fCorrDenom(-1), fSignalMethods(0), fVariables("Pt"), fOwnerSpectrum(kTRUE), + fVisualDebug(kFALSE), fNvars(0), fVars(0x0), fNbins(0x0), @@ -77,12 +82,15 @@ AliDielectronSpectrum::AliDielectronSpectrum(const char* name, const char* title fStepSignal(kTRUE), fStepSignificance(kFALSE), fStepSOB(kFALSE), + fStepMass(kFALSE), + fStepMassWidth(kFALSE), fSignalStep(-1), fCorrNom(-1), fCorrDenom(-1), fSignalMethods(0), fVariables("Pt"), fOwnerSpectrum(kTRUE), + fVisualDebug(kFALSE), fNvars(0), fVars(0x0), fNbins(0x0), @@ -199,6 +207,8 @@ void AliDielectronSpectrum::CreateCFSpectrum() if (fStepSignal) nAddStep+=2; if (fStepSignificance) ++nAddStep; if (fStepSOB) ++nAddStep; + if (fStepMass) ++nAddStep; + if (fStepMassWidth) ++nAddStep; Int_t nStep=nAddStep*(fSignalMethods.GetEntries()); if (fSignalMethods.GetEntries()>1) nStep+=nAddStep; @@ -225,6 +235,12 @@ void AliDielectronSpectrum::CreateCFSpectrum() if (fStepSOB){ fCFSpectrum->SetStepTitle(steps++,(name+" (S/B)").Data()); } + if (fStepMass){ + fCFSpectrum->SetStepTitle(steps++,(name+" (Mass)").Data()); + } + if (fStepMassWidth){ + fCFSpectrum->SetStepTitle(steps++,(name+" (Mass width)").Data()); + } } if (fSignalMethods.GetEntries()>1){ @@ -302,6 +318,18 @@ void AliDielectronSpectrum::ExtractSignalInBins(Int_t variable) for (Int_t iMethod=0; iMethodProcess(&arrPairTypes); + if (fVisualDebug){ + TCanvas *c=(TCanvas*)gROOT->GetListOfCanvases()->FindObject("SpectrumVisualDebug"); + if (!c) c=new TCanvas("SpectrumVisualDebug","SpectrumVisualDebug"); + c->Clear(); + TString binsProc; + for (Int_t ivar=0; ivarDraw("stat"); + binsProc.Append(".png"); + binsProc.Prepend("SpectrumVisualDebug"); + c->Update(); + c->Print(binsProc.Data()); + } const TVectorD &val=signalMethod->GetValues(); const TVectorD &err=signalMethod->GetErrors(); @@ -340,6 +368,18 @@ void AliDielectronSpectrum::ExtractSignalInBins(Int_t variable) fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(3)); ++step; } + + if (fStepMass) { + fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(4)); + fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(4)); + ++step; + } + + if (fStepMassWidth) { + fCFSpectrum->GetGrid(step)->SetElement(fCurrentBins,val(5)); + fCFSpectrum->GetGrid(step)->SetElementError(fCurrentBins,err(5)); + ++step; + } }// end method loop arrPairTypes.Delete(); diff --git a/PWG3/dielectron/AliDielectronSpectrum.h b/PWG3/dielectron/AliDielectronSpectrum.h index 938ca9b80cb..265221f5d3a 100644 --- a/PWG3/dielectron/AliDielectronSpectrum.h +++ b/PWG3/dielectron/AliDielectronSpectrum.h @@ -46,10 +46,14 @@ public: void SetStepForSignal(Bool_t step=kTRUE) { fStepSignal=step; } void SetStepForSignificance(Bool_t step=kTRUE) { fStepSignificance=step; } - void SetStepForSignalOverBackground(Bool_t step=kTRUE) { fStepSOB=step; } - + void SetStepForSignalOverBackground(Bool_t step=kTRUE) { fStepSOB=step; } + void SetStepForMass(Bool_t step=kTRUE) { fStepMass=step; } + void SetStepForMassWidth(Bool_t step=kTRUE) { fStepMassWidth=step; } + void SetNoOwnerSpectrum(Bool_t noOwner=kTRUE) { fOwnerSpectrum=!noOwner; } - + + void SetVisualDebug(Bool_t vis=kTRUE) { fVisualDebug=vis; } + AliCFContainer* GetSpectrumContainer() const { return fCFSpectrum; } AliCFGridSparse* GetCorrectionMatrix() const { return fCFCorrMatrix; } @@ -65,7 +69,9 @@ private: Bool_t fStepSignal; // if to create a step for the signal Bool_t fStepSignificance; // if to create a step for the significance Bool_t fStepSOB; // if to create a step for signal over background - + Bool_t fStepMass; // if to create a setp for the mass + Bool_t fStepMassWidth; // if to create a setp for the mass width + Int_t fSignalStep; // step to use from the signal container Int_t fCorrNom; // Nominator to use from corr matrix container @@ -75,6 +81,8 @@ private: TString fVariables; // variable names as a function of which to extract the signal Bool_t fOwnerSpectrum; // if we own the creted spectrum + + Bool_t fVisualDebug; // if we want to show the fit and print it Int_t fNvars; //! number of variables Int_t *fVars; //! variable numbers translated from fVariables diff --git a/PWG3/dielectron/AliDielectronVarManager.cxx b/PWG3/dielectron/AliDielectronVarManager.cxx index 6f7c178fa79..5e2aae78471 100644 --- a/PWG3/dielectron/AliDielectronVarManager.cxx +++ b/PWG3/dielectron/AliDielectronVarManager.cxx @@ -50,10 +50,23 @@ const char* AliDielectronVarManager::fgkParticleNames[AliDielectronVarManager::k "NclsTRD", "TRDntracklets", "TRDpidQuality", + "TRDpidProb_Electrons", + "TRDpidProb_Pions", "ImpactParXY", "ImpactParZ", "TrackLength", "PdgCode", + + "PdgCodeMother", + + "NumberOfDaughters", + "HaveSameMother", + "ITS_signal", + "SSD1_signal", + "SSD2_signal", + "SDD1_signal", + "SDD2_signal", + "P_InnerParam", "TPC_signal", "TPC_nSigma_Electrons", @@ -61,6 +74,7 @@ const char* AliDielectronVarManager::fgkParticleNames[AliDielectronVarManager::k "TPC_nSigma_Muons", "TPC_nSigma_Kaons", "TPC_nSigma_Protons", + "TOF_nSigma_Electrons", "TOF_nSigma_Pions", "TOF_nSigma_Muons", "TOF_nSigma_Kaons", @@ -70,6 +84,10 @@ const char* AliDielectronVarManager::fgkParticleNames[AliDielectronVarManager::k "DecayLength", "R", "OpeningAngle", + "ThetaHE", + "PhiHE", + "ThetaCS", + "PhiCS", "LegDistance", "LegDistanceXY", "Merr", diff --git a/PWG3/dielectron/AliDielectronVarManager.h b/PWG3/dielectron/AliDielectronVarManager.h index 3a1eeb3b484..2abb87ff2ad 100644 --- a/PWG3/dielectron/AliDielectronVarManager.h +++ b/PWG3/dielectron/AliDielectronVarManager.h @@ -42,8 +42,10 @@ #include #include +#include #include "AliDielectronPair.h" +#include "AliDielectronMC.h" class AliVEvent; @@ -51,7 +53,7 @@ class AliVEvent; class AliDielectronVarManager : public TNamed { public: - + // Particle specific variables enum ValueTypes { kPx = 0, // px @@ -74,13 +76,28 @@ public: kNclsTPC, // number of clusters assigned in the TPC kNFclsTPC, // number of findable clusters in the TPC kTPCsignalN, // number of points used for dEdx + kNclsTRD, // number of clusters assigned in the TRD kTRDntracklets, // number of TRD tracklets used for tracking/PID TODO: correct getter kTRDpidQuality, // number of TRD tracklets used for PID + kTRDprobEle, // TRD electron pid probability + kTRDprobPio, // TRD electron pid probability + kImpactParXY, // Impact parameter in XY plane kImpactParZ, // Impact parameter in Z kTrackLength, // Track length kPdgCode, // PDG code + + kPdgCodeMother, // PDG code of the mother + + kNumberOfDaughters, // number of daughters + kHaveSameMother, // check that particles have the same mother (MC) + kITSsignal, // ITS dE/dx signal + kITSsignalSSD1, // SSD1 dE/dx signal + kITSsignalSSD2, // SSD2 dE/dx signal + kITSsignalSDD1, // SDD1 dE/dx signal + kITSsignalSDD2, // SDD2 dE/dx signal + kPIn, // momentum at inner wall of TPC (if available), used for PID kTPCsignal, // TPC dE/dx signal @@ -90,11 +107,12 @@ public: kTPCnSigmaKao, // number of sigmas to the dE/dx kaon line in the TPC kTPCnSigmaPro, // number of sigmas to the dE/dx proton line in the TPC + kTOFnSigmaEle, // number of sigmas to the pion line in the TOF kTOFnSigmaPio, // number of sigmas to the pion line in the TOF kTOFnSigmaMuo, // number of sigmas to the muon line in the TOF kTOFnSigmaKao, // number of sigmas to the kaon line in the TOF kTOFnSigmaPro, // number of sigmas to the proton line in the TOF - + kParticleMax, // // TODO: kRNClusters ?? // AliDielectronPair specific variables @@ -102,6 +120,13 @@ public: kDecayLength, // decay length kR, // distance to the origin kOpeningAngle, // opening angle + // helicity picture: Z-axis is considered the direction of the mother's 3-momentum vector + kThetaHE, // theta in mother's rest frame in the helicity picture + kPhiHE, // phi in mother's rest frame in the helicity picture + // Collins-Soper picture: Z-axis is considered the direction of the vectorial difference between + // the 3-mom vectors of target and projectile beams + kThetaCS, // theta in mother's rest frame in Collins-Soper picture + kPhiCS, // phi in mother's rest frame in Collins-Soper picture kLegDist, // distance of the legs kLegDistXY, // distance of the legs in XY kMerr, // error of mass calculation @@ -141,14 +166,14 @@ private: static void FillVarVParticle(const AliVParticle *particle, Double_t * const values); static void FillVarESDtrack(const AliESDtrack *particle, Double_t * const values); static void FillVarAODTrack(const AliAODTrack *particle, Double_t * const values); - static void FillVarMCParticle(const AliMCParticle *particle, Double_t * const values); + static void FillVarMCParticle(const AliMCParticle *particle, Double_t * const values); static void FillVarAODMCParticle(const AliAODMCParticle *particle, Double_t * const values); static void FillVarDielectronPair(const AliDielectronPair *pair, Double_t * const values); static void FillVarVEvent(const AliVEvent *event, Double_t * const values); static void FillVarESDEvent(const AliESDEvent *event, Double_t * const values); static void FillVarAODEvent(const AliAODEvent *event, Double_t * const values); static void FillVarMCEvent(const AliMCEvent *event, Double_t * const values); - + static AliESDpid* fgESDpid; // ESD pid object static AliVEvent* fgEvent; // current event pointer static AliKFVertex *fgKFVertex; // kf vertex @@ -202,7 +227,7 @@ inline void AliDielectronVarManager::FillVarVParticle(const AliVParticle *partic values[AliDielectronVarManager::kTheta] = particle->Theta(); values[AliDielectronVarManager::kEta] = particle->Eta(); values[AliDielectronVarManager::kY] = particle->Y(); - + values[AliDielectronVarManager::kE] = particle->E(); values[AliDielectronVarManager::kM] = particle->M(); values[AliDielectronVarManager::kCharge] = particle->Charge(); @@ -219,6 +244,7 @@ inline void AliDielectronVarManager::FillVarESDtrack(const AliESDtrack *particle // Fill common AliVParticle interface information FillVarVParticle(particle, values); + Double_t pidProbs[AliPID::kSPECIES]; // Fill AliESDtrack interface specific information values[AliDielectronVarManager::kNclsITS] = particle->GetNcls(0); // TODO: get rid of the plain numbers values[AliDielectronVarManager::kNclsTPC] = particle->GetNcls(1); // TODO: get rid of the plain numbers @@ -228,11 +254,48 @@ inline void AliDielectronVarManager::FillVarESDtrack(const AliESDtrack *particle values[AliDielectronVarManager::kTRDntracklets] = particle->GetTRDntracklets(); // TODO: GetTRDtracklets/GetTRDntracklets? values[AliDielectronVarManager::kTRDpidQuality] = particle->GetTRDpidQuality(); + //TRD pidProbs + particle->GetTRDpid(pidProbs); + values[AliDielectronVarManager::kTRDprobEle] = pidProbs[AliPID::kElectron]; + values[AliDielectronVarManager::kTRDprobPio] = pidProbs[AliPID::kPion]; + Float_t impactParXY, impactParZ; particle->GetImpactParameters(impactParXY, impactParZ); values[AliDielectronVarManager::kImpactParXY] = impactParXY; values[AliDielectronVarManager::kImpactParZ] = impactParZ; + + values[AliDielectronVarManager::kPdgCode]=0; + values[AliDielectronVarManager::kPdgCodeMother]=0; + + values[AliDielectronVarManager::kNumberOfDaughters]=-999; + + AliDielectronMC *mc=AliDielectronMC::Instance(); + + if (mc->HasMC()){ + if (mc->GetMCTrack(particle)) + values[AliDielectronVarManager::kPdgCode]= + mc->GetMCTrack(particle)->PdgCode(); + + Int_t pdgMother=mc->GetMotherPDG(particle); + if (pdgMother!=-999) + values[AliDielectronVarManager::kPdgCodeMother]=pdgMother; + + values[AliDielectronVarManager::kNumberOfDaughters]=mc->NumberOfDaughters(particle); + } //if(mc->HasMC()) + + + + values[AliDielectronVarManager::kITSsignal] = particle->GetITSsignal(); + + Double_t itsdEdx[4]; + particle->GetITSdEdxSamples(itsdEdx); + + values[AliDielectronVarManager::kITSsignalSSD1] = itsdEdx[0]; + values[AliDielectronVarManager::kITSsignalSSD2] = itsdEdx[1]; + values[AliDielectronVarManager::kITSsignalSDD1] = itsdEdx[2]; + values[AliDielectronVarManager::kITSsignalSDD2] = itsdEdx[3]; + values[AliDielectronVarManager::kTrackLength] = particle->GetIntegratedLength(); //dEdx information Double_t mom = particle->GetP(); @@ -250,6 +313,7 @@ inline void AliDielectronVarManager::FillVarESDtrack(const AliESDtrack *particle values[AliDielectronVarManager::kTPCnSigmaPro]=fgESDpid->NumberOfSigmasTPC(particle,AliPID::kProton); Double_t t0=fgESDpid->GetTOFResponse().GetTimeZero(); + values[AliDielectronVarManager::kTOFnSigmaEle]=fgESDpid->NumberOfSigmasTOF(particle,AliPID::kElectron,t0); values[AliDielectronVarManager::kTOFnSigmaPio]=fgESDpid->NumberOfSigmasTOF(particle,AliPID::kPion,t0); values[AliDielectronVarManager::kTOFnSigmaMuo]=fgESDpid->NumberOfSigmasTOF(particle,AliPID::kMuon,t0); values[AliDielectronVarManager::kTOFnSigmaKao]=fgESDpid->NumberOfSigmasTOF(particle,AliPID::kKaon,t0); @@ -273,6 +337,11 @@ inline void AliDielectronVarManager::FillVarAODTrack(const AliAODTrack *particle values[AliDielectronVarManager::kTRDntracklets] = 0; values[AliDielectronVarManager::kTRDpidQuality] = 0; + //TRD pidProbs + //TODO: set correctly + values[AliDielectronVarManager::kTRDprobEle] = 0; + values[AliDielectronVarManager::kTRDprobPio] = 0; + //TODO: This is only an approximation!!! values[AliDielectronVarManager::kTPCsignalN] = particle->GetTPCClusterMap().CountBits(); @@ -335,8 +404,20 @@ inline void AliDielectronVarManager::FillVarMCParticle(const AliMCParticle *part // Fill common AliVParticle interface information FillVarVParticle(particle, values); + + values[AliDielectronVarManager::kPdgCode]=0; + values[AliDielectronVarManager::kPdgCodeMother]=0; + + AliDielectronMC *mc=AliDielectronMC::Instance(); + // Fill AliMCParticle interface specific information values[AliDielectronVarManager::kPdgCode] = particle->PdgCode(); + + AliMCParticle *mother = mc->GetMCTrackMother(particle); + if (mother) values[AliDielectronVarManager::kPdgCodeMother] = mother->PdgCode(); + + + values[AliDielectronVarManager::kNumberOfDaughters]=mc->NumberOfDaughters(particle); } inline void AliDielectronVarManager::FillVarAODMCParticle(const AliAODMCParticle *particle, Double_t * const values) @@ -348,8 +429,20 @@ inline void AliDielectronVarManager::FillVarAODMCParticle(const AliAODMCParticle // Fill common AliVParticle interface information FillVarVParticle(particle, values); + + values[AliDielectronVarManager::kPdgCode]=0; + values[AliDielectronVarManager::kPdgCodeMother]=0; + + AliDielectronMC *mc=AliDielectronMC::Instance(); + + // Fill AliAODMCParticle interface specific information - values[AliDielectronVarManager::kPdgCode] = particle->GetPdgCode(); + values[AliDielectronVarManager::kPdgCode] = particle->PdgCode(); + + AliVParticle *mother = mc->GetMCTrackMother(particle); + if (mother) values[AliDielectronVarManager::kPdgCodeMother] = mother->PdgCode(); + + values[AliDielectronVarManager::kNumberOfDaughters]=mc->NumberOfDaughters(particle); } inline void AliDielectronVarManager::FillVarDielectronPair(const AliDielectronPair *pair, Double_t * const values) @@ -368,10 +461,25 @@ inline void AliDielectronVarManager::FillVarDielectronPair(const AliDielectronPa values[AliDielectronVarManager::kDecayLength] = kfPair.GetDecayLength(); values[AliDielectronVarManager::kR] = kfPair.GetR(); values[AliDielectronVarManager::kOpeningAngle] = pair->OpeningAngle(); + values[AliDielectronVarManager::kThetaHE] = pair->ThetaPhiCM(kTRUE, kTRUE); + values[AliDielectronVarManager::kPhiHE] = pair->ThetaPhiCM(kTRUE, kFALSE); + values[AliDielectronVarManager::kThetaCS] = pair->ThetaPhiCM(kFALSE, kTRUE); + values[AliDielectronVarManager::kPhiCS] = pair->ThetaPhiCM(kFALSE, kFALSE); values[AliDielectronVarManager::kLegDist] = pair->DistanceDaughters(); values[AliDielectronVarManager::kLegDistXY] = pair->DistanceDaughtersXY(); - values[AliDielectronVarManager::kMerr] = kfPair.GetErrMass()>1e-30&&kfPair.GetErrMass()>1e-30?kfPair.GetErrMass()/kfPair.GetMass():1000000; + values[AliDielectronVarManager::kMerr] = kfPair.GetErrMass()>1e-30&&kfPair.GetMass()>1e-30?kfPair.GetErrMass()/kfPair.GetMass():1000000; values[AliDielectronVarManager::kPairType] = pair->GetType(); + + + + AliDielectronMC *mc=AliDielectronMC::Instance(); + + if (mc->HasMC()){ + Bool_t samemother = mc->HaveSameMother(pair); + values[AliDielectronVarManager::kHaveSameMother] = samemother ; + }//if (mc->HasMC()) + + } @@ -381,13 +489,23 @@ inline void AliDielectronVarManager::FillVarVEvent(const AliVEvent *event, Doubl // Fill event information available for histogramming into an array // const AliVVertex *primVtx = event->GetPrimaryVertex(); + + values[AliDielectronVarManager::kXvPrim] = 0; + values[AliDielectronVarManager::kYvPrim] = 0; + values[AliDielectronVarManager::kZvPrim] = 0; + values[AliDielectronVarManager::kChi2NDF] = 0; + + values[AliDielectronVarManager::kNTrk] = 0; + values[AliDielectronVarManager::kNevents] = 0; //always fill bin 0; + + if (!primVtx) return; + values[AliDielectronVarManager::kXvPrim] = primVtx->GetX(); values[AliDielectronVarManager::kYvPrim] = primVtx->GetY(); values[AliDielectronVarManager::kZvPrim] = primVtx->GetZ(); values[AliDielectronVarManager::kChi2NDF] = primVtx->GetChi2perNDF(); values[AliDielectronVarManager::kNTrk] = event->GetNumberOfTracks(); - values[AliDielectronVarManager::kNevents] = 0; //always fill bin 0; } inline void AliDielectronVarManager::FillVarESDEvent(const AliESDEvent *event, Double_t * const values) @@ -445,6 +563,7 @@ inline void AliDielectronVarManager::InitESDpid(Int_t type) alephParameters[2] = 3.40030e-09; alephParameters[3] = 1.96178e+00; alephParameters[4] = 3.91720e+00; + fgESDpid->GetTOFResponse().SetTimeResolution(80.); // data if (type==1){ @@ -453,12 +572,15 @@ inline void AliDielectronVarManager::InitESDpid(Int_t type) alephParameters[2] = 5.04114e-11; alephParameters[3] = 2.12543e+00; alephParameters[4] = 4.88663e+00; + fgESDpid->GetTOFResponse().SetTimeResolution(130.); + fgESDpid->GetTPCResponse().SetMip(47.9); } fgESDpid->GetTPCResponse().SetBetheBlochParameters( alephParameters[0],alephParameters[1],alephParameters[2], alephParameters[3],alephParameters[4]); + fgESDpid->GetTPCResponse().SetSigma(3.79301e-03, 2.21280e+04); } @@ -468,7 +590,7 @@ inline void AliDielectronVarManager::SetEvent(AliVEvent * const ev) fgEvent = ev; if (fgKFVertex) delete fgKFVertex; fgKFVertex=0x0; - if (ev) fgKFVertex=new AliKFVertex(*ev->GetPrimaryVertex()); + if (ev && ev->GetPrimaryVertex()) fgKFVertex=new AliKFVertex(*ev->GetPrimaryVertex()); } /* inline void AliDielectronVarManager::FillValues(const TParticle *particle, Double_t *values) diff --git a/PWG3/dielectron/AddTaskJPSI.C b/PWG3/dielectron/macros/AddTaskJPSI.C similarity index 95% rename from PWG3/dielectron/AddTaskJPSI.C rename to PWG3/dielectron/macros/AddTaskJPSI.C index 9dc535cd968..fdecb598ee9 100644 --- a/PWG3/dielectron/AddTaskJPSI.C +++ b/PWG3/dielectron/macros/AddTaskJPSI.C @@ -32,9 +32,6 @@ AliAnalysisTask *AddTaskJPSI(const char* config=""){ //---------------------- //create data containers //---------------------- - - //find input container - AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer(); TString containerName = mgr->GetCommonFileName(); containerName += ":PWG3_dielectron"; diff --git a/PWG3/dielectron/macros/AddTaskJPSIFilter.C b/PWG3/dielectron/macros/AddTaskJPSIFilter.C new file mode 100644 index 00000000000..3b48a8f1f95 --- /dev/null +++ b/PWG3/dielectron/macros/AddTaskJPSIFilter.C @@ -0,0 +1,54 @@ +AliAnalysisTask *AddTaskJPSIFilter(){ + //get the current analysis manager + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error("AddTask_jpsi_DielectronFilter", "No analysis manager found."); + return 0; + } + // currently don't accept AOD input + if (mgr->GetInputEventHandler()->IsA()==AliAODInputHandler::Class()) { + Warning("AddTask_jpsi_JPsi","No AOD input supported currently. Not adding the task!"); + return 0; + } + //check for output aod handler + if (!mgr->GetOutputEventHandler()||mgr->GetOutputEventHandler()->IsA()!=AliAODHandler::Class()) { + Warning("AddTask_jpsi_DielectronFilter","No AOD output handler available. Not adding the task!"); + return 0; + } + + //Do we have an MC handler? + Bool_t hasMC=(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()!=0x0); + + //Create task and add it to the analysis manager + AliAnalysisTaskDielectronFilter *task=new AliAnalysisTaskDielectronFilter("jpsi_DielectronFilter"); + + gROOT->LoadMacro("$ALICE_ROOT/PWG3/dielectron/macros/ConfigJpsi2eeFilter.C"); + AliDielectron *jpsi=ConfigJpsi2eeFilter(); + if (!hasMC) task->UsePhysicsSelection(); + task->SetDielectron(jpsi); + mgr->AddTask(task); + + //---------------------- + //create data containers + //---------------------- + + + TString containerName = mgr->GetCommonFileName(); + containerName += ":PWG3_dielectronFilter"; + + //create output container + + AliAnalysisDataContainer *cOutputHist1 = + mgr->CreateContainer("QA", TList::Class(), AliAnalysisManager::kOutputContainer, + containerName.Data()); + + AliAnalysisDataContainer *cOutputHist2 = + mgr->CreateContainer("CF", TList::Class(), AliAnalysisManager::kOutputContainer, + containerName.Data()); + + mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task, 1, cOutputHist1); + mgr->ConnectOutput(task, 2, cOutputHist2); + + return task; +} diff --git a/PWG3/dielectron/macros/ConfigJpsi2eeData.C b/PWG3/dielectron/macros/ConfigJpsi2eeData.C index 35a5a83dba6..17886d7a234 100644 --- a/PWG3/dielectron/macros/ConfigJpsi2eeData.C +++ b/PWG3/dielectron/macros/ConfigJpsi2eeData.C @@ -6,8 +6,14 @@ void SetupTrackCuts(AliDielectron *die, Int_t cutDefinition); void SetupPairCuts(AliDielectron *die, Int_t cutDefinition); AliESDtrackCuts *SetupESDtrackCuts(Int_t cutDefinition); - -TString names=("trackQ+Pt>1GeV+loose_PID;trackQ+Pt>1GeV+tight_PID;trackQ+Pt>.5GeV"); +/* +trackQ+Pt>1GeV+|TPCnSigE|<3 +trackQ+PID1 +trackQ+PID6 +trackQ+PID7 +trackQ+Pt>.5GeV +*/ +TString names=("trackQ+Pt>1GeV+|TPCnSigE|<3;trackQ+PID1;trackQ+PID6;trackQ+PID7;trackQ+Pt>.5GeV"); TObjArray *arrNames=names.Tokenize(";"); const Int_t nDie=arrNames->GetEntries(); @@ -39,7 +45,7 @@ AliDielectron* ConfigJpsi2ee(Int_t cutDefinition) InitHistograms(die,cutDefinition); // the last definition uses no cuts and only the QA histograms should be filled! - if (cutDefinition1+|TPCnSigE|<3","Pt>1+|TPCnSigE|<3"); + pt->AddCut(AliDielectronVarManager::kPt,1.,1e30); + pt->AddCut(AliDielectronVarManager::kTPCnSigmaEle, -3., 3.); + die->GetTrackFilter().AddCuts(pt); + } - AliDielectronVarCuts *pt = new AliDielectronVarCuts("Pt>1","Pt>1"); - pt->AddCut(AliDielectronVarManager::kPt,1.,1e30); - die->GetTrackFilter().AddCuts(pt); - - //loose PID + //PID 1 if (cutDefinition==0){ - //ESD pid cuts (TPC nSigma Electrons) - AliDielectronVarCuts *pid = new AliDielectronVarCuts("TPCpid","TPC nSigma cut |e|<3 pi>2 |P|>2 |K|>2"); - //include - pid->AddCut(AliDielectronVarManager::kTPCnSigmaEle, -3., 3.); - //exclude - pid->AddCut(AliDielectronVarManager::kTPCnSigmaPio, -20., 2., kTRUE); - pid->AddCut(AliDielectronVarManager::kTPCnSigmaPro, -2., 2., kTRUE); - pid->AddCut(AliDielectronVarManager::kTPCnSigmaKao, -2., 2., kTRUE); + AliDielectronPID *pid=new AliDielectronPID("pid1","pid1"); + pid->SetDefaults(1); die->GetTrackFilter().AddCuts(pid); } - //tight PID + //PID 6 if (cutDefinition==1){ - //ESD pid cuts (TPC nSigma Electrons) - AliDielectronVarCuts *pid = new AliDielectronVarCuts("TPCpid","TPC nSigma cut |e|<2 pi>2 |P|>2 |K|>2; |dXY|<200um"); - //include - pid->AddCut(AliDielectronVarManager::kTPCnSigmaEle, -2., 2.); - //exclude - pid->AddCut(AliDielectronVarManager::kTPCnSigmaPio, -2., 2., kTRUE); - pid->AddCut(AliDielectronVarManager::kTPCnSigmaPro, -2., 2., kTRUE); - pid->AddCut(AliDielectronVarManager::kTPCnSigmaKao, -2., 2., kTRUE); - //tighten Impact parameter to 200um - pid->AddCut(AliDielectronVarManager::kImpactParXY, -0.02, 0.02); + AliDielectronPID *pid=new AliDielectronPID("pid6","pid6"); + pid->SetDefaults(6); + die->GetTrackFilter().AddCuts(pid); + } + + //PID 7 + if (cutDefinition==2){ + AliDielectronPID *pid=new AliDielectronPID("pid7","pid7"); + pid->SetDefaults(7); die->GetTrackFilter().AddCuts(pid); } @@ -131,7 +133,7 @@ AliESDtrackCuts *SetupESDtrackCuts(Int_t cutDefinition) esdTrackCuts->SetAcceptKinkDaughters(kFALSE); esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst); - esdTrackCuts->SetMinNClustersTPC(75); + esdTrackCuts->SetMinNClustersTPC(110); esdTrackCuts->SetMaxChi2PerClusterTPC(4); return esdTrackCuts; diff --git a/PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C b/PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C index ca430af6bba..6b897cbd77e 100644 --- a/PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C +++ b/PWG3/dielectron/macros/ConfigJpsi2eeDataAOD.C @@ -49,8 +49,6 @@ void SetupTrackCuts(AliDielectron *die, Int_t cutDefinition) // Setup the track cuts // - //ESD quality cuts - die->GetTrackFilter().AddCuts(SetupESDtrackCuts(cutDefinition)); //QA no CF diff --git a/PWG3/dielectron/macros/ConfigJpsi2eeEff.C b/PWG3/dielectron/macros/ConfigJpsi2eeEff.C new file mode 100644 index 00000000000..53bd71ed3bf --- /dev/null +++ b/PWG3/dielectron/macros/ConfigJpsi2eeEff.C @@ -0,0 +1,254 @@ + +void InitHistograms(AliDielectron *die, Int_t cutDefinition); +void InitCF(AliDielectron* die, Int_t cutDefinition); + +void SetupTrackCuts(AliDielectron *die, Int_t cutDefinition); +void SetupPairCuts(AliDielectron *die, Int_t cutDefinition); + +AliESDtrackCuts *SetupESDtrackCuts(); + +TString names=("no_cuts;track_quality;tq+mcPIDele;tq+4#sigma_dEdx_E"); +TObjArray *arrNames=names.Tokenize(";"); +const Int_t nDie=arrNames->GetEntriesFast(); + +AliDielectron* ConfigJpsi2ee(Int_t cutDefinition=1) +{ + // + // Setup the instance of AliDielectron + // + + // create the actual framework object + TString name=Form("%02d",cutDefinition); + if (cutDefinitionGetEntriesFast()){ + name=arrNames->At(cutDefinition)->GetName(); + } + AliDielectron *die = + new AliDielectron(Form("%s",name.Data()), + Form("Track cuts: %s",name.Data())); + + // cut setup + SetupTrackCuts(die,cutDefinition); + SetupPairCuts(die,cutDefinition); + + // + // histogram setup + // only if an AliDielectronHistos object is attached to the + // dielectron framework histograms will be filled + // + InitHistograms(die,cutDefinition); + InitCF(die,cutDefinition); + + return die; +} + +//______________________________________________________________________________________ +void SetupTrackCuts(AliDielectron *die, Int_t cutDefinition) +{ + // + // Setup the track cuts + // + + //ESD quality cuts + if (cutDefinition>0){ + die->GetTrackFilter().AddCuts(SetupESDtrackCuts()); + } + + //MC pid + if (cutDefinition==2){ + //MC pid + AliDielectronVarCuts *mcpid=new AliDielectronVarCuts("legMCpid","Leg MC pid"); + mcpid->SetCutOnMCtruth(); + mcpid->SetCutType(AliDielectronVarCuts::kAny); + mcpid->AddCut(AliDielectronVarManager::kPdgCode, 11); + mcpid->AddCut(AliDielectronVarManager::kPdgCode, -11); + die->GetTrackFilter().AddCuts(mcpid); + } + + if (cutDefinition>=3){ + //ESD pid cuts (TPC nSigma) + AliDielectronVarCuts *pid = new AliDielectronVarCuts("TPCnSigma","TPC nSigma cut"); + pid->AddCut(AliDielectronVarManager::kTPCnSigmaEle, -4., 4.); + die->GetTrackFilter().AddCuts(pid); + } + +} + +//______________________________________________________________________________________ +void SetupPairCuts(AliDielectron *die, Int_t cutDefinition) +{ + // + // Setup the pair cuts + // + + + if (cutDefinition>0){ + //MC mother in y +-.9 + AliDielectronVarCuts *etaMC=new AliDielectronVarCuts("|MC Y|<.9","|MC Y|<.9"); + etaMC->AddCut(AliDielectronVarManager::kY,-.8,.8); + etaMC->SetCutOnMCtruth(); + die->GetPairFilter().AddCuts(etaMC); + + //reject conversions + AliDielectronVarCuts *openingAngleCut=new AliDielectronVarCuts("OpeningAngle","Opening angle > .35rad"); + openingAngleCut->AddCut(AliDielectronVarManager::kOpeningAngle,.35,4.); + die->GetPairFilter().AddCuts(openingAngleCut); + + } + + if (cutDefinition==3){ + // + //ESD pid cuts (TPC nSigma Pions) + // + AliDielectronPairLegCuts *tpcNSigmaPi = new AliDielectronPairLegCuts("|n#sigma#pi|>2","|n#sigma#pi|>2"); + + AliDielectronVarCuts *pidPiM = new AliDielectronVarCuts("TPCnSigmaPi-","TPC nSigma- pi cut"); + pidPiM->AddCut(AliDielectronVarManager::kTPCnSigmaPio, -10., -2.); + + AliDielectronVarCuts *pidPiP = new AliDielectronVarCuts("TPCnSigmaPi+","TPC nSigma+ pi cut"); + pidPiP->AddCut(AliDielectronVarManager::kTPCnSigmaPio, 2., 10.); + + tpcNSigmaPi->GetLeg1Filter().AddCuts(pidPiM); + tpcNSigmaPi->GetLeg1Filter().AddCuts(pidPiP); + + tpcNSigmaPi->GetLeg2Filter().AddCuts(pidPiM); + tpcNSigmaPi->GetLeg2Filter().AddCuts(pidPiP); + + die->GetPairFilter().AddCuts(tpcNSigmaPi); + + // + //TRD pid quality cuts + // + AliDielectronVarCuts *trdNtrklet = new AliDielectronVarCuts("TRDpidQuality","TRD pid quality"); + trdNtrklet->AddCut(AliDielectronVarManager::kTRDpidQuality, 3.5, 8.); + + AliDielectronPairLegCuts *trdNtrkletAny = new AliDielectronPairLegCuts("TRDntrlt>=4 any","TRDntrlt>=4 any"); + trdNtrkletAny->GetLeg1Filter().AddCuts(trdNtrklet); + trdNtrkletAny->GetLeg2Filter().AddCuts(trdNtrklet); + trdNtrkletAny->SetCutType(AliDielectronPairLegCuts::kAnyLeg); + die->GetPairFilter().AddCuts(trdNtrkletAny); + + AliDielectronPairLegCuts *trdNtrkletBoth = new AliDielectronPairLegCuts("TRDntrlt>=4 both","TRDntrlt>=4 both"); + trdNtrkletBoth->GetLeg1Filter().AddCuts(trdNtrklet); + trdNtrkletBoth->GetLeg2Filter().AddCuts(trdNtrklet); + die->GetPairFilter().AddCuts(trdNtrkletBoth); + } +} + +//______________________________________________________________________________________ +AliESDtrackCuts *SetupESDtrackCuts() +{ + // + // Setup default AliESDtrackCuts + // + AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts; + + esdTrackCuts->SetMaxDCAToVertexZ(3.0); + esdTrackCuts->SetMaxDCAToVertexXY(1); + esdTrackCuts->SetRequireTPCRefit(kTRUE); + esdTrackCuts->SetRequireITSRefit(kTRUE); + esdTrackCuts->SetAcceptKinkDaughters(kFALSE); + esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); + + esdTrackCuts->SetMinNClustersTPC(75); + esdTrackCuts->SetMaxChi2PerClusterTPC(4); + + return esdTrackCuts; +} + +//______________________________________________________________________________________ +void InitHistograms(AliDielectron *die, Int_t cutDefinition) +{ + // + // Initialise the histograms + // + + //Setup histogram classes + AliDielectronHistos *histos= + new AliDielectronHistos(die->GetName(), + die->GetTitle()); + + //Initialise histogram classes + histos->SetReservedWords("Track;Pair"); + + //Event class + histos->AddClass("Event"); + + //Track classes + //to fill also track info from 2nd event loop until 2 + for (Int_t i=0; i<2; ++i){ + histos->AddClass(Form("Track_%s",AliDielectron::TrackClassName(i))); + } + + //Pair classes + // to fill also mixed event histograms loop until 10 + for (Int_t i=0; i<3; ++i){ + histos->AddClass(Form("Pair_%s",AliDielectron::PairClassName(i))); + } + //add histograms to event class + histos->UserHistogram("Event","nEvents","Number of processed events after cuts;Number events", + 1,0.,1.,AliDielectronVarManager::kNevents); + + //add histograms to Track classes + histos->UserHistogram("Track","Pt","Pt;Pt [GeV];#tracks",200,0,20.,AliDielectronVarManager::kPt); + histos->UserHistogram("Track","dEdx_P","dEdx;P [GeV];TPC signal (arb units);#tracks", + 400,1e-2,20.,200,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal,kTRUE); + histos->UserHistogram("Track","TPCnSigma_P","TPC number of sigmas Electrons;P [GeV];TPC number of sigmas Electrons;#tracks", + 400,1e-2,20.,100,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle,kTRUE); + histos->UserHistogram("Track","dXY","dXY;dXY [cm];#tracks",200,-2.,2.,AliDielectronVarManager::kImpactParXY); + + //add histograms to Pair classes + histos->UserHistogram("Pair","InvMass","Inv.Mass;Inv. Mass [GeV];#pairs", + 500,0.,4.,AliDielectronVarManager::kM); + histos->UserHistogram("Pair","Rapidity","Rapidity;Rapidity;#pairs", + 100,-2.,2.,AliDielectronVarManager::kY); + histos->UserHistogram("Pair","OpeningAngle","Opening angle;angle", + 100,0.,3.15,AliDielectronVarManager::kOpeningAngle); + histos->UserHistogram("Pair","InvMass_OpeningAngle","Opening angle:Inv.Mass;Inv. Mass [GeV];angle", + 100,0.,4.,100,0.,3.15,AliDielectronVarManager::kM,AliDielectronVarManager::kOpeningAngle); + + die->SetHistogramManager(histos); +} + +void InitCF(AliDielectron* die, Int_t cutDefinition) +{ + // + // Setupd the CF Manager if needed + // + + AliDielectronCF *cf=new AliDielectronCF(die->GetName(),die->GetTitle()); + + //pair variables + cf->AddVariable(AliDielectronVarManager::kPt,50,0,10); + cf->AddVariable(AliDielectronVarManager::kEta,40,-2,2); + cf->AddVariable(AliDielectronVarManager::kY,40,-2,2); + cf->AddVariable(AliDielectronVarManager::kM,100,0,4); + cf->AddVariable(AliDielectronVarManager::kPairType,10,0,10); + cf->AddVariable(AliDielectronVarManager::kOpeningAngle,10,0,3.15); + //leg variables + cf->AddVariable(AliDielectronVarManager::kEta,40,-2,2,kTRUE); + cf->AddVariable(AliDielectronVarManager::kImpactParXY,50,-.1,.1,kTRUE); + cf->AddVariable(AliDielectronVarManager::kTPCnSigmaEle,12,-3,3,kTRUE); + + //no cuts + //only in this case write MC truth info + if (cutDefinition==0){ + cf->SetStepForMCtruth(); + cf->SetStepForNoCutsMCmotherPid(); + cf->SetStepForAfterAllCuts(kFALSE); + } + + if (cutDefinition==1){ + cf->SetStepForNoCutsMCmotherPid(); + } + + if (cutDefinition>1){ + cf->SetStepsForEachCut(); + } + + if (cutDefinition>2){ + cf->SetStepsForCutsIncreasing(); + } + + die->SetCFManagerPair(cf); + +} diff --git a/PWG3/dielectron/macros/ConfigJpsi2eeFilter.C b/PWG3/dielectron/macros/ConfigJpsi2eeFilter.C new file mode 100644 index 00000000000..543bb75767a --- /dev/null +++ b/PWG3/dielectron/macros/ConfigJpsi2eeFilter.C @@ -0,0 +1,156 @@ + +void InitHistograms(AliDielectron *die); +void InitCF(AliDielectron* die); + +void SetupTrackCuts(AliDielectron *die); +void SetupPairCuts(AliDielectron *die); + +AliESDtrackCuts *SetupESDtrackCuts(); + +AliDielectron* ConfigJpsi2eeFilter() +{ + // + // Setup the instance of AliDielectron + // + + // create the actual framework object + TString name="trackQ+Pt>0.5+60GetTrackFilter().AddCuts(SetupESDtrackCuts()); + + //Pt cut + AliDielectronVarCuts *pt = new AliDielectronVarCuts("Pt>.5+60.5 && 60AddCut(AliDielectronVarManager::kPt,.5,1e30); + pt->AddCut(AliDielectronVarManager::kTPCsignal,60.,100.); + + die->GetTrackFilter().AddCuts(pt); +} + +//______________________________________________________________________________________ +void SetupPairCuts(AliDielectron *die) +{ + // + // Setup the pair cuts + // + + + //Invarian mass selection + AliDielectronVarCuts *invMassCut=new AliDielectronVarCuts("InvMass","2AddCut(AliDielectronVarManager::kM,2.,4.); +// invMassCut->AddCut(AliDielectronVarManager::kPairType,1.); + die->GetPairFilter().AddCuts(invMassCut); + +} + +//______________________________________________________________________________________ +AliESDtrackCuts *SetupESDtrackCuts() +{ + // + // Setup default AliESDtrackCuts + // + AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts; + + esdTrackCuts->SetMaxDCAToVertexZ(3.0); + esdTrackCuts->SetMaxDCAToVertexXY(1.0); + esdTrackCuts->SetRequireTPCRefit(kTRUE); + esdTrackCuts->SetRequireITSRefit(kTRUE); + esdTrackCuts->SetAcceptKinkDaughters(kFALSE); + esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst); + + esdTrackCuts->SetMinNClustersTPC(100); + esdTrackCuts->SetMaxChi2PerClusterTPC(4); + + return esdTrackCuts; +} + + +//______________________________________________________________________________________ +void InitHistograms(AliDielectron *die) +{ + // + // Initialise the histograms + // + +//Setup histogram classes + AliDielectronHistos *histos= + new AliDielectronHistos(die->GetName(), + die->GetTitle()); + + //Initialise histogram classes + histos->SetReservedWords("Track;Pair"); + + //Event class + histos->AddClass("Event"); + + //Track classes + //to fill also track info from 2nd event loop until 2 + for (Int_t i=0; i<2; ++i){ + histos->AddClass(Form("Track_%s",AliDielectron::TrackClassName(i))); + } + + //Pair classes + // to fill also mixed event histograms loop until 10 + for (Int_t i=0; i<3; ++i){ + histos->AddClass(Form("Pair_%s",AliDielectron::PairClassName(i))); + } + + //add histograms to event class + histos->UserHistogram("Event","nEvents","Number of processed events after cuts;Number events", + 1,0.,1.,AliDielectronVarManager::kNevents); + + //add histograms to Track classes + histos->UserHistogram("Track","Pt","Pt;Pt [GeV];#tracks",200,0,20.,AliDielectronVarManager::kPt); + histos->UserHistogram("Track","dEdx_P","dEdx;P [GeV];TPC signal (arb units);#tracks", + 400,0.2,20.,200,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal,kTRUE); + histos->UserHistogram("Track","TPCnSigmaEle_P","TPC number of sigmas Electrons;P [GeV];TPC number of sigmas Electrons;#tracks", + 400,0.2,20.,100,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle,kTRUE); + histos->UserHistogram("Track","TPCnSigmaKao_P","TPC number of sigmas Kaons;P [GeV];TPC number of sigmas Kaons;#tracks", + 400,0.2,20.,100,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaKao,kTRUE); + histos->UserHistogram("Track","TPCnSigmaPio_P","TPC number of sigmas Pions;P [GeV];TPC number of sigmas Pions;#tracks", + 400,0.2,20.,100,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaPio,kTRUE); + + + histos->UserHistogram("Track","TOFnSigmaKao_P","TOF number of sigmas Kaons;P [GeV];TOF number of sigmas Kaons;#tracks", + 400,0.2,20.,100,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFnSigmaKao,kTRUE); + histos->UserHistogram("Track","TOFnSigmaPro_P","TOF number of sigmas Protons;P [GeV];TOF number of sigmas Protons;#tracks", + 400,0.2,20.,100,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFnSigmaPro,kTRUE); + + histos->UserHistogram("Track","dXY","dXY;dXY [cm];#tracks",200,-2.,2.,AliDielectronVarManager::kImpactParXY); + + histos->UserHistogram("Track","TPCnCls","Number of Clusters TPC;TPC number clusteres;#tracks",159,0.,159.,AliDielectronVarManager::kNclsTPC); + + //add histograms to Pair classes + histos->UserHistogram("Pair","InvMass","Inv.Mass;Inv. Mass [GeV];#pairs", + 201,-.01,4.01,AliDielectronVarManager::kM); + histos->UserHistogram("Pair","Rapidity","Rapidity;Rapidity;#pairs", + 100,-2.,2.,AliDielectronVarManager::kY); + histos->UserHistogram("Pair","OpeningAngle","Opening angle;angle", + 100,0.,3.15,AliDielectronVarManager::kOpeningAngle); + histos->UserHistogram("Pair","Chi2/NDF","#Chi^{2}/NDF;#Chi^{2}/NDF", + 100, 0., 20., AliDielectronVarManager::kChi2NDF); + + die->SetHistogramManager(histos); +} diff --git a/PWG3/dielectron/macros/MakeDataReport.C b/PWG3/dielectron/macros/MakeDataReport.C new file mode 100644 index 00000000000..eb9440f234f --- /dev/null +++ b/PWG3/dielectron/macros/MakeDataReport.C @@ -0,0 +1,708 @@ +void SetupStyle(); + +void MakeDataReport(const char* outputFile="JpsiDataReport.pdf", + const char* histos="jpsi_HistosSE.root", + const char* cf="jpsi_CF.root") +{ + // + // Make a pdf file with the efficiency report + // + + SetupStyle(); + + AliDielectronCFdraw d(cf); + d.SetRangeUser("PairType",1,1); + d.SetRangeUser("Y",-.89,.9,"0"); + + + TFile f("jpsi_HistosSE.root"); + + AliDielectronHistos h; + TIter nextHists((TList*)f.Get("Dielectron_Histos")); + + TPaveText pt(.02,.6,.98,.8); + TText *t1=pt.AddText(""); + TText *t2=pt.AddText(""); + + TCanvas *c1=new TCanvas; + + TPDF p(outputFile); + + // + // Invariant mass plots + // + + + // + // Make QA info + // + + t1->SetTitle("QA summary plots for"); + THashList *list=0x0; + while ( (list=(THashList*)nextHists()) ){ + h.SetHistogramList(*list); + t2->SetTitle(list->GetName()); + pt.Draw(); + c1->Update(); + h.Draw(); + c1->Clear(); + } + p.Close(); + delete c1; +} + +void SetupStyle() +{ + const Int_t NCont=255; + + TStyle *st = new TStyle("mystyle","mystyle"); + gROOT->GetStyle("Plain")->Copy((*st)); + st->SetTitleX(0.1); + st->SetTitleW(0.8); + st->SetTitleH(0.08); + st->SetStatX(.9); + st->SetStatY(.9); + st->SetNumberContours(NCont); + st->SetPalette(1,0); + st->SetOptStat("erm"); + st->SetOptFit(0); + st->SetGridColor(kGray+1); + st->SetPadGridX(kTRUE); + st->SetPadGridY(kTRUE); + st->SetPadTickX(kTRUE); + st->SetPadTickY(kTRUE); + st->cd(); + + const Int_t NRGBs = 5; + Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; + Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 }; + Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 }; + Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 }; + + TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); + +} + +void DrawUnbinned(){ + TFile f("jpsi_debug.root"); +// if (!f.IsOpen()) return; + + TTree *t=(TTree*)f.Get("Pair"); +// if (!t) return; + + TCanvas c1; + gPad->SetLogy(); + gStyle->SetOptStat(0); + + TLegend *leg=new TLegend(0.59,.79,.99,.99); + TLine l; + + l.SetLineColor(kGreen-5); + l.SetLineWidth(2); + l.SetLineStyle(2); + leg->SetFillColor(10); + + leg->Clear(); + + + t->SetLineColor(kBlack); + t->Draw("M>>hAll(200,-.01,3.99)","","histe"); + TH1 *hAll=(TH1*)gROOT->FindObject("hAll"); + hAll->SetMinimum(0.1); + hAll->SetTitle(";M [GeV]; yield"); + leg->AddEntry(hAll,"|n#sigma e|<2 + pt>0.3 GeV","l"); + + l.DrawLine(3.097,1,3.097,1e4); + + t->SetLineColor(kOrange-5); + t->Draw("M>>hC11(200,-.01,3.99)","abs(Leg1_ImpactParXY)<.004&&abs(Leg2_ImpactParXY)<.004","histesame"); + hAll=(TH1*)gROOT->FindObject("hC11"); + leg->AddEntry(hAll,"|n#sigma e|<2 + pt>0.3 GeV + |dXY|<40#mum","l"); + + TCut d1_1="abs(Leg1_TPC_nSigma_Electrons)<1"; + TCut d2_1="abs(Leg2_TPC_nSigma_Electrons)<1"; + TCut d_1=d1_1+d2_1; + + t->SetLineColor(kRed); + t->Draw("M>>hC1(200,-.01,3.99)","Leg2_Pt>1&&Leg1_Pt>1","histesame"); + hAll=(TH1*)gROOT->FindObject("hC1"); + leg->AddEntry(hAll,"|n#sigma e|<2 + pt>1 GeV","l"); + + + t->SetLineColor(kGreen); + t->Draw("M>>hC2(200,-.01,3.99)",d_1+"Leg2_Pt>1&&Leg1_Pt>1","histesame"); + hAll=(TH1*)gROOT->FindObject("hC2"); + leg->AddEntry(hAll,"|n#sigma e|<1 + pt>1 GeV","l"); + + t->SetLineColor(kMagenta); + t->Draw("M>>hC3(200,-.01,3.99)","Leg1_Pt>2&&Leg2_Pt>2","histesame"); + hAll=(TH1*)gROOT->FindObject("hC3"); + leg->AddEntry(hAll,"|n#sigma e|<2 + pt>2 GeV","l"); + + t->SetLineColor(kBlue); + t->Draw("M>>hC4(200,-.01,3.99)","Leg1_Pt>3&&Leg2_Pt>3","histesame"); + hAll=(TH1*)gROOT->FindObject("hC4"); + leg->AddEntry(hAll,"|n#sigma e|<2 + pt>3 GeV","l"); + + leg->Draw(); +} + +/* + Double_t alephParameters[5]; + // simulation + alephParameters[0] = 2.15898e+00/50.; + alephParameters[1] = 1.75295e+01; + alephParameters[2] = 3.40030e-09; + alephParameters[3] = 1.96178e+00; + alephParameters[4] = 3.91720e+00; + Color_t color=kRed; + + TF1 *foProton = new TF1("foProton", "50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20); + TF1 *foPion = new TF1("foPion", "50*AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20); + TF1 *foElec = new TF1("foElec", "50*AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20); + TF1 *foKaon = new TF1("foKaon", "50*AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20); + TF1 *foMuon = new TF1("foMuon", "50*AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20); + // + foProton->SetParameters(alephParameters); + foPion->SetParameters(alephParameters); + foElec->SetParameters(alephParameters); + foKaon->SetParameters(alephParameters); + foMuon->SetParameters(alephParameters); + // + foProton->SetLineColor(color); + foPion->SetLineColor(color); + foElec->SetLineColor(color); + foKaon->SetLineColor(color); + foMuon->SetLineColor(color); + // + Int_t lineWidth=1; + foProton->SetLineWidth(lineWidth); + foPion->SetLineWidth(lineWidth); + foElec->SetLineWidth(lineWidth); + foKaon->SetLineWidth(lineWidth); + foMuon->SetLineWidth(lineWidth); + + // + foProton->SetNpx(200); + foPion->SetNpx(200); + foElec->SetNpx(200); + foKaon->SetNpx(200); + foMuon->SetNpx(200); + // + foProton->Draw("same"); + foPion->Draw("same"); + foElec->Draw("same"); + foKaon->Draw("same"); + foMuon->Draw("same"); + + + + + + // data + Double_t res=5.7e-2; + alephParameters[0] = 0.0283086; + alephParameters[1] = 2.63394e+01; + alephParameters[2] = 5.04114e-11; + alephParameters[3] = 2.12543e+00; + alephParameters[4] = 4.88663e+00; + Color_t color=kRed; + + + + TF1 *foDataProton = new TF1("foDataProton", "50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])",0.05,20); + TF1 *foDataProtonP = new TF1("foDataProtonP",Form( "50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])*(1+%f)",res),0.05,20); + TF1 *foDataProtonM = new TF1("foDataProtonM", Form("50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])*(1-%f)",res),0.05,20); + + TF1 *foDataPion = new TF1("foDataPion", "50*AliExternalTrackParam::BetheBlochAleph(x/0.13957,[0],[1],[2],[3],[4])",0.05,20); + TF1 *foDataPionP = new TF1("foDataPionP",Form( "50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])*(1+%f)",res),0.05,20); + TF1 *foDataPionM = new TF1("foDataPionM", Form("50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])*(1-%f)",res),0.05,20); + + TF1 *foDataElec = new TF1("foDataElec", "50*AliExternalTrackParam::BetheBlochAleph(x/0.000511,[0],[1],[2],[3],[4])",0.05,20); + TF1 *foDataElecP = new TF1("foDataElecP",Form( "50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])*(1+%f)",res),0.05,20); + TF1 *foDataElecM = new TF1("foDataElecM", Form("50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])*(1-%f)",res),0.05,20); + + TF1 *foDataKaon = new TF1("foDataKaon", "50*AliExternalTrackParam::BetheBlochAleph(x/0.493677,[0],[1],[2],[3],[4])",0.05,20); + TF1 *foDataKaonP = new TF1("foDataKaonP",Form( "50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])*(1+%f)",res),0.05,20); + TF1 *foDataKaonM = new TF1("foDataKaonM", Form("50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])*(1-%f)",res),0.05,20); + + TF1 *foDataMuon = new TF1("foDataMuon", "50*AliExternalTrackParam::BetheBlochAleph(x/0.105658,[0],[1],[2],[3],[4])",0.05,20); + TF1 *foDataMuonP = new TF1("foDataMuonP",Form( "50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])*(1+%f)",res),0.05,20); + TF1 *foDataMuonM = new TF1("foDataMuonM", Form("50*AliExternalTrackParam::BetheBlochAleph(x/0.93827,[0],[1],[2],[3],[4])*(1-%f)",res),0.05,20); + + // + foDataProton->SetParameters(alephParameters); + foDataProtonP->SetParameters(alephParameters); + foDataProtonM->SetParameters(alephParameters); + foDataPion->SetParameters(alephParameters); + foDataPionP->SetParameters(alephParameters); + foDataPionM->SetParameters(alephParameters); + foDataElec->SetParameters(alephParameters); + foDataElecP->SetParameters(alephParameters); + foDataElecM->SetParameters(alephParameters); + foDataKaon->SetParameters(alephParameters); + foDataKaonP->SetParameters(alephParameters); + foDataKaonM->SetParameters(alephParameters); + foDataMuon->SetParameters(alephParameters); + foDataMuonP->SetParameters(alephParameters); + foDataMuonM->SetParameters(alephParameters); + // + foDataProton->SetLineColor(color); + foDataProtonP->SetLineColor(color-4); + foDataProtonM->SetLineColor(color-4); + foDataPion->SetLineColor(color); + foDataPionP->SetLineColor(color-4); + foDataPionM->SetLineColor(color-4); + foDataElec->SetLineColor(color); + foDataElecP->SetLineColor(color-4); + foDataElecM->SetLineColor(color-4); + foDataKaon->SetLineColor(color); + foDataKaonP->SetLineColor(color-4); + foDataKaonM->SetLineColor(color-4); + foDataMuon->SetLineColor(color); + foDataMuonP->SetLineColor(color-4); + foDataMuonM->SetLineColor(color-4); + // + Int_t lineWidth=1; + foDataProton->SetLineWidth(lineWidth); + foDataProtonP->SetLineWidth(lineWidth); + foDataProtonM->SetLineWidth(lineWidth); + foDataPion->SetLineWidth(lineWidth); + foDataPionP->SetLineWidth(lineWidth); + foDataPionM->SetLineWidth(lineWidth); + foDataElec->SetLineWidth(lineWidth); + foDataElecP->SetLineWidth(lineWidth); + foDataElecM->SetLineWidth(lineWidth); + foDataKaon->SetLineWidth(lineWidth); + foDataKaonP->SetLineWidth(lineWidth); + foDataKaonM->SetLineWidth(lineWidth); + foDataMuon->SetLineWidth(lineWidth); + foDataMuonP->SetLineWidth(lineWidth); + foDataMuonM->SetLineWidth(lineWidth); + + // + foDataProtonP->SetLineStyle(2); + foDataProtonM->SetLineStyle(2); + foDataPionP->SetLineStyle(2); + foDataPionM->SetLineStyle(2); + foDataElecP->SetLineStyle(2); + foDataElecM->SetLineStyle(2); + foDataKaonP->SetLineStyle(2); + foDataKaonM->SetLineStyle(2); + foDataMuonP->SetLineStyle(2); + foDataMuonM->SetLineStyle(2); + + // + foDataProton->SetNpx(200); + foDataProtonP->SetNpx(200); + foDataProtonM->SetNpx(200); + foDataPion->SetNpx(200); + foDataPionP->SetNpx(200); + foDataPionM->SetNpx(200); + foDataElec->SetNpx(200); + foDataKaon->SetNpx(200); + foDataMuon->SetNpx(200); + // + foDataProton->Draw("same"); + foDataProtonP->Draw("same"); + foDataProtonM->Draw("same"); + foDataPion->Draw("same"); + foDataElec->Draw("same"); + foDataKaon->Draw("same"); + foDataMuon->Draw("same"); + + + + + +{ + + Int_t baseColors[5]={kRed, kGreen+1, kAzure-4, kMagenta, kCyan+1}; + Int_t sigmaColorOffset=1; + +Int_t baseColors[5]={kRed, kGreen+1, kAzure-4, kMagenta, kCyan+1}; + Int_t sigmaColorOffset=0; + + Double_t sigmas[5]={2,2,2,2,2}; + Double_t masses[5]; + + for (Int_t i=0; iSetParameters(alephParameters); + fBetheP[i]->SetParameters(alephParameters); + fBetheM[i]->SetParameters(alephParameters); + + fBethe[i]->SetLineColor(baseColors[i]); + fBetheP[i]->SetLineColor(baseColors[i]-sigmaColorOffset); + fBetheM[i]->SetLineColor(baseColors[i]-sigmaColorOffset); + + fBethe[i]->SetLineWidth(lineWidth); + fBetheP[i]->SetLineWidth(lineWidth); + fBetheM[i]->SetLineWidth(lineWidth); + + fBetheP[i]->SetLineStyle(2); + fBetheM[i]->SetLineStyle(2); + + fBethe[i]->SetNpx(200); + fBetheP[i]->SetNpx(200); + fBetheM[i]->SetNpx(200); +} + +for (Int_t i=0; i<5; ++i){ + fBethe[i]->Draw("same"); + fBetheP[i]->Draw("same"); + fBetheM[i]->Draw("same"); +} +} + +*/ + +/* + +c->SetLineColor(kRed); c->Draw("M","Leg1_Pt>1.2","same"); +c->SetLineColor(kGreen); c->Draw("M","Leg1_Pt>1.5","same"); +c->SetLineColor(kBlue); c->Draw("M","Leg1_Pt>2","same"); +c->SetLineColor(kMagenta); c->Draw("M","Leg1_Pt>2.5","same"); +c->SetLineColor(kGreen-4); c->Draw("M","Leg1_Pt>3","same"); + +c->SetLineColor(kRed-2); c->Draw("M","Leg1_Pt>1.2&&Leg2_Pt>.5","same"); +c->SetLineColor(kRed-4); c->Draw("M","Leg1_Pt>1.2&&Leg2_Pt>.7","same"); +c->SetLineColor(kRed-6); c->Draw("M","Leg1_Pt>1.2&&Leg2_Pt>1","same"); +c->SetLineColor(kRed-8); c->Draw("M","Leg1_Pt>1.2&&Leg2_Pt>1.2","same"); + + +c=Pair + + +c->SetAlias("cutE","Leg1_TPC_nSigma_Electrons>-2&&(Leg2_TPC_nSigma_Electrons)>-2"); +c->SetAlias("cutE","abs(Leg1_TPC_nSigma_Electrons)<2&&abs(Leg2_TPC_nSigma_Electrons)<2"); + +c->SetAlias("cutPi","Leg1_TPC_nSigma_Pions>2&&Leg2_TPC_nSigma_Pions>2"); +c->SetAlias("cutP","abs(Leg1_TPC_nSigma_Protons)>2&&abs(Leg2_TPC_nSigma_Protons)>2"); +c->SetAlias("cutK","abs(Leg1_TPC_nSigma_Kaons)>2&&abs(Leg2_TPC_nSigma_Kaons)>2"); +c->SetAlias("cutMu","abs(Leg1_TPC_nSigma_Muons)>1.5&&abs(Leg2_TPC_nSigma_Muons)>1.5"); + +c->SetAlias("cutE","abs(Leg1_TPC_nSigma_Electrons)<3&&abs(Leg2_TPC_nSigma_Electrons)<3"); +c->SetAlias("cutPi","abs(Leg1_TPC_nSigma_Pions)>2&&abs(Leg2_TPC_nSigma_Pions)>2"); +c->SetAlias("cutP","abs(Leg1_TPC_nSigma_Protons+.5)>2&&abs(Leg2_TPC_nSigma_Protons+.5)>2"); +c->SetAlias("pid","cutPi&&cutE"); + + + +c->SetAlias("pid","cutPi&&cutP&&cutK&&cutE"); +c->SetAlias("pid","cutPi&&cutP&&cutK&&cutMu"); +c->Draw("M>>h(50,1.99,3.99)","pid&&PairType==1&&Leg2_Pt>1","e"); +hist=h +TObjArray arr; +arr.AddAt(hist,1); + +AliDielectronSignalFunc fun +fun.SetDefaults(1); +fun.Process(&arr); +fun.Draw("samestat"); + + +c->SetLineColor(kBlack); +c->Draw("M>>h(50,1.99,3.99)","cutPi&&PairType==1","e"); + +c->SetLineColor(kRed); +c->Draw("M>>h2(50,1.99,3.99)","cutPi&&cutK&&PairType==1","esame"); + +c->SetLineColor(kBlue); +c->Draw("M>>h3(50,1.99,3.99)","cutPi&&cutP&&PairType==1","esame"); + +c->SetLineColor(kGreen); +c->Draw("M>>h4(50,1.99,3.99)","cutPi&&cutMu&&PairType==1","esame"); + + + + + + + +c->SetLineColor(kBlack); +c->Draw("Leg1_Pt>>hPt(50,1.99,3.99)","cutPi&&PairType==1","goff"); +c->Draw("Leg2_Pt>>hPt+","cutPi&&PairType==1","e"); + +c->SetLineColor(kRed); +c->Draw("Leg1_Pt>>hPt2(50,1.99,3.99)","cutPi&&cutK&&PairType==1","goff"); +c->Draw("Leg2_Pt>>hPt2+","cutPi&&cutK&&PairType==1","esame"); + +c->SetLineColor(kBlue); +c->Draw("Leg1_Pt>>hPt3(50,1.99,3.99)","cutPi&&cutP&&PairType==1","goff"); +c->Draw("Leg2_Pt>>hPt3+","cutPi&&cutP&&PairType==1","esame"); + +c->SetLineColor(kGreen); +c->Draw("Leg1_Pt>>hPt4(50,1.99,3.99)","cutPi&&cutLeg1_Ptu&&PairType==1","goff"); +c->Draw("Leg1_Pt>>hPt4+","cutPi&&cutMu&&PairType==1","esame"); + + + + +c->Draw("M>>hM5(100,1.99,3.99)","Leg1_TPC_signal>65&&Leg2_TPC_signal>65&&Leg2_Pt>1.3&&PairType==1","e"); + + +c=Pair +c->Draw("M>>hM5(100,1.99,3.99)","Leg1_TPC_signal>65&&Leg2_TPC_signal>65&&Leg2_Pt>1.&&PairType==1","e") +hist=hM5 +AliDielectronSignalFunc fun +fun.SetDefaults(1); +fun.SetFitRange(2,4) +fun.Process(hM) +fun.Draw("samestat"); + +c->Draw("M>>hM5(50,1.99,3.99)","Leg1_TPC_signal>65&&Leg2_TPC_signal>65&&Leg2_Pt>1.1&&PairType==1","e") + + + + +c->SetAlias("cut","Leg1_TPC_signal>65&&Leg2_TPC_signal>65&&Leg2_Pt>1.&&PairType==1") +c->SetAlias("cut","Leg2_Pt>1&&PairType==1") +c->SetAlias("cut","Leg1_TPC_signal>65&&Leg2_TPC_signal>65&&cutP&&cutK&&PairType==1") +c->SetAlias("cut","cutP&&PairType==1") +c->SetAlias("cut","PairType==1&&abs(Leg1_TOF_nSigma_Protons)>2&&abs(Leg2_TOF_nSigma_Protons)>2") + +c->SetAlias("cut","abs(Leg2_TOF_nSigma_Protons)>3&&abs(Leg2_TOF_nSigma_Pions)>") +c->SetAlias("cut","abs(Leg2_TOF_nSigma_Protons)>7-5.5/3*Leg2_P_InnerParam") +// histos +AliDielectronHistos h("h","h"); +h.AddClass("TPCsignal"); +h.UserHistogram("TPCsignal","sigTPC","TPC signal;P [GeV];TPC signal [arb. Units]",400,.3,40,400,0.,200.,0,0,kTRUE,kFALSE) +h.GetHistogram("TPCsignal","sigTPC")->SetDirectory(gDirectory) + +h.UserHistogram("TPCsignal","nSigE","TPC n #sigma Electrons;P [GeV];TPC n #sigma Electrons",400,.3,40.,400,-4.,4.,0,0,kTRUE,kFALSE) +h.GetHistogram("TPCsignal","nSigE")->SetDirectory(gDirectory) +h.UserHistogram("TPCsignal","nSigMu","TPC n #sigma Muons;P [GeV];TPC n #sigma Muons",400,.3,40.,400,-4.,4.,0,0,kTRUE,kFALSE) +h.GetHistogram("TPCsignal","nSigMu")->SetDirectory(gDirectory) +h.UserHistogram("TPCsignal","nSigPi","TPC n #sigma Pions;P [GeV];TPC n #sigma Pions",400,.3,40.,400,-4.,4.,0,0,kTRUE,kFALSE) +h.GetHistogram("TPCsignal","nSigPi")->SetDirectory(gDirectory) +h.UserHistogram("TPCsignal","nSigK","TPC n #sigma Kaons;P [GeV];TPC n #sigma Kaons",400,.3,40,400,-4,4,0,0,kTRUE,kFALSE) +h.GetHistogram("TPCsignal","nSigK")->SetDirectory(gDirectory) +h.UserHistogram("TPCsignal","nSigP","TPC n #sigma Protons;P [GeV];TPC n #sigma Protons",400,.3,40.,400,-4,4.,0,0,kTRUE,kFALSE) +h.GetHistogram("TPCsignal","nSigP")->SetDirectory(gDirectory) + + +c->Draw("Leg1_TPC_signal:Leg1_P_InnerParam>>sigTPC","cut","colz") +c->Draw("Leg2_TPC_signal:Leg2_P_InnerParam>>+sigTPC","cut","colz") + + +c->Draw("Leg1_TPC_nSigma_Electrons:Leg1_P_InnerParam>>nSigE","cut","colz") +c->Draw("Leg2_TPC_nSigma_Electrons:Leg2_P_InnerParam>>+nSigE","cut","colz") + +c->Draw("Leg1_TPC_nSigma_Muos:Leg1_P_InnerParam>>nSigMu","cut","goff") +c->Draw("Leg2_TPC_nSigma_Muons:Leg1_P_InnerParam>>nSigMu","cut","goff") +c->Draw("Leg2_TPC_nSigma_Muons:Leg2_P_InnerParam>>+nSigMu","cut","colz") + +c->Draw("Leg1_TPC_nSigma_Pions:Leg1_P_InnerParam>>nSigPi","cut","goff") +c->Draw("Leg2_TPC_nSigma_Pions:Leg2_P_InnerParam>>+nSigPi","cut","colz") + +c->Draw("Leg1_TPC_nSigma_Kaons:Leg1_P_InnerParam>>nSigK","cut","goff") +c->Draw("Leg2_TPC_nSigma_Kaons:Leg2_P_InnerParam>>+nSigK","cut","colz") + +c->Draw("Leg1_TPC_nSigma_Protons+.5:Leg1_P_InnerParam>>nSigP","cut","goff") +c->Draw("Leg2_TPC_nSigma_Protons+.5:Leg2_P_InnerParam>>+nSigP","cut","colz") + + + + +c->Draw("Leg1_TOF_nSigma_Protons:Leg1_P_InnerParam>>nSigP","cut","goff") +c->Draw("Leg2_TOF_nSigma_Protons:Leg2_P_InnerParam>>+nSigP","cut","colz") + + +Pair->SetScanField(0) +Pair->Scan("EventInFile:File.GetString()","","colsize=1 col=3.d:100.s") + + +AliDielectronSignalFunc sig; +sig.SetDefaults(1); + +//WooJins cuts: +c->SetAlias("cutE","abs(Leg1_TPC_nSigma_Electrons)<3&&abs(Leg2_TPC_nSigma_Electrons)<3"); +c->SetAlias("cutPi","abs(Leg1_TPC_nSigma_Pions)>2&&abs(Leg2_TPC_nSigma_Pions)>2"); +// c->SetAlias("cutPi","Leg1_TPC_signal>65&&Leg2_TPC_signal>65"); +c->SetAlias("cutP","abs(Leg1_TPC_nSigma_Protons)>2&&abs(Leg2_TPC_nSigma_Protons)>2"); +c->SetAlias("cutK","abs(Leg1_TPC_nSigma_Kaons)>2&&abs(Leg2_TPC_nSigma_Kaons)>2"); +c->SetAlias("pid","cutE&&cutPi&&cutP&&cutK"); + + +c->SetAlias("cutAdd","PairType==1&&abs(Leg1_ImpactParXY)<.02&&abs(Leg2_ImpactParXY)<.02&&Leg2_Pt>1.") +c->Draw("M>>hM(100,2,4)","cutAdd&&pid","e"); +h.Rebin(); +h.Rebin(); +h.Rebin(); +sig.Process(hM); +sig.Draw("samestat"); + + + + +//test +c->SetAlias("cutE","abs(Leg1_TPC_nSigma_Electrons)<3&&abs(Leg2_TPC_nSigma_Electrons)<2"); +c->SetAlias("cutPi","abs(Leg1_TPC_nSigma_Pions)>1&&abs(Leg2_TPC_nSigma_Pions)>2"); +c->SetAlias("cutPi","Leg1_TPC_signal>67&&Leg2_TPC_signal>67"); +c->SetAlias("cutP","abs(Leg1_TPC_nSigma_Protons)>2&&abs(Leg2_TPC_nSigma_Protons)>2"); +c->SetAlias("cutK","abs(Leg1_TPC_nSigma_Kaons)>2&&abs(Leg2_TPC_nSigma_Kaons)>2"); + + +c->SetAlias("pid","cutE&&cutPi&&cutP&&cutK"); +c->SetAlias("cutAdd","PairType==1&&abs(Leg1_ImpactParXY)<.03&&abs(Leg2_ImpactParXY)<.03&&Leg2_Pt>1") +c->SetAlias("cut","cutAdd&&pid") + +c->Draw("M>>hM(50,1.99,3.99)","cut","e"); +sig.Process(hM); +sig.Draw("samestat"); + + +c->SetAlias("cut","PairType==1") +c->Draw("Leg1_TPC_signal:Leg1_P_InnerParam>>sigTPC","cut","colz") +c->Draw("Leg2_TPC_signal:Leg2_P_InnerParam>>+sigTPC","cut","colz") + +//// +c->SetAlias("cutE","abs(Leg1_TPC_nSigma_Electrons)<2&&abs(Leg2_TPC_nSigma_Electrons)<2"); +c->SetAlias("cutPi","abs(Leg1_TPC_nSigma_Pions)>2&&abs(Leg2_TPC_nSigma_Pions)>2"); +c->SetAlias("cutPi2","Leg1_TPC_signal>65&&Leg2_TPC_signal>65"); +c->SetAlias("cutPi3","abs(Leg1_TPC_nSigma_Pions)>2.5&&abs(Leg2_TPC_nSigma_Pions)>2.5"); +c->SetAlias("cutP","abs(Leg1_TPC_nSigma_Protons)>2&&abs(Leg2_TPC_nSigma_Protons)>2"); +c->SetAlias("cutP2","abs(Leg1_TPC_nSigma_Protons)>1.5&&abs(Leg2_TPC_nSigma_Protons)>1.5"); +c->SetAlias("cutK","abs(Leg1_TPC_nSigma_Kaons)>2&&abs(Leg2_TPC_nSigma_Kaons)>2"); +c->SetAlias("cutdXY","abs(Leg1_ImpactParXY)<.02&&abs(Leg2_ImpactParXY)<.02") +c->SetAlias("cutPt","Leg2_Pt>1.") +c->SetAlias("cutPt2","Leg2_Pt>1.2") + +//---- +c->SetMarkerSize(0.7); + +c->SetMarkerStyle(20); +c->SetLineColor(kRed); +c->SetMarkerColor(kRed); +c->Draw("M>>hM(100,2,4)","cutPi&&cutPt","e"); + +c->SetMarkerStyle(21); +c->SetLineColor(kRed-1); +c->SetMarkerColor(kRed-2); +c->Draw("M>>hM2(100,2,4)","cutPi2&&cutPt","esame"); + +c->SetMarkerStyle(22); +c->SetLineColor(kRed-2); +c->SetMarkerColor(kRed-2); +c->Draw("M>>hM3(100,2,4)","cutPi3&&cutPt","esame"); + +//---- +c->SetMarkerStyle(20); +c->SetLineColor(kBlue); +c->SetMarkerColor(kBlue); +c->Draw("M>>hM4(100,2,4)","cutPi&&cutPt&&cutP","esame"); + +c->SetMarkerStyle(21); +c->SetLineColor(kBlue-1); +c->SetMarkerColor(kBlue-1); +c->Draw("M>>hM5(100,2,4)","cutPi2&&cutPt&&cutP","esame"); + +c->SetMarkerStyle(22); +c->SetLineColor(kBlue-2); +c->SetMarkerColor(kBlue-2); +c->Draw("M>>hM6(100,2,4)","cutPi3&&cutPt&&cutP","esame"); + +//---- + +c->SetMarkerStyle(20); +c->SetLineColor(kGreen); +c->SetMarkerColor(kGreen); +c->Draw("M>>hM7(100,2,4)","cutPi&&cutPt&&cutP2","esame"); + +c->SetMarkerStyle(21); +c->SetLineColor(kGreen-1); +c->SetMarkerColor(kGreen-1); +c->Draw("M>>hM8(100,2,4)","cutPi2&&cutPt&&cutP2","esame"); + +c->SetMarkerStyle(22); +c->SetLineColor(kGreen-2); +c->SetMarkerColor(kGreen-2); +c->Draw("M>>hM9(100,2,4)","cutPi3&&cutPt&&cutP2","esame"); + +//---- + + +c->SetMarkerStyle(20); +c->SetLineColor(kMagentha); +c->SetMarkerColor(kMagentha); +c->Draw("M>>hM7(100,2,4)","cutPi&&cutPt&&cutP2","esame"); + +c->SetMarkerStyle(21); +c->SetLineColor(kMagentha-1); +c->SetMarkerColor(kMagentha-1); +c->Draw("M>>hM8(100,2,4)","cutPi2&&cutPt&&cutP2","esame"); + +c->SetMarkerStyle(22); +c->SetLineColor(kMagentha-2); +c->SetMarkerColor(kMagentha-2); +c->Draw("M>>hM9(100,2,4)","cutPi3&&cutPt&&cutP2","esame"); + + +c->SetLineColor(kBlack); +c->SetMarkerColor(kBlue); +c->Draw("M>>hM4(100,2,4)","cutE&&cutPi&&cutK&&cutP&&cutdXY&&cutPt","esame"); + + + + +// +c->SetAlias("cutE","Leg1_TPC_nSigma_Electrons>-1.5&&Leg2_TPC_nSigma_Electrons>-1.5"); +// c->SetAlias("cutE","Leg1_TPC_signal>60&&Leg2_TPC_signal>60"); +c->SetAlias("cutP","abs(Leg1_TPC_nSigma_Protons)>3&&abs(Leg2_TPC_nSigma_Protons)>3") + +c->SetAlias("cutAdd","PairType==1&&abs(Leg1_ImpactParXY)<.03&&abs(Leg2_ImpactParXY)<.03&&Leg2_Pt>0") +c->SetAlias("cut","Leg2_Pt>1&&cutE&&cutP") + +c->Draw("M>>hM(50,1.99,3.99)","cut","e"); + +c->SetAlias("cutAdd","PairType==1&&abs(Leg1_ImpactParXY)<.03&&abs(Leg2_ImpactParXY)<.03&&Leg2_Pt>.8") +c->Draw("M>>hM2(50,1.99,3.99)","cut","esame"); + +c->SetAlias("cutAdd","PairType==1&&abs(Leg1_ImpactParXY)<.03&&abs(Leg2_ImpactParXY)<.03&&Leg2_Pt>1") +c->Draw("M>>hM3(50,1.99,3.99)","cut","esame"); + +c->SetAlias("cutAdd","PairType==1&&abs(Leg1_ImpactParXY)<.03&&abs(Leg2_ImpactParXY)<.03&&Leg2_Pt>1.2") +c->Draw("M>>hM4(50,1.99,3.99)","cut","esame"); + + +c->Draw("Leg1_TPC_signal:Leg1_P_InnerParam>>sigTPC","cut","goff") +c->Draw("Leg2_TPC_signal:Leg2_P_InnerParam>>+sigTPC","cut","goff") +c1->Modified();c1->Update() + + +c->SetAlias("cutE","Leg1_TPC_nSigma_Electrons>-1.5&&Leg2_TPC_nSigma_Electrons>-1.5"); +c->SetAlias("cut","Leg2_P_InnerParam>1.5&&cutE") +c->SetMarkerStyle(21); +c->Draw("M>>hM(50,1.99,3.99)","cut","e"); + +c->SetAlias("cutE","Leg1_TPC_nSigma_Electrons>-1.5+.8&&Leg2_TPC_nSigma_Electrons>-1.5+.8"); +c->SetAlias("cut","Leg2_P_InnerParam>1.5&&cutE") +c->SetMarkerStyle(20); +c->SetMarkerColor(kRed); +c->Draw("M>>hM2(50,1.99,3.99)","cut","esame"); + +*/ + diff --git a/PWG3/dielectron/macros/MakeEfficiencyReport.C b/PWG3/dielectron/macros/MakeEfficiencyReport.C new file mode 100644 index 00000000000..3b381c670eb --- /dev/null +++ b/PWG3/dielectron/macros/MakeEfficiencyReport.C @@ -0,0 +1,135 @@ +void SetupStyle(); + +void MakeEfficiencyReport(const char* outputFile="JpsiEffReport.pdf", + const char* histos="jpsi_HistosSE.root", + const char* cf="jpsi_CF.root") +{ + // + // Make a pdf file with the efficiency report + // + + SetupStyle(); + + AliDielectronCFdraw d(cf); + d.SetRangeUser("PairType",1,1); + d.SetRangeUser("Y",-.89,.9,"0"); + + + TFile f("jpsi_HistosSE.root"); + + AliDielectronHistos h; + TIter nextHists((TList*)f.Get("Dielectron_Histos")); + + TPaveText pt(.02,.6,.98,.8); + TText *t1=pt.AddText(""); + TText *t2=pt.AddText(""); + + TCanvas *c1=new TCanvas; + + TPDF p(outputFile); + + // + // Efficiency plots + // + t1->SetTitle("Efficiency plots"); + t2->SetTitle("Pair"); + pt.Draw(); + c1->Update(); + +// d.DrawEfficiency("Pt","1;2;4;10;16;8",0); + d.DrawEfficiency("Pt","4;16;22",0); + c1->Update(); + + + // + // Efficiency as a function of the impact parameter + // +// c1->Clear(); + Int_t refStep=16; + const Int_t impSteps=5; + Double_t impMax[impSteps]={.05,.04,.03,.02,.01}; + TString title=d.GetCFContainer()->GetStepTitle(refStep); +// d.DrawEfficiency("Pt:Y","1;2;4;10;16;8",0); + d.DrawEfficiency("Pt",Form("%d",refStep),0,"sameleg2"); +// c1->Update(); + + for (Int_t i=0; iSetStepTitle(refStep, (title+Form(" |leg dXY|<%.2f",impMax[i])).Data()); + d.DrawEfficiency("Pt",Form("%d",refStep),0,"same+leg2"); + } + + d.UnsetRangeUser("Leg1_ImpactParXY"); + d.UnsetRangeUser("Leg2_ImpactParXY"); + d.GetCFContainer()->SetStepTitle(refStep, title.Data()); + ((TH1*)gPad->GetListOfPrimitives()->FindObject("eff_16/00"))->SetMaximum(1.); + + c1->Update(); + + + c1->Clear(); + d.SetRangeUser("Y",-.89,.9); +// d.DrawEfficiency("Pt:Y","1;2;4;10;16;8",0,"colz2"); + d.DrawEfficiency("Pt:Y","4;16;22",0,"colz2"); + c1->Update(); + c1->Clear(); + + // + // Inv mass plots + // +// t1->SetTitle("Invariant Mass plots"); +// t2->SetTitle(""); +// pt.Draw(); +// c1->Update(); + + + // + // Make QA info + // + + t1->SetTitle("QA summary plots for"); + THashList *list=0x0; + while ( (list=(THashList*)nextHists()) ){ + h.SetHistogramList(*list); + t2->SetTitle(list->GetName()); + pt.Draw(); + c1->Update(); + h.Draw(); + c1->Clear(); + } + p.Close(); + delete c1; +} + +void SetupStyle() +{ + const Int_t NCont=255; + + TStyle *st = new TStyle("mystyle","mystyle"); + gROOT->GetStyle("Plain")->Copy((*st)); + st->SetTitleX(0.1); + st->SetTitleW(0.8); + st->SetTitleH(0.08); + st->SetStatX(.9); + st->SetStatY(.9); + st->SetNumberContours(NCont); + st->SetPalette(1,0); + st->SetOptStat("erm"); + st->SetOptFit(0); + st->SetGridColor(kGray+1); + st->SetPadGridX(kTRUE); + st->SetPadGridY(kTRUE); + st->SetPadTickX(kTRUE); + st->SetPadTickY(kTRUE); + st->cd(); + + const Int_t NRGBs = 5; + Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 }; + Double_t red[NRGBs] = { 0.00, 0.00, 0.87, 1.00, 0.51 }; + Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 }; + Double_t blue[NRGBs] = { 0.51, 1.00, 0.12, 0.00, 0.00 }; + + TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont); + +} \ No newline at end of file -- 2.43.5