#include <AliAODHandler.h>
#include <AliAnalysisManager.h>
#include <AliVEvent.h>
+#include <AliInputEventHandler.h>
+#include <AliESDInputHandler.h>
#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
//_________________________________________________________________________________
AliAnalysisTaskDielectronFilter::AliAnalysisTaskDielectronFilter(const char *name) :
- AliAnalysisTaskSE(name),
- fDielectron(0)
+ AliAnalysisTaskSE(name),
+ fDielectron(0),
+ fSelectPhysics(kTRUE)
{
//
// Constructor
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");
//
// Main loop. Called for every event
//
-
+
if (!fDielectron) return;
+ AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
+ AliESDInputHandler *esdHandler=0x0;
+ if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(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;j<obj->GetEntriesFast();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;it<aod->GetNumberOfTracks();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<AliAODHandler*>
- ((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();
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; }
AliDielectron *fDielectron; // J/psi framework object
+ Bool_t fSelectPhysics; // Whether to use physics selection
+
AliAnalysisTaskDielectronFilter(const AliAnalysisTaskDielectronFilter &c);
AliAnalysisTaskDielectronFilter& operator= (const AliAnalysisTaskDielectronFilter &c);
fTrackFilter("TrackFilter"),
fPairFilter("PairFilter"),
fPdgMother(443),
+ fPdgLeg1(11),
+ fPdgLeg2(11),
fHistos(0x0),
fPairCandidates(new TObjArray(10)),
fCfManagerPair(0x0),
fTrackFilter("TrackFilter"),
fPairFilter("PairFilter"),
fPdgMother(443),
+ fPdgLeg1(11),
+ fPdgLeg2(11),
fHistos(0x0),
fPairCandidates(new TObjArray(10)),
fCfManagerPair(0x0),
// Initialise objects
//
if (fCfManagerPair) fCfManagerPair->InitialiseContainer(fPairFilter);
-
+ if (fDebugTree) fDebugTree->SetDielectron(this);
}
//________________________________________________________________
Int_t end=ntrack2;
if (arr1==arr2) end=itrack1;
for (Int_t itrack2=0; itrack2<end; ++itrack2){
- //create the pair TODO: change hardcoded PID of tracks
- candidate->SetTracks(static_cast<AliVTrack*>(fTracks[arr1].UncheckedAt(itrack1)), 11,
- static_cast<AliVTrack*>(fTracks[arr2].UncheckedAt(itrack2)), 11);
+ //create the pair
+ candidate->SetTracks(static_cast<AliVTrack*>(fTracks[arr1].UncheckedAt(itrack1)), fPdgLeg1,
+ static_cast<AliVTrack*>(fTracks[arr2].UncheckedAt(itrack2)), fPdgLeg2);
candidate->SetType(pairIndex);
candidate->SetLabel(AliDielectronMC::Instance()->GetLabelMotherWithPdg(candidate,fPdgMother));
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)?
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
AliDielectron(const AliDielectron &c);
AliDielectron &operator=(const AliDielectron &c);
- ClassDef(AliDielectron,2);
+ ClassDef(AliDielectron,3);
};
inline void AliDielectron::InitPairCandidateArrays()
//Pure MC truth
if (fStepForMCtruth){
- fCfContainer->SetStepTitle(step++,"MC truth");
+ fCfContainer->SetStepTitle(step++,"MC truth");
}
//before cuts (MC truth)
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
}
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));
}
}
// 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);
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;
}
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()){
// }
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 );
#include <AliESDEvent.h>
#include <AliVTrack.h>
+#include "AliDielectron.h"
+#include "AliDielectronMC.h"
#include "AliDielectronPair.h"
#include "AliDielectronDebugTree.h"
fFileName("jpsi_debug.root"),
fNVars(0),
fNVarsLeg(0),
- fStreamer(0x0)
+ fStreamer(0x0),
+ fDielectron(0x0)
{
//
// Default Constructor
fFileName("jpsi_debug.root"),
fNVars(0),
fNVarsLeg(0),
- fStreamer(0x0)
+ fStreamer(0x0),
+ fDielectron(0x0)
{
//
// Named Constructor
//Get File and event information
TObjString fileName;
Int_t eventInFile=-1;
+ Int_t runNumber=-1;
TTree *t=man->GetTree();
if (t) {
if (han){
AliESDEvent *ev=dynamic_cast<AliESDEvent*>(han->GetEvent());
eventInFile=ev->GetEventNumberInFile();
+ runNumber=ev->GetRunNumber();
}
if (!fStreamer) fStreamer=new TTreeSRedirector(fFileName.Data());
(*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];
class TTreeSRedirector;
class AliDielectronPair;
+class AliDielectron;
class AliDielectronDebugTree : public TNamed {
public:
void Fill(AliDielectronPair *pair);
+ void SetDielectron(AliDielectron * const dielectron) { fDielectron=dielectron; }
+
void DeleteStreamer();
void WriteTree();
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);
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();
else if (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()->IsA()==AliAODHandler::Class()) type=kAOD;
AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
-
+
fgInstance=new AliDielectronMC(type);
fgInstance->SetHasMC(mcHandler!=0x0);
//
// 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;
//
// connect stack object from the mc handler
//
+ fMCEvent=0x0;
AliMCEventHandler* mcHandler = dynamic_cast<AliMCEventHandler*> (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; }
}
//____________________________________________________________
-AliMCParticle* AliDielectronMC::GetMCTrack(AliESDtrack* _track)
+AliMCParticle* AliDielectronMC::GetMCTrack( const AliESDtrack* _track)
{
//
// return MC track
}
//____________________________________________________________
-TParticle* AliDielectronMC::GetMCTrackFromStack(AliESDtrack* _track)
+TParticle* AliDielectronMC::GetMCTrackFromStack(const AliESDtrack* _track)
{
//
// return MC track from stack
}
//____________________________________________________________
-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<AliMCParticle *>(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<AliMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
+ return mcmother;
+}
+
+//____________________________________________________________
+AliAODMCParticle* AliDielectronMC::GetMCTrackMother(const AliAODMCParticle* _particle){
+ //
+ // return MC track mother
+ //
+ AliAODMCParticle* mcmother = dynamic_cast<AliAODMCParticle *>(fMCEvent->GetTrack(_particle->GetMother()));
+ return mcmother;
+}
+
+//____________________________________________________________
+TParticle* AliDielectronMC::GetMCTrackMotherFromStack(const AliESDtrack* _track)
{
//
// return MC track mother from stack
}
//____________________________________________________________
-Int_t AliDielectronMC::GetMCPID(AliESDtrack* _track)
+Int_t AliDielectronMC::GetMCPID(const AliESDtrack* _track)
{
//
// return PDG code of the track from the MC truth info
}
//____________________________________________________________
-Int_t AliDielectronMC::GetMCPIDFromStack(AliESDtrack* _track)
+Int_t AliDielectronMC::GetMCPIDFromStack(const AliESDtrack* _track)
{
//
// return MC PDG code from stack
}
//____________________________________________________________
-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
}
//____________________________________________________________
-Int_t AliDielectronMC::GetMotherPDGFromStack(AliESDtrack* _track)
+Int_t AliDielectronMC::GetMotherPDGFromStack(const AliESDtrack* _track)
{
//
// return PDG code of the mother track from stack
}
//____________________________________________________________
-Int_t AliDielectronMC::GetMCProcess(AliESDtrack* _track)
+Int_t AliDielectronMC::GetMCProcess(const AliESDtrack* _track)
{
//
// return process number of the track
}
//____________________________________________________________
-Int_t AliDielectronMC::GetMCProcessFromStack(AliESDtrack* _track)
+Int_t AliDielectronMC::GetMCProcessFromStack(const AliESDtrack* _track)
{
//
// return process number of the 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
}
//____________________________________________________________
-Int_t AliDielectronMC::GetMCProcessMotherFromStack(AliESDtrack* _track)
+Int_t AliDielectronMC::GetMCProcessMotherFromStack(const AliESDtrack* _track)
{
//
// return process number of the mother of the track
//
if (!fMCEvent) return kFALSE;
+ if (!particle) return kFALSE;
if (particle->IsA()==AliMCParticle::Class()){
return IsMCMotherToEEesd(static_cast<const AliMCParticle*>(particle),pdgMother);
if ((ilast-ifirst)!=1) return kFALSE;
AliMCParticle *firstD=static_cast<AliMCParticle*>(GetMCTrackFromMCEvent(ifirst));
AliMCParticle *secondD=static_cast<AliMCParticle*>(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;
AliAODMCParticle *firstD=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(ifirst));
AliAODMCParticle *secondD=static_cast<AliAODMCParticle*>(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;
// 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<AliMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
// 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<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle1->GetLabel()));
AliAODMCParticle *mcPart2=static_cast<AliAODMCParticle*>(GetMCTrackFromMCEvent(particle2->GetLabel()));
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<AliMCParticle*>(GetMCTrackFromMCEvent(daughter1->GetLabel()));
+ AliMCParticle *mcDaughter2=static_cast<AliMCParticle*>(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;
+
+}
+
+
+
+
+
+
+
+
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();
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);
--- /dev/null
+/*************************************************************************
+* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+* *
+* Author: The ALICE Off-line Project. *
+* Contributors are mentioned in the code where appropriate. *
+* *
+* Permission to use, copy, modify and distribute this software and its *
+* documentation strictly for non-commercial purposes is hereby granted *
+* without fee, provided that the above copyright notice appears in all *
+* copies and that both the copyright notice and this permission notice *
+* appear in the supporting documentation. The authors make no claims *
+* about the suitability of this software for any purpose. It is *
+* provided "as is" without express or implied warranty. *
+**************************************************************************/
+
+///////////////////////////////////////////////////////////////////////////
+// Dielectron PID //
+// //
+// //
+/*
+Detailed description
+
+
+*/
+// //
+///////////////////////////////////////////////////////////////////////////
+
+#include <TMath.h>
+#include <TF1.h>
+
+#include <AliVParticle.h>
+#include <AliLog.h>
+#include <AliESDtrack.h>
+
+#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<kNmaxPID; ++icut){
+ fDetType[icut]=kTPC;
+ fPartType[icut]=AliPID::kPion;
+ fNsigmaLow[icut]=0;
+ fNsigmaUp[icut]=0;
+ fPmin[icut]=0;
+ fPmax[icut]=0;
+ fExclude[icut]=kFALSE;
+ fFunUpperCut[icut]=0x0;
+ fFunLowerCut[icut]=0x0;
+ }
+}
+
+//______________________________________________
+AliDielectronPID::AliDielectronPID(const char* name, const char* title) :
+ AliAnalysisCuts(name, title),
+ fNcuts(0),
+ fRequirePIDbit(kTRUE),
+ fESDpid(0x0)
+{
+ //
+ // Named Constructor
+ //
+ for (Int_t icut=0; icut<kNmaxPID; ++icut){
+ fDetType[icut]=kTPC;
+ fPartType[icut]=AliPID::kPion;
+ fNsigmaLow[icut]=0;
+ fNsigmaUp[icut]=0;
+ fPmin[icut]=0;
+ fPmax[icut]=0;
+ fExclude[icut]=kFALSE;
+ fFunUpperCut[icut]=0x0;
+ fFunLowerCut[icut]=0x0;
+ }
+}
+
+//______________________________________________
+AliDielectronPID::~AliDielectronPID()
+{
+ //
+ // Default Destructor
+ //
+}
+
+//______________________________________________
+void AliDielectronPID::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*/)
+{
+ //
+ // Add a pid nsigma cut
+ // use response of detector 'det' in the momentum range ('pMin') to ['pMax']
+ // use a sigma band betwee 'nSigmaLow' and 'nSigmaUp'
+ // if nSigmaUp==-99999. then nSigmaLow will be uesd as a symmetric band +- nSigmaLow
+ // specify whether to 'exclude' the given band
+ //
+
+ if (fNcuts==kNmaxPID){
+ AliError(Form("only %d pid cut ranges allowed",kNmaxPID));
+ }
+ if (TMath::Abs(nSigmaUp+99999.)<1e-20){
+ nSigmaUp=TMath::Abs(nSigmaLow);
+ nSigmaLow=-1.*nSigmaUp;
+ }
+ fDetType[fNcuts]=det;
+ fPartType[fNcuts]=type;
+ fNsigmaLow[fNcuts]=nSigmaLow;
+ fNsigmaUp[fNcuts]=nSigmaUp;
+ fPmin[fNcuts]=pMin;
+ fPmax[fNcuts]=pMax;
+ fExclude[fNcuts]=exclude;
+ ++fNcuts;
+
+}
+
+//______________________________________________
+void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, Double_t nSigmaLow, TF1 * const funUp,
+ Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/)
+{
+ //
+ // cut using a TF1 as upper cut
+ //
+ if (funUp==0x0){
+ AliError("A valid function is required for the upper cut. Not adding the cut!");
+ return;
+ }
+ fFunUpperCut[fNcuts]=funUp;
+ AddCut(det,type,nSigmaLow,0.,pMin,pMax,exclude);
+}
+
+//______________________________________________
+void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, Double_t nSigmaUp,
+ Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/)
+{
+ //
+ // cut using a TF1 as lower cut
+ //
+ if (funLow==0x0){
+ AliError("A valid function is required for the lower cut. Not adding the cut!");
+ return;
+ }
+ fFunLowerCut[fNcuts]=funLow;
+ AddCut(det,type,0.,nSigmaUp,pMin,pMax,exclude);
+}
+
+//______________________________________________
+void AliDielectronPID::AddCut(DetType det, AliPID::EParticleType type, TF1 * const funLow, TF1 * const funUp,
+ Double_t pMin/*=0*/, Double_t pMax/*=0*/, Bool_t exclude/*=kFALSE*/)
+{
+ //
+ // cut using a TF1 as lower and upper cut
+ //
+ if ( (funUp==0x0) || (funLow==0x0) ){
+ AliError("A valid function is required for upper and lower cut. Not adding the cut!");
+ return;
+ }
+ fFunUpperCut[fNcuts]=funUp;
+ fFunLowerCut[fNcuts]=funLow;
+ AddCut(det,type,0.,0.,pMin,pMax,exclude);
+}
+
+//______________________________________________
+Bool_t AliDielectronPID::IsSelected(TObject* track)
+{
+ //
+ // perform PID cuts
+ //
+
+ //loop over all cuts
+ AliVParticle *part=static_cast<AliVParticle*>(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; icut<fNcuts; ++icut){
+ Double_t pMin=fPmin[icut];
+ Double_t pMax=fPmax[icut];
+
+ // test momentum range. In case pMin==pMax use all momenta
+ if ( (TMath::Abs(pMin-pMax)>1e-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<AliESDtrack*>(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<AliAODTrack*>(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<AliESDtrack*>(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<AliAODTrack*>(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<AliESDtrack*>(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<AliAODTrack*>(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 0<p<inf
+ // - exclude mu,K,pi,p 0<p<2
+
+ AddCut(kTPC,AliPID::kElectron,2.);
+ AddCut(kTPC,AliPID::kMuon,-2.,2.,0.,2.,kTRUE);
+ AddCut(kTPC,AliPID::kPion,-2.,2.,0.,2.,kTRUE);
+ AddCut(kTPC,AliPID::kKaon,-2.,2.,0.,2.,kTRUE);
+ AddCut(kTPC,AliPID::kProton,-2.,2.,0.,2.,kTRUE);
+
+ } else if (def==2) {
+ // include 2sigma e TPC
+ // 3sigma bands TOF
+ // - exclude K,P
+ AddCut(kTPC,AliPID::kElectron,2.);
+ AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,0.,kTRUE);
+ AddCut(kTOF,AliPID::kProton,-3.,3.,0.,0.,kTRUE);
+
+ } else if (def==3) {
+ // include 2sigma e TPC
+ // 3sigma bands TOF
+ // - exclude K 0<p<1
+ // - exclude P 0<p<2
+ AddCut(kTPC,AliPID::kElectron,2);
+ AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
+ AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
+ AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
+
+ } else if (def==4) {
+ // include 2sigma e TPC
+ // 3sigma bands TOF
+ // - exclude K 0<p<1
+ // - exclude P 0<p<2
+ AddCut(kTPC,AliPID::kElectron,2.);
+ AddCut(kTOF,AliPID::kKaon,-3.,3.,0.,1.,kTRUE);
+ AddCut(kTOF,AliPID::kProton,-6.,6.,0.,1.,kTRUE);
+ AddCut(kTOF,AliPID::kProton,-3.,3.,1.,2.,kTRUE);
+ fRequirePIDbit=kFALSE;
+ } else if (def==5) {
+ AddCut(kTPC,AliPID::kElectron,-0.5,3);
+ AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
+ } else if (def==6) {
+ // lower cut TPC: parametrisation by HFE
+ // upper cut TPC: 3 sigma
+ // TOF ele band 3sigma 0<p<1.5GeV
+ TF1 *lowerCut=new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 100);
+ lowerCut->SetParameters(-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<p<1.5GeV
+ AddCut(kTPC,AliPID::kElectron,10.);
+ AddCut(kTOF,AliPID::kElectron,-3,3,0,1.5);
+ }
+}
+
--- /dev/null
+#ifndef ALIDIELECTRONPID_H
+#define ALIDIELECTRONPID_H
+
+/* Copyright(c) 1998-2009, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice */
+
+//#############################################################
+//# #
+//# Class AliDielectronPID #
+//# #
+//# Authors: #
+//# Anton Andronic, GSI / A.Andronic@gsi.de #
+//# Ionut C. Arsene, GSI / I.C.Arsene@gsi.de #
+//# Julian Book, Uni Ffm / Julian.Book@cern.ch #
+//# Frederick Kramer, Uni Ffm, / Frederick.Kramer@cern.ch #
+//# Magnus Mager, CERN / Magnus.Mager@cern.ch #
+//# WooJin J. Park, GSI / W.J.Park@gsi.de #
+//# Jens Wiechula, Uni HD / Jens.Wiechula@cern.ch #
+//# #
+//#############################################################
+
+#include <AliPID.h>
+#include <AliESDpid.h>
+#include <AliAODTrack.h>
+#include <AliAODPid.h>
+
+#include <AliAnalysisCuts.h>
+
+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
#include "AliDielectronPair.h"
#include "AliVTrack.h"
+#include "AliPID.h"
ClassImp(AliDielectronPair)
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;
}
}
+//______________________________________________
+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<AliVParticle*>(fRefD1.GetObject());
+ AliVParticle *d2 = dynamic_cast<AliVParticle*>(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));
+ }
+}
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; }
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; }
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<AliVParticle*>(fRefD1.GetObject()); }
AliVParticle* GetSecondDaughter() const { return dynamic_cast<AliVParticle*>(fRefD2.GetObject()); }
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;
}
class AliDielectronPairLegCuts : public AliAnalysisCuts {
public:
- enum CutType { kBothLegs=0, kAnyLeg };
+ enum CutType { kBothLegs=0, kAnyLeg, kMixLegs };
AliDielectronPairLegCuts();
AliDielectronPairLegCuts(const char* name, const char* title);
AliDielectronSignalBase::AliDielectronSignalBase() :
TNamed(),
- fValues(4),
- fErrors(4),
- fIntMin(2.99),
+ fValues(6),
+ fErrors(6),
+ fIntMin(3.01),
fIntMax(3.19)
{
//
//______________________________________________
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)
{
//
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;
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;}
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);}
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)
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;}
fFitMin(2.5),
fFitMax(4),
fParM(-1),
- fParMres(-1)
+ fParMres(-1),
+ fSignalHist(0x0)
{
//
// Default Constructor
fFitMin(2.5),
fFitMax(4),
fParM(-1),
- fParMres(-1)
+ fParMres(-1),
+ fSignalHist(0x0)
{
//
// Named Constructor
return;
}
+ //set current signal hist pointer
+ fSignalHist=hist;
+
//initial the fit function to its default parameters
if (fVInitParam.GetMatrixArray()) fSigBack->SetParameters(fVInitParam.GetMatrixArray());
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;
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");
}
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();
+
}
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);
#include <TMath.h>
#include <TVectorT.h>
#include <TH1.h>
+#include <TROOT.h>
+#include <TCanvas.h>
#include <AliCFContainer.h>
#include <AliCFGridSparse.h>
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),
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),
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;
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){
for (Int_t iMethod=0; iMethod<fSignalMethods.GetEntries(); ++iMethod){
AliDielectronSignalBase *signalMethod=(AliDielectronSignalBase*)fSignalMethods.At(iMethod);
signalMethod->Process(&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; ivar<fNvars; ++ivar) binsProc+=Form("_%02d",fCurrentBins[ivar]);
+ signalMethod->Draw("stat");
+ binsProc.Append(".png");
+ binsProc.Prepend("SpectrumVisualDebug");
+ c->Update();
+ c->Print(binsProc.Data());
+ }
const TVectorD &val=signalMethod->GetValues();
const TVectorD &err=signalMethod->GetErrors();
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();
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; }
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
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
"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",
"TPC_nSigma_Muons",
"TPC_nSigma_Kaons",
"TPC_nSigma_Protons",
+ "TOF_nSigma_Electrons",
"TOF_nSigma_Pions",
"TOF_nSigma_Muons",
"TOF_nSigma_Kaons",
"DecayLength",
"R",
"OpeningAngle",
+ "ThetaHE",
+ "PhiHE",
+ "ThetaCS",
+ "PhiCS",
"LegDistance",
"LegDistanceXY",
"Merr",
#include <AliExternalTrackParam.h>
#include <AliESDpid.h>
+#include <AliPID.h>
#include "AliDielectronPair.h"
+#include "AliDielectronMC.h"
class AliVEvent;
class AliDielectronVarManager : public TNamed {
public:
-
+
// Particle specific variables
enum ValueTypes {
kPx = 0, // px
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
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
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
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
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();
// 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
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();
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);
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();
// 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)
// 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)
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())
+
+
}
// 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)
alephParameters[2] = 3.40030e-09;
alephParameters[3] = 1.96178e+00;
alephParameters[4] = 3.91720e+00;
+ fgESDpid->GetTOFResponse().SetTimeResolution(80.);
// data
if (type==1){
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);
}
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)
//----------------------
//create data containers
//----------------------
-
- //find input container
- AliAnalysisDataContainer *cinput = mgr->GetCommonInputContainer();
TString containerName = mgr->GetCommonFileName();
containerName += ":PWG3_dielectron";
--- /dev/null
+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;
+}
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();
InitHistograms(die,cutDefinition);
// the last definition uses no cuts and only the QA histograms should be filled!
- if (cutDefinition<nDie-1) InitCF(die,cutDefinition);
+// if (cutDefinition<nDie-1) InitCF(die,cutDefinition);
return die;
}
return;
}
+ // pt + 3 sigma ele TPC
+ if (cutDefinition==0){
+ AliDielectronVarCuts *pt = new AliDielectronVarCuts("Pt>1+|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);
}
esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
- esdTrackCuts->SetMinNClustersTPC(75);
+ esdTrackCuts->SetMinNClustersTPC(110);
esdTrackCuts->SetMaxChi2PerClusterTPC(4);
return esdTrackCuts;
// Setup the track cuts
//
- //ESD quality cuts
- die->GetTrackFilter().AddCuts(SetupESDtrackCuts(cutDefinition));
//QA no CF
--- /dev/null
+
+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 (cutDefinition<arrNames->GetEntriesFast()){
+ 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);
+
+}
--- /dev/null
+
+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+60<dEdx<100";
+ AliDielectron *die =
+ new AliDielectron(Form("%s",name.Data()),
+ Form("Track cuts: %s",name.Data()));
+
+ // cut setup
+ SetupTrackCuts(die);
+ SetupPairCuts(die);
+
+ //
+ // QA histogram setup
+ //
+ InitHistograms(die);
+
+ return die;
+}
+
+//______________________________________________________________________________________
+void SetupTrackCuts(AliDielectron *die)
+{
+ //
+ // Setup the track cuts
+ //
+
+ //ESD quality cuts
+ die->GetTrackFilter().AddCuts(SetupESDtrackCuts());
+
+ //Pt cut
+ AliDielectronVarCuts *pt = new AliDielectronVarCuts("Pt>.5+60<dEdx<100","Pt>.5 && 60<dEdx<100");
+ pt->AddCut(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","2<M<4");
+ invMassCut->AddCut(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);
+}
--- /dev/null
+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; i<AliPID::kSPECIES; ++i) masses[i]=AliPID::ParticleMass(i);
+
+ Double_t res=7.e-2;
+ Double_t alephParameters[5];
+
+ alephParameters[0] = 0.0283086;
+ alephParameters[1] = 2.63394e+01;
+ alephParameters[2] = 5.04114e-11;
+ alephParameters[3] = 2.12543e+00;
+ alephParameters[4] = 4.88663e+00;
+ Double_t mip=49.2;
+
+ Color_t color=kRed;
+ Int_t lineWidth=2;
+
+TF1 *fBethe[5];
+TF1 *fBetheP[5];
+TF1 *fBetheM[5];
+
+for (Int_t i=0; i<5; ++i){
+ fBethe[i] = new TF1(Form("fBethe%d",i), Form("%f*AliExternalTrackParam::BetheBlochAleph(x/%f,[0],[1],[2],[3],[4])",mip,masses[i]),0.05,20);
+ fBetheP[i] = new TF1(Form("fBethe%d",i), Form("%f*AliExternalTrackParam::BetheBlochAleph(x/%f,[0],[1],[2],[3],[4])*(1+%f*%f)",mip,masses[i],res,sigmas[i]),0.05,20);
+ fBetheM[i] = new TF1(Form("fBethe%d",i), Form("%f*AliExternalTrackParam::BetheBlochAleph(x/%f,[0],[1],[2],[3],[4])*(1-%f*%f)",mip,masses[i],res,sigmas[i]),0.05,20);
+
+ fBethe[i]->SetParameters(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");
+
+*/
+
--- /dev/null
+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; i<impSteps;++i){
+ d.SetRangeUser("Leg1_ImpactParXY",-1*impMax[i],impMax[i]);
+ d.SetRangeUser("Leg2_ImpactParXY",-1*impMax[i],impMax[i]);
+ d.GetCFContainer()->SetStepTitle(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
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff