From eac83c86b22667cac09a840d5b9d7857b58b17e0 Mon Sep 17 00:00:00 2001 From: wiechula Date: Thu, 10 May 2012 06:03:44 +0000 Subject: [PATCH] o macros and classes for Raa of 2010 data (Ionut) --- .../RaaPbPb2010/AddTask_iarsene_dst.C | 1011 +++++++++++++ .../AliAnalysisTaskCorrelationTree.cxx | 711 +++++++++ .../AliAnalysisTaskCorrelationTree.h | 104 ++ .../AliCorrelationReducedEvent.cxx | 502 +++++++ .../RaaPbPb2010/AliCorrelationReducedEvent.h | 458 ++++++ .../macrosJPSI/RaaPbPb2010/DstCommonMacros.C | 1302 +++++++++++++++++ .../RaaPbPb2010/RunDstPbPbJpsiAnalysis.C | 928 ++++++++++++ 7 files changed, 5016 insertions(+) create mode 100644 PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AddTask_iarsene_dst.C create mode 100644 PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliAnalysisTaskCorrelationTree.cxx create mode 100644 PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliAnalysisTaskCorrelationTree.h create mode 100644 PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliCorrelationReducedEvent.cxx create mode 100644 PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliCorrelationReducedEvent.h create mode 100644 PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/DstCommonMacros.C create mode 100644 PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/RunDstPbPbJpsiAnalysis.C diff --git a/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AddTask_iarsene_dst.C b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AddTask_iarsene_dst.C new file mode 100644 index 00000000000..4d9d19f5aed --- /dev/null +++ b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AddTask_iarsene_dst.C @@ -0,0 +1,1011 @@ +//__________________________________________________________________________________________ +AliAnalysisTask *AddTask_iarsene_dst(){ + //get the current analysis manager + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error("AddTask_iarsene_dst", "No analysis manager found."); + return 0; + } + + //Do we have an MC handler? + Bool_t hasMC=(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()!=0x0); + + //Get the current train configuration + //TString trainConfig=gSystem->Getenv("CONFIG_FILE"); + + gROOT->LoadMacro("AliCorrelationReducedEvent.cxx+g"); + gROOT->LoadMacro("AliAnalysisTaskCorrelationTree.cxx+g"); + + //create task and add it to the manager + AliAnalysisTaskCorrelationTree *task=new AliAnalysisTaskCorrelationTree("DSTTreeMaker"); + //if(trainConfig.Contains("PbPb")) task->SetTriggerMask(AliVEvent::kMB+AliVEvent::kCentral+AliVEvent::kSemiCentral); + //if(trainConfig=="pp") task->SetRejectPileup(); + if(!hasMC) task->UsePhysicsSelection(kFALSE); + mgr->AddTask(task); + + + task->SetEventFilter(CreateEventFilter()); + task->SetTrackFilter(CreateGlobalTrackFilter()); + task->SetFlowTrackFilter(CreateFlowTrackFilter()); + task->SetK0sPionCuts(CreateK0sPionCuts()); + task->SetLambdaProtonCuts(CreateLambdaProtonCuts()); + task->SetLambdaPionCuts(CreateLambdaPionCuts()); + task->SetK0sCuts(CreateK0sCuts()); + task->SetK0sMassRange(0.44,0.55); + task->SetLambdaMassRange(1.090,1.14); + task->SetLambdaCuts(CreateLambdaCuts()); + task->SetV0Histograms(CreateV0Histograms()); + + task->AddDielectron(ConfigDielectron(0)); // J/psi -> e+e- + task->AddDielectron(ConfigDielectron(1)); // phi -> K+K- + + //create output container + AliAnalysisDataContainer *coutput1 = + mgr->CreateContainer("diele_defaultTree", + TChain::Class(), + AliAnalysisManager::kExchangeContainer, + "diele_default"); + + AliAnalysisDataContainer *cOutputHist1 = + mgr->CreateContainer("qaHistos", + TList::Class(), + AliAnalysisManager::kOutputContainer, + "dst_qaHistos.root"); + + AliAnalysisDataContainer *cOutputHist2 = + mgr->CreateContainer("dstTree", + TTree::Class(), + AliAnalysisManager::kOutputContainer, + "dstTree.root"); + + AliAnalysisDataContainer *cOutputHist3 = + mgr->CreateContainer("friendTree", + TTree::Class(), + AliAnalysisManager::kOutputContainer, + "dstTree_friend.root"); + cout << "output containers: " << endl + << cOutputHist1 << endl + << cOutputHist2 << endl + << cOutputHist3 << endl; + + mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task, 0, coutput1 ); + mgr->ConnectOutput(task, 1, cOutputHist1); + //mgr->ConnectOutput(task, 2, cOutputHist2); + //mgr->ConnectOutput(task, 3, cOutputHist3); + mgr->ConnectOutput(task, 2, cOutputHist3); + + return task; +} + + +//_______________________________________________________________________________________________ +AliAnalysisCuts* CreateEventFilter() { + // + // Event wise cuts + // + AliDielectronEventCuts *eventCuts=new AliDielectronEventCuts("eventCuts","Vertex Track && |vtxZ|<10 && ncontrib>0"); + //eventCuts->SetRequireVertex(); + //eventCuts->SetMinVtxContributors(1); + //eventCuts->SetVertexZ(-10.,10.); + return eventCuts; +} + + +//_______________________________________________________________________________________________ +AliAnalysisCuts* CreateGlobalTrackFilter() { + // + // Cuts for tracks to be written in the dst + // + AliDielectronCutGroup* cuts = new AliDielectronCutGroup("cuts","cuts",AliDielectronCutGroup::kCompAND); + // general ESD cuts --------------------------------- + AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts; + // basic track quality cuts (basicQ) + //esdTrackCuts->SetMaxDCAToVertexZ(10.0); + //esdTrackCuts->SetMaxDCAToVertexXY(3.0); + //esdTrackCuts->SetEtaRange( -0.9 , 0.9 ); + //esdTrackCuts->SetAcceptKinkDaughters(kFALSE); + // esdTrackCuts->SetRequireITSRefit(kTRUE); + //esdTrackCuts->SetRequireTPCRefit(kTRUE); + esdTrackCuts->SetPRange(.1,1e30); + //esdTrackCuts->SetMinNClustersTPC(60); + //esdTrackCuts->SetMaxChi2PerClusterTPC(4); + // esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); + // esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst); + cuts->AddCut(esdTrackCuts); + + return cuts; +} + + +//_______________________________________________________________________________________________ +AliAnalysisCuts* CreateFlowTrackFilter() { + // + // Cuts for tracks to be used for the event plane q-vector + // These cuts are applied aditionally to the global track cuts + // + AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts; + esdTrackCuts->SetPtRange(0.2,2.0); + esdTrackCuts->SetEtaRange(-0.8, 0.8); + esdTrackCuts->SetMinNClustersTPC(70); + return esdTrackCuts; +} + + +//_______________________________________________________________________________________________ +AliAnalysisCuts* CreateK0sPionCuts() { + // + // Cuts on the K0s pions (tracking cuts, pid cuts, ...) + // + AliESDtrackCuts *pionCuts = new AliESDtrackCuts; + pionCuts->SetPtRange(0.15,100.0); + return pionCuts; +} + + +//_______________________________________________________________________________________________ +AliAnalysisCuts* CreateLambdaPionCuts() { + // + // Cuts on the Lambda pions (tracking cuts, pid cuts, ...) + // + AliESDtrackCuts *pionCuts = new AliESDtrackCuts; + pionCuts->SetPtRange(0.15,100.0); + return pionCuts; +} + + +//_______________________________________________________________________________________________ +AliAnalysisCuts* CreateLambdaProtonCuts() { + // + // Cuts on the Lambda protons (tracking cuts, pid cuts, ...) + // + AliESDtrackCuts *protonCuts = new AliESDtrackCuts; + protonCuts->SetPtRange(0.15,100.0); + return protonCuts; +} + + +//______________________________________________________________________________________ +AliESDv0Cuts* CreateK0sCuts() { + // + // Cuts on the V0s with K0s hypothesis + // + + AliESDv0Cuts *cuts=new AliESDv0Cuts(); + //cuts->SetMinDcaPosToVertex(-1); + //cuts->SetMinDcaNegToVertex(-1); + //cuts->SetMaxChi2(10); + //cuts->SetMaxDcaV0Daughters(0.3); + //cuts->SetMinRadius(3.0); + //cuts->SetMaxRadius(90.0); + //cuts->SetMinCosinePointingAngle(0.9); + //cuts->SetRequireOnFlyStatus(kTRUE); + //cuts->SetMaxDcaV0ToVertex(0.5); + //cuts->SetPRange(0.,1.0e10); + //cuts->SetPtRange(0.,1.0e10); + return cuts; +} + + +//______________________________________________________________________________________ +AliESDv0Cuts* CreateLambdaCuts() { + // + // Cuts on the V0s with Lambda hypothesis + // + + AliESDv0Cuts *cuts=new AliESDv0Cuts(); + //cuts->SetMinDcaPosToVertex(-1); + //cuts->SetMinDcaNegToVertex(-1); + //cuts->SetMaxChi2(10); + //cuts->SetMaxDcaV0Daughters(0.3); + //cuts->SetMinRadius(3.0); + //cuts->SetMaxRadius(90.0); + //cuts->SetMinCosinePointingAngle(0.9); + //cuts->SetRequireOnFlyStatus(kTRUE); + //cuts->SetMaxDcaV0ToVertex(0.5); + //cuts->SetPRange(0.,1.0e10); + //cuts->SetPtRange(0.,1.0e10); + return cuts; +} + + +//______________________________________________________________________________________ +AliDielectronHistos* CreateV0Histograms() +{ + // + // Initialise the histograms + // + + //Setup histogram Manager + AliDielectronHistos *histos=new AliDielectronHistos("V0Histograms",""); + + //add histograms to Track classes --------------------------------------------------------------------------- + histos->SetReservedWords("V0Track"); + + //Track classes + histos->AddClass("V0Track_Pos"); histos->AddClass("V0Track_Neg"); + + // kinematic acceptance histograms + histos->UserHistogram("V0Track","Pt","Pt;Pt (GeV/c);#tracks",400,0,20.,AliDielectronVarManager::kPt); + histos->UserHistogram("V0Track","Eta_P","Eta_P;#eta;P (GeV/c)",40, -1.0, +1.0, 200,0,10., AliDielectronVarManager::kEta, AliDielectronVarManager::kP); + histos->UserHistogram("V0Track","Eta_TRDPhi","Eta Phi Map; Eta; Phi;#tracks", + 40,-1,1,200,-3.15,3.15,AliDielectronVarManager::kEta,AliDielectronVarManager::kTRDphi); + + // DCA diagnostic histograms + histos->UserHistogram("V0Track","dXY","dXY;dXY (cm);#tracks",500,-1.,1.,AliDielectronVarManager::kImpactParXY); + histos->UserHistogram("V0Track","dZ","dZ;dZ (cm);#tracks",600,-3.,3.,AliDielectronVarManager::kImpactParZ); + histos->UserHistogram("V0Track","dZ_dXY","dZ_dXY;dZ (cm);dXY (cm)",600,-3.,3.,500,-1.,1.,AliDielectronVarManager::kImpactParZ, AliDielectronVarManager::kImpactParXY); + histos->UserHistogram("V0Track","P_dXY","P_dXY;P (GeV/c);dXY (cm)",200,0.,10.,500,-1.,1.,AliDielectronVarManager::kP, AliDielectronVarManager::kImpactParXY); + histos->UserHistogram("V0Track","P_dZ","P_dZ;P (GeV/c);dZ (cm)",200,0.,10.,600,-3.,3.,AliDielectronVarManager::kP, AliDielectronVarManager::kImpactParZ); + + // ITS + histos->UserHistogram("V0Track","ITSsignal","ITS signal;ITS dE/dx",1000,0.0,1000.0, AliDielectronVarManager::kITSsignal); + histos->UserHistogram("V0Track","ITSsignal_P","ITS signal vs P;P (GeV/c); ITS dE/dx",300, 0.0, 3.0, 200,0.0,1000.0, AliDielectronVarManager::kP, AliDielectronVarManager::kITSsignal); + + // TPC + histos->UserHistogram("V0Track","TPCchi2","#chi^{2}/cluster;#chi^{2}/cluster;#tracks",200,0.0,5.0, AliDielectronVarManager::kTPCchi2Cl); + histos->UserHistogram("V0Track","Pt_TPCchi2","#chi^{2}/cluster vs pt;p_{T} (GeV/c);#chi^{2}/cluster",400,0.,20., 200,0.0,5.0, + AliDielectronVarManager::kPt,AliDielectronVarManager::kTPCchi2Cl); + histos->UserHistogram("V0Track","TPCnCls","Number of Clusters TPC;TPC number clusters;#tracks",160,-0.5,159.5, AliDielectronVarManager::kNclsTPC); + histos->UserHistogram("V0Track","Phi_TPCnCls","Number of Clusters TPC vs #phi;#phi (rad.);TPC number clusters",200,0,6.285, 160,-0.5,159.5, + AliDielectronVarManager::kPhi,AliDielectronVarManager::kNclsTPC); + histos->UserHistogram("V0Track","Pin_TPCnCls","Number of Clusters TPC vs TPC inner param P;p (GeV/c);TPC number clusters",400,0,20., 160,-0.5,159.5, + AliDielectronVarManager::kPIn, AliDielectronVarManager::kNclsTPC); + histos->UserHistogram("V0Track","TPCnCls_dEdx","Number of Clusters TPC vs dE/dx; TPC number clusters; dE/dx (arb.units)", 160,-0.5,159.5, 200,0,200., + AliDielectronVarManager::kNclsTPC, AliDielectronVarManager::kTPCsignal); + histos->UserHistogram("V0Track","dEdx_P","dEdx;P [GeV];TPC signal (arb units);#tracks", + 400,0.2,20.,200,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal); + histos->UserHistogram("V0Track","dEdx_Eta","dEdx;#eta;TPC signal (arb units);#tracks", + 110,-1.1,1.1,200,0.,200.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTPCsignal); + histos->UserHistogram("V0Track","TPCnSigmaEle_P","TPC number of sigmas Electrons;P [GeV];TPC number of sigmas Electrons;#tracks", + 400,0.2,20.,200,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle); + histos->UserHistogram("V0Track","TPCnSigmaEle_Eta","TPC number of sigmas Electrons;#eta;TPC number of sigmas Electrons;#tracks", + 110,-1.1,1.1,200,-5.,5.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTPCnSigmaEle); + histos->UserHistogram("V0Track","TPCnSigmaPio_P","TPC number of sigmas Pions;P [GeV];TPC number of sigmas Pions;#tracks", + 400,0.2,20.,200,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaPio); + histos->UserHistogram("V0Track","TPCnSigmaPio_Eta","TPC number of sigmas Pions;#eta;TPC number of sigmas Pions;#tracks", + 110,-1.1,1.1,200,-5.,5.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTPCnSigmaPio); + histos->UserHistogram("V0Track","TPCnSigmaKao_P","TPC number of sigmas Kaons;P [GeV];TPC number of sigmas Kaons;#tracks", + 400,0.2,20.,200,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaKao); + histos->UserHistogram("V0Track","TPCnSigmaKao_Eta","TPC number of sigmas Kaons;#eta;TPC number of sigmas Kaons;#tracks", + 110,-1.1,1.1,200,-5.,5.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTPCnSigmaKao); + histos->UserHistogram("V0Track","TPCnSigmaPro_P","TPC number of sigmas Protons;P [GeV];TPC number of sigmas Protons;#tracks", + 400,0.2,20.,200,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaPro); + histos->UserHistogram("V0Track","TPCnSigmaPro_Eta","TPC number of sigmas Protons;#eta;TPC number of sigmas Protons;#tracks", + 110,-1.1,1.1,200,-5.,5.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTPCnSigmaPro); + histos->UserHistogram("V0Track","POut","Track POut;P (GeV/c)",100,0.0,10.0, AliDielectronVarManager::kPOut); + + // TOF + histos->UserHistogram("V0Track","Pin_TOFsignal","TOF signal vs TPC inner param P; P [GeV]; TOF signal", + 120,0.0,6.,1200,0.,120000.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFsignal); + histos->UserHistogram("V0Track","Pin_TOFbeta","TOF #beta vs TPC inner param P; P [GeV]; TOF #beta", + 120,0.0,6.,120,0.,1.2,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFbeta); + histos->UserHistogram("V0Track","TOFnSigmaPio_Pin","TOF number of sigmas Pions;P [GeV];TOF number of sigmas Pions;#tracks", + 400,0.2,20.,200,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFnSigmaPio); + histos->UserHistogram("V0Track","TOFnSigmaPio_eta","TOF number of sigmas Pions vs #eta;#eta;TOF number of sigmas Pions;#tracks", + 100,-1.2,1.2,200,-10.,10.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTOFnSigmaPio); + histos->UserHistogram("V0Track","TOFnSigmaKao_Pin","TOF number of sigmas Kaons;P [GeV];TOF number of sigmas Kaons;#tracks", + 400,0.2,20.,200,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFnSigmaKao); + histos->UserHistogram("V0Track","TOFnSigmaKao_eta","TOF number of sigmas Kaons vs #eta;#eta;TOF number of sigmas Kaons;#tracks", + 100,-1.2,1.2,200,-10.,10.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTOFnSigmaKao); + histos->UserHistogram("V0Track","TOFnSigmaPro_Pin","TOF number of sigmas Protons;P [GeV];TOF number of sigmas Protons;#tracks", + 400,0.2,20.,200,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFnSigmaPro); + histos->UserHistogram("V0Track","TOFnSigmaPro_eta","TOF number of sigmas Protons vs #eta;#eta;TOF number of sigmas Protons;#tracks", + 100,-1.2,1.2,200,-10.,10.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTOFnSigmaPro); + return histos; +} + + +//____________________________________________________________________________________________ +AliDielectron* ConfigDielectron(Int_t cutDefinition) +{ + // + // Setup the instance of an AliDielectron + // + const Char_t* dieleNames[2] = {"JpsiToEE", "PhiToKK"}; + + AliDielectron *die = + new AliDielectron(dieleNames[cutDefinition], + Form("Track cuts: %s", dieleNames)); + //die->SetEstimatorFilename("$TRAIN_ROOT/util/dielectron/dielectron/estimators.root"); + if(cutDefinition==0) { // J/psi -> e+e- + die->SetLegPdg(11,11); + die->SetMotherPdg(443); + } + if(cutDefinition==1) { // phi -> K+K- + die->SetLegPdg(321,321); + die->SetMotherPdg(333); + } + + // cut setup + SetupDielectronTrackCuts(die,cutDefinition); + SetupDielectronPairCuts(die,cutDefinition); + + // Set MC signals + SetDielectronMCSignals(die); + + // + // histogram setup + // only if an AliDielectronHistos object is attached to the + // dielectron framework histograms will be filled + // + InitDielectronHistograms(die,cutDefinition); + + //setup eta correction + SetEtaCorrection(); + + return die; +} + +//______________________________________________________________________________________ +void SetupDielectronTrackCuts(AliDielectron *die, Int_t cutDefinition) +{ + // + // Setup the track cuts + // + + AliDielectronCutGroup* cuts = new AliDielectronCutGroup("cuts","cuts",AliDielectronCutGroup::kCompAND); + die->GetTrackFilter().AddCuts(cuts); + //ESD quality cuts + AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts; + // basic track quality cuts (basicQ) + esdTrackCuts->SetMaxDCAToVertexZ(3.0); + esdTrackCuts->SetMaxDCAToVertexXY(1.0); + if(cutDefinition==1) { + esdTrackCuts->SetMaxDCAToVertexZ(0.3); + esdTrackCuts->SetMaxDCAToVertexXY(0.3); + } + if(cutDefinition==1) esdTrackCuts->SetEtaRange(-1.0,1.0); + if(cutDefinition==0) esdTrackCuts->SetEtaRange(-0.9,0.9); + esdTrackCuts->SetAcceptKinkDaughters(kFALSE); + //esdTrackCuts->SetRequireITSRefit(kTRUE); + if(cutDefinition==0) esdTrackCuts->SetRequireTPCRefit(kTRUE); + if(cutDefinition==0) esdTrackCuts->SetPRange(.9,1e30); // for J/psi + if(cutDefinition==1) esdTrackCuts->SetPRange(.1,1e30); // for phi + if(cutDefinition==0) esdTrackCuts->SetMinNClustersTPC(60); + if(cutDefinition==0) esdTrackCuts->SetMaxChi2PerClusterTPC(4); + //esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny); + cuts->AddCut(esdTrackCuts); + + Bool_t hasMC=(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()!=0x0); + + AliDielectronPID *pid = new AliDielectronPID("PID","PID cut"); + if(cutDefinition==0) { // J/psi->ee, phi->ee + pid->AddCut(AliDielectronPID::kTPC,AliPID::kElectron,-3.0, 4.5, 0.0, 0.0, kFALSE, AliDielectronPID::kRequire); // TPC 3-sigma inclusion for electron + pid->AddCut(AliDielectronPID::kTPC,AliPID::kProton, -100.0, 2.0, 0.0, 10.0, kTRUE, AliDielectronPID::kRequire); // TPC 3.0-sigma exclusion for proton +// pid->AddCut(AliDielectronPID::kTPC,AliPID::kKaon, -3.0, 3.0, 0.0, 10.0, kTRUE, AliDielectronPID::kRequire); // TPC 3.0-sigma exclusion for kaon +// pid->AddCut(AliDielectronPID::kTPC,AliPID::kPion, -3.0, 3.0, 0.0, 10.0, kTRUE, AliDielectronPID::kRequire); // TPC 3.0-sigma exclusion for pion + } + if(cutDefinition==1) { // phi -> KK + pid->AddCut(AliDielectronPID::kTOF,AliPID::kKaon, -3.5, 3.5, 0.6, 100.0, kFALSE, AliDielectronPID::kRequire, AliDielectronVarManager::kPIn); // TOF 3-sigma inclusion for kaon + pid->AddCut(AliDielectronPID::kTPC,AliPID::kKaon, -3.0, 3.5, 0.3, 100.0, kFALSE, AliDielectronPID::kRequire, AliDielectronVarManager::kPIn); // TPC 3-sigma inclusion for kaon + /* pid->AddCut(AliDielectronPID::kTPC,AliPID::kProton, -2.0, 2.0, 0.0, 0.0, kTRUE, AliDielectronPID::kRequire); // TPC 2.0-sigma exclusion for proton + pid->AddCut(AliDielectronPID::kTPC,AliPID::kPion, -2.0, 2.0, 0.0, 0.0, kTRUE, AliDielectronPID::kRequire); // TPC 2.0-sigma exclusion for pion + pid->AddCut(AliDielectronPID::kTOF,AliPID::kProton, -2.0, 2.0, 0.6, 3.0, kTRUE, AliDielectronPID::kIfAvailable); // TOF 2.0-sigma exclusion for proton + pid->AddCut(AliDielectronPID::kTOF,AliPID::kPion, -2.0, 2.0, 0.6, 3.0, kTRUE, AliDielectronPID::kIfAvailable); // TOF 2.0-sigma exclusion for pion +*/ + } + + // ++++++++++++++++++++++++++++++++++++++++ + // shifts for the nSigma electrons + TGraph* nSigmaCorrection = new TGraph(); + // b period (LHC10f7a) + nSigmaCorrection->SetPoint(0, 114000., -0.23-(0.0)); + nSigmaCorrection->SetPoint(1, 117222., -0.23-(0.0)); + // c period (LHC10f7a) + nSigmaCorrection->SetPoint(2, 119000., -0.23-(0.0)); + nSigmaCorrection->SetPoint(3, 120829., -0.23-(0.0)); + // d period (LHC10f7a) + nSigmaCorrection->SetPoint(4, 121000., -0.353-(-0.111)); + nSigmaCorrection->SetPoint(5, 127000., -0.353-(-0.111)); + // e period (LHC10f7a) + nSigmaCorrection->SetPoint(6, 127500., -0.351-(-0.116)); + nSigmaCorrection->SetPoint(7, 130900., -0.351-(-0.116)); + // LHC11a period (LHC11b10b) + nSigmaCorrection->SetPoint(8, 146680., -0.70-(+0.381)); // + nSigmaCorrection->SetPoint(9, 146870., -0.70-(+0.381)); // + if(hasMC) + pid->SetCorrGraph(nSigmaCorrection); + + cuts->AddCut(pid); +} + +//______________________________________________________________________________________ +void SetupDielectronPairCuts(AliDielectron *die, Int_t cutDefinition) +{ + // + // Setup the pair cuts + // + + //Invariant mass selection + AliDielectronVarCuts *cuts=new AliDielectronVarCuts("|y|<0.9","|Y|<.9"); + if(cutDefinition==0) cuts->AddCut(AliDielectronVarManager::kY,-0.9,0.9); + if(cutDefinition==1) cuts->AddCut(AliDielectronVarManager::kY,-1.0,1.0); + if(cutDefinition==1) cuts->AddCut(AliDielectronVarManager::kM, 0.95,1.1); + die->GetPairFilter().AddCuts(cuts); + + // Configure prefilter to remove conversions + AliDielectronVarCuts *gammaCut=new AliDielectronVarCuts("gammaCut","gammaCut"); + gammaCut->AddCut(AliDielectronVarManager::kM,0.,.05); + if(cutDefinition==0) { + die->GetPairPreFilter().AddCuts(gammaCut); + die->SetPreFilterUnlikeOnly(); + } +} + +//______________________________________________________________________________________ +void InitDielectronHistograms(AliDielectron *die, Int_t cutDefinition) +{ + // + // Initialise the histograms + // + + //Setup histogram Manager + AliDielectronHistos *histos=new AliDielectronHistos(die->GetName(),die->GetTitle()); + + //Initialise histogram classes + histos->SetReservedWords("Track;Pair"); + + //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))); + } + + //legs from pair + for (Int_t i=0; i<3; ++i){ + histos->AddClass(Form("Track_Legs_%s",AliDielectron::PairClassName(i))); + } + + //add histograms to event class + if (cutDefinition==0) { + histos->AddClass("Event"); + histos->UserHistogram("Event","VtxZ","Vertex Z;Z[cm]",500,-25.,25.,AliDielectronVarManager::kZvPrim); + histos->UserHistogram("Event","NTrk","NTrk",500,-0.5,499.5,AliDielectronVarManager::kNTrk); + histos->UserHistogram("Event","VtxZ_NaccTrcklts","VtxZ_NaccTrcklts;VtxZ;NaccTrcklts", + 240,-12.,12.,200,-0.5,199.5,AliDielectronVarManager::kZvPrim, AliDielectronVarManager::kNaccTrcklts); + + histos->UserHistogram("Event","MultV0A","VZERO mult",2000,0.,2000.,AliDielectronVarManager::kMultV0A); + histos->UserHistogram("Event","MultV0C","VZERO mult",2000,0.,2000.,AliDielectronVarManager::kMultV0C); + histos->UserHistogram("Event","MultV0","VZERO mult",2000,0.,2000.,AliDielectronVarManager::kMultV0); + histos->UserHistogram("Event","AdcV0A","VZERO adc",10000,0.,100000.,AliDielectronVarManager::kAdcV0A); + histos->UserHistogram("Event","AdcV0C","VZERO adc",10000,0.,100000.,AliDielectronVarManager::kAdcV0C); + histos->UserHistogram("Event","AdcV0","VZERO adc",10000,0.,100000.,AliDielectronVarManager::kAdcV0); + histos->UserHistogram("Event","MultV0A_MultV0C","MultV0A_MultV0C;MultV0A;MultV0C", + 2000,0.0,2000.,2000,0.0,2000,AliDielectronVarManager::kMultV0A, AliDielectronVarManager::kMultV0C); + histos->UserHistogram("Event","AdcV0A_AdcV0C","AdcV0A_AdcV0C;AdcV0A;AdcV0C", + 1000,0.0,20000.,1000,0.0,20000,AliDielectronVarManager::kAdcV0A, AliDielectronVarManager::kAdcV0C); + histos->UserHistogram("Event","MultV0A_AdcV0A","MultV0A_AdcV0A;MultV0A;AdcV0A", + 2000,0.0,2000.,2000,0.0,20000,AliDielectronVarManager::kMultV0A, AliDielectronVarManager::kAdcV0A); + histos->UserHistogram("Event","MultV0C_AdcV0C","MultV0C_AdcV0C;MultV0C;AdcV0C", + 2000,0.0,2000.,2000,0.0,20000,AliDielectronVarManager::kMultV0C, AliDielectronVarManager::kAdcV0C); + } + + //add histograms to Track classes --------------------------------------------------------------------------- + // kinematic acceptance histograms + histos->UserHistogram("Track","Pt","Pt;Pt (GeV/c);#tracks",400,0,20.,AliDielectronVarManager::kPt); + histos->UserHistogram("Track","Eta_P","Eta_P;#eta;P (GeV/c)",40, -1.0, +1.0, 200,0,10., AliDielectronVarManager::kEta, AliDielectronVarManager::kP); + histos->UserHistogram("Track","Eta_Phi","Eta Phi Map; Eta; Phi;#tracks", + 40,-1,1,200,0,6.285,AliDielectronVarManager::kEta,AliDielectronVarManager::kPhi); + + // DCA diagnostic histograms + histos->UserHistogram("Track","dXY","dXY;dXY (cm);#tracks",500,-1.,1.,AliDielectronVarManager::kImpactParXY); + histos->UserHistogram("Track","dZ","dZ;dZ (cm);#tracks",600,-3.,3.,AliDielectronVarManager::kImpactParZ); + histos->UserHistogram("Track","dZ_dXY","dZ_dXY;dZ (cm);dXY (cm)",600,-3.,3.,500,-1.,1.,AliDielectronVarManager::kImpactParZ, AliDielectronVarManager::kImpactParXY); + histos->UserHistogram("Track","P_dXY","P_dXY;P (GeV/c);dXY (cm)",200,0.,10.,500,-1.,1.,AliDielectronVarManager::kP, AliDielectronVarManager::kImpactParXY); + histos->UserHistogram("Track","P_dZ","P_dZ;P (GeV/c);dZ (cm)",200,0.,10.,600,-3.,3.,AliDielectronVarManager::kP, AliDielectronVarManager::kImpactParZ); + + // ITS + histos->UserHistogram("Track","ITSsignal","ITS signal;ITS dE/dx",1000,0.0,1000.0, AliDielectronVarManager::kITSsignal); + histos->UserHistogram("Track","ITSsignal_P","ITS signal vs P;P (GeV/c); ITS dE/dx",300, 0.0, 3.0, 400,0.0,1000.0, AliDielectronVarManager::kP, AliDielectronVarManager::kITSsignal); + + // TPC + histos->UserHistogram("Track","TPCchi2","#chi^{2}/cluster;#chi^{2}/cluster;#tracks",200,0.0,5.0, AliDielectronVarManager::kTPCchi2Cl); + histos->UserHistogram("Track","Pt_TPCchi2","#chi^{2}/cluster vs pt;p_{T} (GeV/c);#chi^{2}/cluster",400,0.,20., 200,0.0,5.0, + AliDielectronVarManager::kPt,AliDielectronVarManager::kTPCchi2Cl); + histos->UserHistogram("Track","TPCnCls","Number of Clusters TPC;TPC number clusters;#tracks",160,-0.5,159.5, AliDielectronVarManager::kNclsTPC); + histos->UserHistogram("Track","Phi_TPCnCls","Number of Clusters TPC vs #phi;#phi (rad.);TPC number clusters",200,0,6.285, 160,-0.5,159.5, + AliDielectronVarManager::kPhi,AliDielectronVarManager::kNclsTPC); + histos->UserHistogram("Track","Pin_TPCnCls","Number of Clusters TPC vs TPC inner param P;p (GeV/c);TPC number clusters",400,0,20., 160,-0.5,159.5, + AliDielectronVarManager::kPIn, AliDielectronVarManager::kNclsTPC); + histos->UserHistogram("Track","TPCnCls_dEdx","Number of Clusters TPC vs dE/dx; TPC number clusters; dE/dx (arb.units)", 160,-0.5,159.5, 200,0,200., + AliDielectronVarManager::kNclsTPC, AliDielectronVarManager::kTPCsignal); + histos->UserHistogram("Track","dEdx_P","dEdx;P [GeV];TPC signal (arb units);#tracks", + 400,0.2,20.,200,0.,200.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCsignal); + histos->UserHistogram("Track","TPCnSigmaEle_P","TPC number of sigmas Electrons;P [GeV];TPC number of sigmas Electrons;#tracks", + 400,0.2,20.,200,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaEle); + histos->UserHistogram("Track","TPCnSigmaPio_P","TPC number of sigmas Pions;P [GeV];TPC number of sigmas Pions;#tracks", + 400,0.2,20.,200,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaPio); + histos->UserHistogram("Track","TPCnSigmaKao_P","TPC number of sigmas Kaons;P [GeV];TPC number of sigmas Kaons;#tracks", + 400,0.2,20.,200,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaKao); + histos->UserHistogram("Track","TPCnSigmaPro_P","TPC number of sigmas Protons;P [GeV];TPC number of sigmas Protons;#tracks", + 400,0.2,20.,200,-5.,5.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTPCnSigmaPro); + histos->UserHistogram("Track","dEdx_Eta","dEdx;#eta;TPC signal (arb units);#tracks", + 100,-1.0,1.0,200,0.,200.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTPCsignal); + histos->UserHistogram("Track","TPCnSigmaEle_Eta","TPC number of sigmas Electrons;#eta;TPC number of sigmas Electrons;#tracks", + 100,-1.0,1.0,200,-5.,5.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTPCnSigmaEle); + histos->UserHistogram("Track","TPCnSigmaPio_Eta","TPC number of sigmas Pions;#eta;TPC number of sigmas Pions;#tracks", + 100,-1.0,1.0,200,-5.,5.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTPCnSigmaPio); + histos->UserHistogram("Track","TPCnSigmaKao_Eta","TPC number of sigmas Kaons;#eta;TPC number of sigmas Kaons;#tracks", + 100,-1.0,1.0,200,-5.,5.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTPCnSigmaKao); + histos->UserHistogram("Track","TPCnSigmaPro_Eta","TPC number of sigmas Protons;#eta;TPC number of sigmas Protons;#tracks", + 100,-1.0,1.0,200,-5.,5.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTPCnSigmaPro); + + // TOF + histos->UserHistogram("Track","Pin_TOFsignal","TOF signal vs TPC inner param P; P [GeV]; TOF signal", + 120,0.0,6.,1200,0.,120000.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFsignal); + histos->UserHistogram("Track","Pin_TOFbeta","TOF #beta vs TPC inner param P; P [GeV]; TOF #beta", + 120,0.0,6.,120,0.,1.2,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFbeta); + histos->UserHistogram("Track","TOFnSigmaPio_Pin","TOF number of sigmas Pions;P [GeV];TOF number of sigmas Pions;#tracks", + 400,0.2,20.,200,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFnSigmaPio); + histos->UserHistogram("Track","TOFnSigmaPio_eta","TOF number of sigmas Pions vs #eta;#eta;TOF number of sigmas Pions;#tracks", + 100,-1.2,1.2,200,-10.,10.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTOFnSigmaPio); + histos->UserHistogram("Track","TOFnSigmaKao_Pin","TOF number of sigmas Kaons;P [GeV];TOF number of sigmas Kaons;#tracks", + 400,0.2,20.,200,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFnSigmaKao); + histos->UserHistogram("Track","TOFnSigmaKao_eta","TOF number of sigmas Kaons vs #eta;#eta;TOF number of sigmas Kaons;#tracks", + 100,-1.2,1.2,200,-10.,10.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTOFnSigmaKao); + histos->UserHistogram("Track","TOFnSigmaPro_Pin","TOF number of sigmas Protons;P [GeV];TOF number of sigmas Protons;#tracks", + 400,0.2,20.,200,-10.,10.,AliDielectronVarManager::kPIn,AliDielectronVarManager::kTOFnSigmaPro); + histos->UserHistogram("Track","TOFnSigmaPro_eta","TOF number of sigmas Protons vs #eta;#eta;TOF number of sigmas Protons;#tracks", + 100,-1.2,1.2,200,-10.,10.,AliDielectronVarManager::kEta,AliDielectronVarManager::kTOFnSigmaPro); + + //add histograms to Pair classes + Int_t nMassBins = 10000; + Double_t massRange[2] = {-0.005, 100.005}; + if(cutDefinition==1) { // Phi -> K+K- + nMassBins = 4000; + massRange[0] = 0.8; massRange[1] = 1.2; + } + histos->UserHistogram("Pair","InvMass","Inv.Mass;Inv. Mass [GeV];#pairs", + nMassBins, massRange[0], massRange[1], AliDielectronVarManager::kM); + histos->UserHistogram("Pair","InvMass_Pt","Inv.Mass vs pt;Inv. Mass [GeV]; pt (GeV/c);#pairs", + nMassBins/10, massRange[0], massRange[1], 100, 0.0, 50.0, AliDielectronVarManager::kM, AliDielectronVarManager::kPt); + histos->UserHistogram("Pair","InvMass_Cent","Inv.Mass vs centrality;Inv. Mass [GeV]; centrality;#pairs", + nMassBins/10, massRange[0], massRange[1], 20, 0.0, 100.0, AliDielectronVarManager::kM, AliDielectronVarManager::kCentrality); + histos->UserHistogram("Pair","InvMass_Rap","Inv.Mass vs rapidity;Inv. Mass [GeV]; rapidity;#pairs", + nMassBins/10, massRange[0], massRange[1], 100, -1.0, 1.0, AliDielectronVarManager::kM, AliDielectronVarManager::kY); + histos->UserHistogram("Pair","InvMass_OpAngle","Inv.Mass vs opening angle;Inv. Mass [GeV]; op. angle (rad.);#pairs", + nMassBins/10, massRange[0], massRange[1], 100, 0.0, 4.0, AliDielectronVarManager::kM, AliDielectronVarManager::kOpeningAngle); + histos->UserHistogram("Pair","InvMass_ThetaHE","Inv.Mass vs cos #theta^{*}_{HE};Inv. Mass [GeV]; cos #theta^{*}_{HE};#pairs", + nMassBins/10, massRange[0], massRange[1], 100, -1.0, 1.0, AliDielectronVarManager::kM, AliDielectronVarManager::kThetaHE); + histos->UserHistogram("Pair","InvMass_ThetaCS","Inv.Mass vs cos #theta^{*}_{CS};Inv. Mass [GeV]; cos #theta^{*}_{CS};#pairs", + nMassBins/10, massRange[0], massRange[1], 100, -1.0, 1.0, AliDielectronVarManager::kM, AliDielectronVarManager::kThetaCS); + histos->UserHistogram("Pair","Pt","Transverse momentum; P_{t} [GeV/c];#pairs", + 1000, 0.0, 50.0, AliDielectronVarManager::kPt); + histos->UserHistogram("Pair","Pt_Y","Transverse momentum vs rapidity; rapidity; P_{t} [GeV/c];#pairs", + 100, -1.0, 1.0, 200, 0.0, 50.0, AliDielectronVarManager::kY, AliDielectronVarManager::kPt); + histos->UserHistogram("Pair","Pt_ThetaCS","Transverse momentum vs cos #theta^{*}_{CS}; cos #theta^{*}_{CS}; P_{t} [GeV/c];#pairs", + 100, -1.0, 1.0, 200, 0.0, 50.0, AliDielectronVarManager::kThetaCS, AliDielectronVarManager::kPt); + histos->UserHistogram("Pair","Pt_ThetaHE","Transverse momentum vs cos #theta^{*}_{HE}; cos #theta^{*}_{HE}; P_{t} [GeV/c];#pairs", + 100, -1.0, 1.0, 200, 0.0, 50.0, AliDielectronVarManager::kThetaHE, AliDielectronVarManager::kPt); + histos->UserHistogram("Pair","Pt_Phi","Transverse momentum vs phi; #phi (rad.); P_{t} [GeV/c];#pairs", + 200, -6.3, 6.3, 200, 0.0, 50.0, AliDielectronVarManager::kPhi, AliDielectronVarManager::kPt); + histos->UserHistogram("Pair","Phi","#phi; #phi (rad.);#pairs", + 1000, -6.3, 6.3, AliDielectronVarManager::kPhi); + histos->UserHistogram("Pair","Rapidity","Rapidity;Rapidity;#pairs", + 100,-1.,1.,AliDielectronVarManager::kY); + histos->UserHistogram("Pair","Rapidity_OpAngle","Rapidity vs opening angle;Rapidity; Op.Angle (rad.);#pairs", + 100,-1.,1.,200, 0.0, 4.0, AliDielectronVarManager::kY, AliDielectronVarManager::kOpeningAngle); + histos->UserHistogram("Pair","Rapidity_ThetaCS","Rapidity vs cos #theta^{*}_{CS};Rapidity; cos #theta^{*}_{CS};#pairs", + 100,-1.,1.,200, -1.0, 1.0, AliDielectronVarManager::kY, AliDielectronVarManager::kThetaCS); + histos->UserHistogram("Pair","Rapidity_ThetaHE","Rapidity vs cos #theta^{*}_{HE};Rapidity; cos #theta^{*}_{HE};#pairs", + 100,-1.,1.,200, -1.0, 1.0, AliDielectronVarManager::kY, AliDielectronVarManager::kThetaHE); + histos->UserHistogram("Pair","OpeningAngle","Opening angle;angle", + 100,0.,3.15,AliDielectronVarManager::kOpeningAngle); + histos->UserHistogram("Pair","Radius","Radius;R[cm]", + 1000,0.,300.,AliDielectronVarManager::kR); + histos->UserHistogram("Pair","ThetaHE","cos #theta^{*}_{HE}; cos #theta^{*}_{HE}", + 220,-1.1,1.1,AliDielectronVarManager::kThetaHE); + histos->UserHistogram("Pair","PhiHE","#varphi^{*}_{HE};#varphi^{*}_{HE} (rad.)", + 160,-3.2,3.2,AliDielectronVarManager::kPhiHE); + histos->UserHistogram("Pair","ThetaCS","cos #theta^{*}_{CS}; cos #theta^{*}_{CS}", + 220,-1.1,1.1,AliDielectronVarManager::kThetaCS); + histos->UserHistogram("Pair","PhiCS","#varphi^{*}_{CS};#varphi^{*}_{CS} (rad.)", + 160,-3.2,3.2,AliDielectronVarManager::kPhiCS); + histos->UserHistogram("Pair","Lxy","Pseudo-proper time;Lxy (cm.)", + 1000,0.0,2.0,AliDielectronVarManager::kPseudoProperTime); + histos->UserHistogram("Pair","Chi2NDF","Pair #chi^{2}/NDF; #chi^{2}/NDF", + 1000,0.0,10.0,AliDielectronVarManager::kChi2NDF); + + + die->SetHistogramManager(histos); +} + + +//_____________________________________________________________________________________________ +void SetDielectronMCSignals(AliDielectron *die) +{ + // J/psi sources (direct + feed down from charm higher states) + // 0 + AliDielectronSignalMC* promptJpsi = new AliDielectronSignalMC("promptJpsi","Prompt J/psi"); // prompt J/psi (not from beauty decays) + promptJpsi->SetLegPDGs(11,-11); + promptJpsi->SetMotherPDGs(443,443); + promptJpsi->SetGrandMotherPDGs(503,503,kTRUE,kTRUE); // not from beauty hadrons + promptJpsi->SetMothersRelation(AliDielectronSignalMC::kSame); + promptJpsi->SetFillPureMCStep(kTRUE); + promptJpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + promptJpsi->SetCheckBothChargesLegs(kTRUE,kTRUE); + promptJpsi->SetCheckBothChargesMothers(kTRUE,kTRUE); + promptJpsi->SetCheckBothChargesGrandMothers(kTRUE,kTRUE); + die->AddSignalMC(promptJpsi); + + // 1 + AliDielectronSignalMC* beautyJpsi = new AliDielectronSignalMC("beautyJpsi","beauty hadron -> J/psi"); // J/psi->e+e- from beauty hadron decays + beautyJpsi->SetLegPDGs(11,-11); + beautyJpsi->SetMotherPDGs(443,443); + beautyJpsi->SetMothersRelation(AliDielectronSignalMC::kSame); + beautyJpsi->SetGrandMotherPDGs(503,503); + beautyJpsi->SetFillPureMCStep(kTRUE); + beautyJpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + beautyJpsi->SetCheckBothChargesLegs(kTRUE,kTRUE); + beautyJpsi->SetCheckBothChargesMothers(kTRUE,kTRUE); + beautyJpsi->SetCheckBothChargesGrandMothers(kTRUE,kTRUE); + die->AddSignalMC(beautyJpsi); + + // 2 + AliDielectronSignalMC* beautyMesonJpsi = new AliDielectronSignalMC("beautyMesonJpsi","beauty meson -> J/psi"); // J/psi->e+e- from beauty meson decays + beautyMesonJpsi->SetLegPDGs(11,-11); + beautyMesonJpsi->SetMotherPDGs(443,443); + beautyMesonJpsi->SetMothersRelation(AliDielectronSignalMC::kSame); + beautyMesonJpsi->SetGrandMotherPDGs(500,500); + beautyMesonJpsi->SetFillPureMCStep(kTRUE); + beautyMesonJpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + beautyMesonJpsi->SetCheckBothChargesLegs(kTRUE,kTRUE); + beautyMesonJpsi->SetCheckBothChargesMothers(kTRUE,kTRUE); + beautyMesonJpsi->SetCheckBothChargesGrandMothers(kTRUE,kTRUE); + die->AddSignalMC(beautyMesonJpsi); + + // 3 + AliDielectronSignalMC* openBeautyMesonJpsi = new AliDielectronSignalMC("openBeautyMesonJpsi","open beauty meson -> J/psi"); // J/psi->e+e- from open beauty meson decays + openBeautyMesonJpsi->SetLegPDGs(11,-11); + openBeautyMesonJpsi->SetMotherPDGs(443,443); + openBeautyMesonJpsi->SetMothersRelation(AliDielectronSignalMC::kSame); + openBeautyMesonJpsi->SetGrandMotherPDGs(501,501); + openBeautyMesonJpsi->SetFillPureMCStep(kTRUE); + openBeautyMesonJpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + openBeautyMesonJpsi->SetCheckBothChargesLegs(kTRUE,kTRUE); + openBeautyMesonJpsi->SetCheckBothChargesMothers(kTRUE,kTRUE); + openBeautyMesonJpsi->SetCheckBothChargesGrandMothers(kTRUE,kTRUE); + die->AddSignalMC(openBeautyMesonJpsi); + + // 4 + AliDielectronSignalMC* chic0Jpsi = new AliDielectronSignalMC("chic0Jpsi","chic0 -> J/psi + X"); // J/psi->e+e- from chic0 decays + chic0Jpsi->SetLegPDGs(11,-11); + chic0Jpsi->SetMotherPDGs(443,443); + chic0Jpsi->SetMothersRelation(AliDielectronSignalMC::kSame); + chic0Jpsi->SetGrandMotherPDGs(10441,10441); + chic0Jpsi->SetFillPureMCStep(kTRUE); + chic0Jpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + chic0Jpsi->SetCheckBothChargesLegs(kTRUE,kTRUE); + chic0Jpsi->SetCheckBothChargesMothers(kTRUE,kTRUE); + chic0Jpsi->SetCheckBothChargesGrandMothers(kTRUE,kTRUE); + die->AddSignalMC(chic0Jpsi); + + // 5 + AliDielectronSignalMC* chic1Jpsi = new AliDielectronSignalMC("chic1Jpsi","chic1 -> J/psi + X"); // J/psi->e+e- from chic1 decays + chic1Jpsi->SetLegPDGs(11,-11); + chic1Jpsi->SetMotherPDGs(443,443); + chic1Jpsi->SetMothersRelation(AliDielectronSignalMC::kSame); + chic1Jpsi->SetGrandMotherPDGs(20443,20443); + chic1Jpsi->SetFillPureMCStep(kTRUE); + chic1Jpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + chic1Jpsi->SetCheckBothChargesLegs(kTRUE,kTRUE); + chic1Jpsi->SetCheckBothChargesMothers(kTRUE,kTRUE); + chic1Jpsi->SetCheckBothChargesGrandMothers(kTRUE,kTRUE); + die->AddSignalMC(chic1Jpsi); + + // 6 + AliDielectronSignalMC* chic2Jpsi = new AliDielectronSignalMC("chic2Jpsi","chic2 -> J/psi + X"); // J/psi->e+e- from chic2 decays + chic2Jpsi->SetLegPDGs(11,-11); + chic2Jpsi->SetMotherPDGs(443,443); + chic2Jpsi->SetMothersRelation(AliDielectronSignalMC::kSame); + chic2Jpsi->SetGrandMotherPDGs(445,445); + chic2Jpsi->SetFillPureMCStep(kTRUE); + chic2Jpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + chic2Jpsi->SetCheckBothChargesLegs(kTRUE,kTRUE); + chic2Jpsi->SetCheckBothChargesMothers(kTRUE,kTRUE); + chic2Jpsi->SetCheckBothChargesGrandMothers(kTRUE,kTRUE); + die->AddSignalMC(chic2Jpsi); + + // 7 + AliDielectronSignalMC* psiPrimeJpsi = new AliDielectronSignalMC("psiPrimeJpsi","psi(2S) -> J/psi + X"); // J/psi->e+e- from psi(2S) decays + psiPrimeJpsi->SetLegPDGs(11,-11); + psiPrimeJpsi->SetMotherPDGs(443,443); + psiPrimeJpsi->SetMothersRelation(AliDielectronSignalMC::kSame); + psiPrimeJpsi->SetGrandMotherPDGs(445,445); + psiPrimeJpsi->SetFillPureMCStep(kTRUE); + psiPrimeJpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + psiPrimeJpsi->SetCheckBothChargesLegs(kTRUE,kTRUE); + psiPrimeJpsi->SetCheckBothChargesMothers(kTRUE,kTRUE); + psiPrimeJpsi->SetCheckBothChargesGrandMothers(kTRUE,kTRUE); + die->AddSignalMC(psiPrimeJpsi); + + // physical backgrounds (electrons pairs from other sources and their combinatorics) + // 8 + AliDielectronSignalMC* diEleContinuum = new AliDielectronSignalMC("diEleContinuum","di-electron continuum"); // all di-electrons originating in the collision + diEleContinuum->SetLegPDGs(11,-11); + diEleContinuum->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleContinuum->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleContinuum->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleContinuum); + + // 9 + AliDielectronSignalMC* diEleCharm = new AliDielectronSignalMC("diEleCharm","di-electrons from charm"); // dielectrons originating from charm hadrons (not neccessary from same mother) + diEleCharm->SetLegPDGs(11,-11); + diEleCharm->SetMotherPDGs(403,403); + diEleCharm->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleCharm->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleCharm->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleCharm->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleCharm); + + // 10 + AliDielectronSignalMC* diEleOpenCharm = new AliDielectronSignalMC("diEleOpenCharm","di-electrons from open charm"); // dielectrons originating from open charm hadrons + diEleOpenCharm->SetLegPDGs(11,-11); + diEleOpenCharm->SetMotherPDGs(402,402); + diEleOpenCharm->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleOpenCharm->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleOpenCharm->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleOpenCharm->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleOpenCharm); + + // 11 + AliDielectronSignalMC* diEleOpenCharmJpsi = new AliDielectronSignalMC("diEleOpenCharmJpsi","1 leg from open charm + 1 jpsi leg"); // 1 leg from open charm + 1 leg from jpsi + diEleOpenCharmJpsi->SetLegPDGs(11,-11); + diEleOpenCharmJpsi->SetMotherPDGs(402,443); + diEleOpenCharmJpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleOpenCharmJpsi->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleOpenCharmJpsi->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleOpenCharmJpsi->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleOpenCharmJpsi); + + // 12 + AliDielectronSignalMC* diEleBeauty = new AliDielectronSignalMC("diEleBeauty","di-electrons from beauty"); // dielectrons originating from beauty hadrons (not neccessary from same mother) + diEleBeauty->SetLegPDGs(11,-11); + diEleBeauty->SetMotherPDGs(503,503); + diEleBeauty->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleBeauty->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleBeauty->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleBeauty->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleBeauty); + + // 13 + AliDielectronSignalMC* diEleOpenBeauty = new AliDielectronSignalMC("diEleOpenBeauty","di-electrons from open beauty"); // dielectrons originating from open beauty hadrons + diEleOpenBeauty->SetLegPDGs(11,-11); + diEleOpenBeauty->SetMotherPDGs(502,502); + diEleOpenBeauty->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleOpenBeauty->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleOpenBeauty->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleOpenBeauty->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleOpenBeauty); + + // 14 + AliDielectronSignalMC* diEleBeautyJpsi = new AliDielectronSignalMC("diEleBeautyJpsi","1 leg from beauty + 1 jpsi leg"); // 1 leg from beauty + 1 leg from jpsi + diEleBeautyJpsi->SetLegPDGs(11,-11); + diEleBeautyJpsi->SetMotherPDGs(503,443); + diEleBeautyJpsi->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleBeautyJpsi->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleBeautyJpsi->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleBeautyJpsi->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleBeautyJpsi); + + // 15 + AliDielectronSignalMC* diEleBeautyOpenCharm = new AliDielectronSignalMC("diEleBeautyOpenCharm","1 leg from beauty + 1 leg from open charm"); // 1 leg from open charm + 1 leg from beauty + diEleBeautyOpenCharm->SetLegPDGs(11,-11); + diEleBeautyOpenCharm->SetMotherPDGs(503,403); + diEleBeautyOpenCharm->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleBeautyOpenCharm->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleBeautyOpenCharm->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleBeautyOpenCharm->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleBeautyOpenCharm); + + // 16 + AliDielectronSignalMC* diEleStrange = new AliDielectronSignalMC("diEleStrange","di-electrons from strange particles"); // dielectrons originating from strange hadrons (not neccessary from same mother) + diEleStrange->SetLegPDGs(11,-11); + diEleStrange->SetMotherPDGs(300,300); + diEleStrange->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleStrange->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleStrange->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleStrange->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleStrange); + + // 17 + AliDielectronSignalMC* diEleStrangeCharm = new AliDielectronSignalMC("diEleStrangeCharm","1 leg from strange + 1 leg from charm"); // 1 leg from strange + 1 leg from charm + diEleStrangeCharm->SetLegPDGs(11,-11); + diEleStrangeCharm->SetMotherPDGs(300,403); + diEleStrangeCharm->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleStrangeCharm->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleStrangeCharm->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleStrangeCharm->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleStrangeCharm); + + // 18 + AliDielectronSignalMC* diEleStrangeBeauty = new AliDielectronSignalMC("diEleStrangeBeauty","1 leg from strange + 1 leg from beauty"); // 1 leg from strange + 1 leg from beauty + diEleStrangeBeauty->SetLegPDGs(11,-11); + diEleStrangeBeauty->SetMotherPDGs(300,503); + diEleStrangeBeauty->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleStrangeBeauty->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleStrangeBeauty->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleStrangeBeauty->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleStrangeBeauty); + + // 19 + AliDielectronSignalMC* diEleLight1 = new AliDielectronSignalMC("diEleLight1","1 leg from light particles + 1 leg from everything"); // 1 leg from light hadrons (100+ PYTHIA codes) + 1 leg from everything + diEleLight1->SetLegPDGs(11,-11); + diEleLight1->SetMotherPDGs(100,0); + diEleLight1->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleLight1->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleLight1->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleLight1->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleLight1); + + // 20 + AliDielectronSignalMC* diEleLight1Charm = new AliDielectronSignalMC("diEleLight1Charm","1 leg from light particles + 1 leg from charm"); // 1 leg from light hadrons (100+ PYTHIA codes) + 1 leg from charm + diEleLight1Charm->SetLegPDGs(11,-11); + diEleLight1Charm->SetMotherPDGs(100,403); + diEleLight1Charm->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleLight1Charm->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleLight1Charm->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleLight1Charm->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleLight1Charm); + + // 21 + AliDielectronSignalMC* diEleLight1Beauty = new AliDielectronSignalMC("diEleLight1Beauty","1 leg from light particles + 1 leg from beauty"); // 1 leg from light hadrons (100+ PYTHIA codes) + 1 leg from beauty + diEleLight1Beauty->SetLegPDGs(11,-11); + diEleLight1Beauty->SetMotherPDGs(100,503); + diEleLight1Beauty->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleLight1Beauty->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleLight1Beauty->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleLight1Beauty->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleLight1Beauty); + + // 22 + AliDielectronSignalMC* diEleLight2 = new AliDielectronSignalMC("diEleLight2","1 leg from light particles + 1 leg from everything"); // 1 leg from light hadrons (200+ PYTHIA codes) + 1 leg from everything + diEleLight2->SetLegPDGs(11,-11); + diEleLight2->SetMotherPDGs(200,0); + diEleLight2->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleLight2->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleLight2->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleLight2->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleLight2); + + // 23 + AliDielectronSignalMC* diEleLight2Charm = new AliDielectronSignalMC("diEleLight2Charm","1 leg from light particles + 1 leg from charm"); // 1 leg from light hadrons (200+ PYTHIA codes) + 1 leg from charm + diEleLight2Charm->SetLegPDGs(11,-11); + diEleLight2Charm->SetMotherPDGs(100,403); + diEleLight2Charm->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleLight2Charm->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleLight2Charm->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleLight2Charm->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleLight2Charm); + + // 24 + AliDielectronSignalMC* diEleLight2Beauty = new AliDielectronSignalMC("diEleLight2Beauty","1 leg from light particles + 1 leg from beauty"); // 1 leg from light hadrons (100+ PYTHIA codes) + 1 leg from beauty + diEleLight2Beauty->SetLegPDGs(11,-11); + diEleLight2Beauty->SetMotherPDGs(100,503); + diEleLight2Beauty->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kFinalState); + diEleLight2Beauty->SetCheckBothChargesLegs(kTRUE,kTRUE); + diEleLight2Beauty->SetCheckBothChargesMothers(kTRUE,kTRUE); + diEleLight2Beauty->SetFillPureMCStep(kTRUE); + die->AddSignalMC(diEleLight2Beauty); + + // background from secondary electrons + // 25 + AliDielectronSignalMC* secondaryElectrons = new AliDielectronSignalMC("secondaryElectrons","Secondary electrons"); // all di-electrons from secondary electrons (interaction with detector) + secondaryElectrons->SetLegPDGs(11,-11); + secondaryElectrons->SetLegSources(AliDielectronSignalMC::kSecondary, AliDielectronSignalMC::kSecondary); + secondaryElectrons->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(secondaryElectrons); + + // 26 + AliDielectronSignalMC* primarySecElePairs = new AliDielectronSignalMC("primarySecElePairs","Primary+Secondary electron pairs"); // primary-secondary pairs + primarySecElePairs->SetLegPDGs(11,-11); + primarySecElePairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + primarySecElePairs->SetLegSources(AliDielectronSignalMC::kFinalState, AliDielectronSignalMC::kSecondary); + die->AddSignalMC(primarySecElePairs); + + // 27 + AliDielectronSignalMC* conversionElePairs = new AliDielectronSignalMC("conversionElePairs","conversion electron pairs"); // pairs made from conversion (may be also from 2 different conversions) + conversionElePairs->SetLegPDGs(11,-11); + conversionElePairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + conversionElePairs->SetLegSources(AliDielectronSignalMC::kSecondary, AliDielectronSignalMC::kSecondary); + conversionElePairs->SetMotherPDGs(22,22); + die->AddSignalMC(conversionElePairs); + + // misidentification + // 28 + AliDielectronSignalMC* allEleMisIdPairs = new AliDielectronSignalMC("allEleMisIdPairs","all electron+misid. pairs"); // one true electron + a mis-id electron (all sources included) + allEleMisIdPairs->SetLegPDGs(11,11,kFALSE,kTRUE); + allEleMisIdPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(allEleMisIdPairs); + + // 29 + AliDielectronSignalMC* allMisIdMisIdPairs = new AliDielectronSignalMC("allMisIdMisIdPairs","all misid.+misid. pairs"); // mis-id + mis-id + allMisIdMisIdPairs->SetLegPDGs(11,11,kTRUE,kTRUE); + allMisIdMisIdPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(allMisIdMisIdPairs); + + // 30 + AliDielectronSignalMC* elePionPairs = new AliDielectronSignalMC("elePionPairs","electron+pion pairs"); // true electron + mis-id pion + elePionPairs->SetLegPDGs(11,211); + elePionPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(elePionPairs); + + // 31 + AliDielectronSignalMC* eleProtonPairs = new AliDielectronSignalMC("eleProtonPairs","Electron+proton pairs"); // true electron + mis-id proton + eleProtonPairs->SetLegPDGs(11,2212); + eleProtonPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(eleProtonPairs); + + // 32 + AliDielectronSignalMC* eleKaonPairs = new AliDielectronSignalMC("eleKaonPairs","electron+kaon pairs"); // true electron + mis-id kaon + eleKaonPairs->SetLegPDGs(11,321); + eleKaonPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(eleKaonPairs); + + // 33 + AliDielectronSignalMC* eleProtonPairs = new AliDielectronSignalMC("eleProtonPairs","Electron+proton pairs"); // true electron + mis-id proton + eleProtonPairs->SetLegPDGs(11,2212); + eleProtonPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(eleProtonPairs); + + // 34 + AliDielectronSignalMC* piPiPairs = new AliDielectronSignalMC("piPiPairs","pion+pion pairs"); // mis-id pion + mis-id pion + piPiPairs->SetLegPDGs(211,211); + piPiPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(piPiPairs); + + // 35 + AliDielectronSignalMC* piKaonPairs = new AliDielectronSignalMC("piKaonPairs","pion+kaon pairs"); // mis-id pion + mis-id kaon + piKaonPairs->SetLegPDGs(211,321); + piKaonPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(piKaonPairs); + + // 36 + AliDielectronSignalMC* piProtonPairs = new AliDielectronSignalMC("piProtonPairs","pion+proton pairs"); // mis-id pion + mis-id proton + piProtonPairs->SetLegPDGs(211,2212); + piProtonPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(piProtonPairs); + + // 37 + AliDielectronSignalMC* kaonKaonPairs = new AliDielectronSignalMC("kaonKaonPairs","kaon+kaon pairs"); // mis-id kaon + mis-id kaon + kaonKaonPairs->SetLegPDGs(321,321); + kaonKaonPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(kaonKaonPairs); + + // 38 + AliDielectronSignalMC* muonAllPairs = new AliDielectronSignalMC("muonAllPairs","muon+everything pairs"); // mis-id muon + something else (electron, pion, kaon, proton) + muonAllPairs->SetLegPDGs(13,13,kFALSE,kTRUE); + muonAllPairs->SetCheckBothChargesLegs(kTRUE,kTRUE); + die->AddSignalMC(muonAllPairs); +} + +//______________________________________________________________________________________ +void SetEtaCorrection() +{ +// +// Eta equalization for the TPC response +// + if (AliDielectronPID::GetEtaCorrFunction()) return; + + TString list=gSystem->Getenv("LIST"); + + TFile f("$TRAIN_ROOT/jpsi_JPSI/EtaCorrMaps.root"); + if (!f.IsOpen()) return; + TList *keys=f.GetListOfKeys(); + + for (Int_t i=0; iGetEntries(); ++i){ + TString kName=keys->At(i)->GetName(); + TPRegexp reg(kName); + if (reg.MatchB(list)){ + printf("Using Eta Correction Function: %s\n",kName.Data()); + AliDielectronPID::SetEtaCorrFunction((TF1*)f.Get(kName.Data())); + } + } +} diff --git a/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliAnalysisTaskCorrelationTree.cxx b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliAnalysisTaskCorrelationTree.cxx new file mode 100644 index 00000000000..04587d5a6c7 --- /dev/null +++ b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliAnalysisTaskCorrelationTree.cxx @@ -0,0 +1,711 @@ +/************************************************************************* +* 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. * +**************************************************************************/ + +/////////////////////////////////////////////////////////////////////////// +// // +// Analysis task for creating a reduced data tree // +// // +/////////////////////////////////////////////////////////////////////////// + +#include +using namespace std; + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include "AliDielectron.h" +#include "AliDielectronHistos.h" +#include "AliDielectronMC.h" +#include "AliDielectronVarManager.h" +#include "AliFlowTrackCuts.h" +#include "AliFlowBayesianPID.h" + +#include "AliCorrelationReducedEvent.h" +#include "AliAnalysisTaskCorrelationTree.h" + +ClassImp(AliAnalysisTaskCorrelationTree) + + +//_________________________________________________________________________________ +AliAnalysisTaskCorrelationTree::AliAnalysisTaskCorrelationTree() : + AliAnalysisTaskSE(), + fListDielectron(), + fListHistos(), + fSelectPhysics(kFALSE), + fTriggerMask(AliVEvent::kAny), + fRejectPileup(kFALSE), + fEventFilter(0x0), + fTrackFilter(0x0), + fFlowTrackFilter(0x0), + fK0sCuts(0x0), + fLambdaCuts(0x0), + fK0sPionCuts(0x0), + fLambdaProtonCuts(0x0), + fLambdaPionCuts(0x0), + fK0sMassRange(), + fLambdaMassRange(), + fV0Histos(0x0), + fTreeFile(0x0), + fTree(0x0), + fFriendTreeFile(0x0), + fFriendTree(0x0), + fReducedEvent(0x0), + fReducedEventFriend(0x0), + fFlowTrackCuts(0x0) +{ + // + // Constructor + // +} + +//_________________________________________________________________________________ +AliAnalysisTaskCorrelationTree::AliAnalysisTaskCorrelationTree(const char *name) : + AliAnalysisTaskSE(name), + fListDielectron(), + fListHistos(), + fSelectPhysics(kFALSE), + fTriggerMask(AliVEvent::kAny), + fRejectPileup(kFALSE), + fEventFilter(0x0), + fTrackFilter(0x0), + fFlowTrackFilter(0x0), + fK0sCuts(0x0), + fLambdaCuts(0x0), + fK0sPionCuts(0x0), + fLambdaProtonCuts(0x0), + fLambdaPionCuts(0x0), + fK0sMassRange(), + fLambdaMassRange(), + fV0Histos(0x0), + fTreeFile(0x0), + fTree(0x0), + fFriendTreeFile(0x0), + fFriendTree(0x0), + fReducedEvent(0x0), + fReducedEventFriend(0x0), + fFlowTrackCuts(0x0) +{ + // + // Constructor + // + fK0sMassRange[0] = 0.4; fK0sMassRange[1] = 0.6; + fLambdaMassRange[0] = 1.08; fLambdaMassRange[1] = 1.15; + + DefineInput(0,TChain::Class()); + DefineOutput(1, TList::Class()); // QA histograms + //DefineOutput(2, TTree::Class()); // reduced information tree + //DefineOutput(3, TTree::Class()); // reduced information tree with friends + DefineOutput(2, TTree::Class()); // reduced information tree with friends + + fListHistos.SetName("QAhistograms"); + fListDielectron.SetOwner(); + fListHistos.SetOwner(kFALSE); +} + + +//_________________________________________________________________________________ +void AliAnalysisTaskCorrelationTree::UserCreateOutputObjects() +{ + // + // Add all histogram manager histogram lists to the output TList + // + + if (!fListHistos.IsEmpty() || fTree || fFriendTree) return; //already initialised + + TIter nextDie(&fListDielectron); + AliDielectron *die=0; + while ( (die=static_cast(nextDie())) ){ + die->Init(); + if (die->GetHistogramList()) fListHistos.Add(const_cast(die->GetHistogramList())); + } + if(fV0Histos) fListHistos.Add(const_cast(fV0Histos->GetHistogramList())); + + fFriendTree = new TTree("DstFriendTree","Reduced ESD information"); + fReducedEventFriend = new AliCorrelationReducedEventFriend(); + fFriendTree->Branch("Event",&fReducedEventFriend,16000,99); + + fTreeFile = new TFile("dstTree.root", "RECREATE"); + fTree = new TTree("DstTree","Reduced ESD information"); + fReducedEvent = new AliCorrelationReducedEvent(); + fTree->Branch("Event",&fReducedEvent,16000,99); + + fFlowTrackCuts = new AliFlowTrackCuts("flow cuts"); + fFlowTrackCuts->SetPID(AliPID::kPion, AliFlowTrackCuts::kTOFbayesian); + + PostData(1, &fListHistos); + //PostData(2, fTree); + //PostData(3, fFriendTree); + PostData(2, fFriendTree); +} + +//_________________________________________________________________________________ +void AliAnalysisTaskCorrelationTree::UserExec(Option_t *option) +{ + // + // Main loop. Called for every event + // + //cout << "Event" << endl; + option = option; + AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager(); + Bool_t isESD=man->GetInputEventHandler()->IsA()==AliESDInputHandler::Class(); + Bool_t isAOD=man->GetInputEventHandler()->IsA()==AliAODInputHandler::Class(); + + AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler()); + if (!inputHandler) return; + + if ( inputHandler->GetPIDResponse() ){ + AliDielectronVarManager::SetPIDResponse( inputHandler->GetPIDResponse() ); + } else { + AliFatal("This task needs the PID response attached to the input event handler!"); + } + + fFlowTrackCuts->SetEvent(InputEvent()); + + // Was event selected ? + UInt_t isSelected = AliVEvent::kAny; + if(fSelectPhysics && inputHandler){ + if((isESD && inputHandler->GetEventSelection()) || isAOD){ + isSelected = inputHandler->IsEventSelected(); + isSelected&=fTriggerMask; + } + } + + if(isSelected==0) { + return; + } + + // fill event histograms before event filter + Double_t values[AliDielectronVarManager::kNMaxValues]={0}; + AliDielectronVarManager::Fill(InputEvent(),values); + + TIter nextDie(&fListDielectron); + AliDielectron *die=0; + while ( (die=static_cast(nextDie())) ){ + AliDielectronHistos *h=die->GetHistoManager(); + if (h){ + if (h->GetHistogramList()->FindObject("Event_noCuts")) + h->FillClass("Event_noCuts",AliDielectronVarManager::kNMaxValues,values); + } + } + nextDie.Reset(); + + //event filter + if (fEventFilter) { + if (!fEventFilter->IsSelected(InputEvent())) return; + } + + //pileup + if (fRejectPileup){ + if (InputEvent()->IsPileupFromSPD(3,0.8,3.,2.,5.)) return; + } + + //cout << "Event selected" << endl; + //bz for AliKF + Double_t bz = InputEvent()->GetMagneticField(); + AliKFParticle::SetField( bz ); + + //Process event in all AliDielectron instances + fReducedEvent->ClearEvent(); + fReducedEventFriend->ClearEvent(); + FillEventInfo(); + FillV0PairInfo(); + + Short_t idie=0; + while((die=static_cast(nextDie()))){ + die->Process(InputEvent()); + FillDielectronPairInfo(die, idie); + ++idie; + } + nextDie.Reset(); + + FillTrackInfo(); + FillFriendEventInfo(); // Q-vector calculation + + fTree->Fill(); + fFriendTree->Fill(); + + // if there are candidate pairs, add the information to the reduced tree + PostData(1, &fListHistos); + //PostData(2, fTree); + //PostData(3, fFriendTree); + PostData(2, fFriendTree); +} + + +//_________________________________________________________________________________ +void AliAnalysisTaskCorrelationTree::FillEventInfo() +{ + // + // fill reduced event information + // + AliVEvent* event = InputEvent(); + // Was event selected ? + AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager(); + Bool_t isESD = (event->IsA()==AliESDEvent::Class()); + Bool_t isAOD = (event->IsA()==AliAODEvent::Class()); + + AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler()); + UInt_t isSelected = AliVEvent::kAny; + if(inputHandler){ + if((isESD && inputHandler->GetEventSelection()) || isAOD){ + isSelected = inputHandler->IsEventSelected(); + isSelected&=fTriggerMask; + } + } + + Double_t values[AliDielectronVarManager::kNMaxValues]; + AliDielectronVarManager::Fill(event, values); + + fReducedEvent->fRunNo = event->GetRunNumber(); + fReducedEvent->fBC = event->GetBunchCrossNumber(); + fReducedEvent->fTriggerMask = event->GetTriggerMask(); + fReducedEvent->fIsPhysicsSelection = (isSelected!=0 ? kTRUE : kFALSE); + AliVVertex* eventVtx = 0x0; + if(isESD) eventVtx = const_cast((static_cast(event))->GetPrimaryVertexTracks()); + if(isAOD) eventVtx = const_cast((static_cast(event))->GetPrimaryVertex()); + if(eventVtx) { + fReducedEvent->fVtx[0] = (isESD ? ((AliESDVertex*)eventVtx)->GetXv() : ((AliAODVertex*)eventVtx)->GetX()); + fReducedEvent->fVtx[1] = (isESD ? ((AliESDVertex*)eventVtx)->GetYv() : ((AliAODVertex*)eventVtx)->GetY()); + fReducedEvent->fVtx[2] = (isESD ? ((AliESDVertex*)eventVtx)->GetZv() : ((AliAODVertex*)eventVtx)->GetZ()); + fReducedEvent->fNVtxContributors = eventVtx->GetNContributors(); + } + if(isESD) { + eventVtx = const_cast((static_cast(event))->GetPrimaryVertexTPC()); + if(eventVtx) { + fReducedEvent->fVtxTPC[0] = ((AliESDVertex*)eventVtx)->GetXv(); + fReducedEvent->fVtxTPC[1] = ((AliESDVertex*)eventVtx)->GetYv(); + fReducedEvent->fVtxTPC[2] = ((AliESDVertex*)eventVtx)->GetZv(); + fReducedEvent->fNVtxTPCContributors = eventVtx->GetNContributors(); + } + } + + AliCentrality *centrality = event->GetCentrality(); + if(centrality) { + fReducedEvent->fCentrality[0] = centrality->GetCentralityPercentile("V0M"); + fReducedEvent->fCentrality[1] = centrality->GetCentralityPercentile("CL1"); + fReducedEvent->fCentrality[2] = centrality->GetCentralityPercentile("TRK"); + fReducedEvent->fCentrality[3] = centrality->GetCentralityPercentile("ZEMvsZDC"); + fReducedEvent->fCentQuality = centrality->GetQuality(); + } + + fReducedEvent->fNtracks[0] = event->GetNumberOfTracks(); + fReducedEvent->fSPDntracklets = values[AliDielectronVarManager::kNaccTrckltsEsd10Corr]; + + AliVVZERO* vzero = event->GetVZEROData(); + for(Int_t i=0;i<64;++i) + fReducedEvent->fVZEROMult[i] = vzero->GetMultiplicity(i); + + AliESDZDC* zdc = (isESD ? (static_cast(event))->GetESDZDC() : 0x0); + if(zdc) { + for(Int_t i=0; i<4; ++i) fReducedEvent->fZDCnEnergy[i] = zdc->GetZN1TowerEnergy()[i]; + for(Int_t i=4; i<8; ++i) fReducedEvent->fZDCnEnergy[i] = zdc->GetZN2TowerEnergy()[i-5]; + } + + // EMCAL/PHOS clusters + FillCaloClusters(); + + // TODO FMD multiplicities + +} + + +//_________________________________________________________________________________ +void AliAnalysisTaskCorrelationTree::FillCaloClusters() { + // + // Fill info about the calorimeter clusters + // + AliVEvent* event = InputEvent(); + Int_t nclusters = event->GetNumberOfCaloClusters(); + + fReducedEvent->fNCaloClusters = 0; + for(Int_t iclus=0; iclusGetCaloCluster(iclus); + + TClonesArray& clusters = *(fReducedEvent->fCaloClusters); + AliCorrelationReducedCaloCluster *reducedCluster=new(clusters[fReducedEvent->fNCaloClusters]) AliCorrelationReducedCaloCluster(); + + reducedCluster->fType = (cluster->IsEMCAL() ? AliCorrelationReducedCaloCluster::kEMCAL : AliCorrelationReducedCaloCluster::kPHOS); + reducedCluster->fEnergy = cluster->E(); + reducedCluster->fTrackDx = cluster->GetTrackDx(); + reducedCluster->fTrackDz = cluster->GetTrackDz(); + fReducedEvent->fNCaloClusters += 1; + } // end loop over clusters +} + + +//_________________________________________________________________________________ +void AliAnalysisTaskCorrelationTree::FillFriendEventInfo() { + // + // Fill event info into the friend tree + // + // Add here calculated Q-vector components from all detectors + for(Int_t idet=0; idetGetQvector(fReducedEventFriend->fQvector[idet], idet); + for(Int_t ih=0; ihfEventPlaneStatus[idet][ih] = AliCorrelationReducedEventFriend::kRaw; + } +} + + +//_________________________________________________________________________________ +void AliAnalysisTaskCorrelationTree::FillTrackInfo() +{ + // + // fill reduced track information + // + AliVEvent* event = InputEvent(); + Bool_t isESD = (event->IsA()==AliESDEvent::Class()); + Bool_t isAOD = (event->IsA()==AliAODEvent::Class()); + + Int_t ntracks=event->GetNumberOfTracks(); + for(Int_t itrack=0; itrackGetTrack(itrack); + //apply track cuts + if(fTrackFilter && !fTrackFilter->IsSelected(particle)) continue; + + TClonesArray& tracks = *(fReducedEvent->fTracks); + AliCorrelationReducedTrack *reducedParticle=new(tracks[fReducedEvent->fNtracks[1]]) AliCorrelationReducedTrack(); + + Double_t values[AliDielectronVarManager::kNMaxValues]; + AliDielectronVarManager::Fill(particle, values); + reducedParticle->fStatus = (ULong_t)values[AliDielectronVarManager::kTrackStatus]; + reducedParticle->fGlobalPhi = values[AliDielectronVarManager::kPhi]; + reducedParticle->fGlobalPt = values[AliDielectronVarManager::kPt]*values[AliDielectronVarManager::kCharge]; + reducedParticle->fGlobalEta = values[AliDielectronVarManager::kEta]; + reducedParticle->fMomentumInner = values[AliDielectronVarManager::kPIn]; + reducedParticle->fDCA[0] = values[AliDielectronVarManager::kImpactParXY]; + reducedParticle->fDCA[1] = values[AliDielectronVarManager::kImpactParZ]; + + reducedParticle->fITSclusterMap = values[AliDielectronVarManager::kITSclusterMap]; + reducedParticle->fITSsignal = values[AliDielectronVarManager::kITSsignal]; + + reducedParticle->fTPCNcls = (UChar_t)values[AliDielectronVarManager::kNclsTPC]; + reducedParticle->fTPCNclsF = (UChar_t)values[AliDielectronVarManager::kNFclsTPC]; + reducedParticle->fTPCNclsIter1 = (UChar_t)values[AliDielectronVarManager::kNclsTPCiter1]; + reducedParticle->fTPCsignal = values[AliDielectronVarManager::kTPCsignal]; + reducedParticle->fTPCnSig[0] = values[AliDielectronVarManager::kTPCnSigmaEle]; + reducedParticle->fTPCnSig[1] = values[AliDielectronVarManager::kTPCnSigmaPio]; + reducedParticle->fTPCnSig[2] = values[AliDielectronVarManager::kTPCnSigmaKao]; + reducedParticle->fTPCnSig[3] = values[AliDielectronVarManager::kTPCnSigmaPro]; + + reducedParticle->fTOFbeta = values[AliDielectronVarManager::kTOFbeta]; + reducedParticle->fTOFnSig[0] = values[AliDielectronVarManager::kTOFnSigmaEle]; + reducedParticle->fTOFnSig[1] = values[AliDielectronVarManager::kTOFnSigmaPio]; + reducedParticle->fTOFnSig[2] = values[AliDielectronVarManager::kTOFnSigmaKao]; + reducedParticle->fTOFnSig[3] = values[AliDielectronVarManager::kTOFnSigmaPro]; + + reducedParticle->fTRDpid[0] = values[AliDielectronVarManager::kTRDprobEle]; + reducedParticle->fTRDpid[1] = values[AliDielectronVarManager::kTRDprobPio]; + + if(fFlowTrackFilter) { + // switch on the first bit if this particle should be used for the event plane + if(fFlowTrackFilter->IsSelected(particle)) reducedParticle->fFlags |= (1<<0); + } + + if(isESD){ + AliESDtrack *track=static_cast(particle); + reducedParticle->fTrackId = (UShort_t)track->GetID(); + reducedParticle->fTPCCrossedRows = (UChar_t)track->GetTPCCrossedRows(); + reducedParticle->fTPCClusterMap = EncodeTPCClusterMap(track); + const AliExternalTrackParam* tpcInner = track->GetTPCInnerParam(); + reducedParticle->fTPCPhi = (tpcInner ? tpcInner->Phi() : 0.0); + reducedParticle->fTPCPt = (tpcInner ? tpcInner->Pt() : 0.0); + reducedParticle->fTPCEta = (tpcInner ? tpcInner->Eta() : 0.0); + reducedParticle->fTRDntracklets[0] = track->GetTRDntracklets(); + reducedParticle->fTRDntracklets[1] = track->GetTRDntrackletsPID(); + fFlowTrackCuts->IsSelected(track); + AliFlowBayesianPID* bayesResponse = fFlowTrackCuts->GetBayesianResponse(); + reducedParticle->fBayesPID[0] = bayesResponse->GetProb()[AliPID::kPion]; + reducedParticle->fBayesPID[1] = bayesResponse->GetProb()[AliPID::kKaon]; + reducedParticle->fBayesPID[2] = bayesResponse->GetProb()[AliPID::kProton]; + if(track->IsEMCAL()) reducedParticle->fCaloClusterId = track->GetEMCALcluster(); + if(track->IsPHOS()) reducedParticle->fCaloClusterId = track->GetPHOScluster(); + } + if(isAOD) { + AliAODTrack *track=static_cast(particle); + reducedParticle->fTrackId = track->GetID(); + if(track->IsEMCAL()) reducedParticle->fCaloClusterId = track->GetEMCALcluster(); + if(track->IsPHOS()) reducedParticle->fCaloClusterId = track->GetPHOScluster(); + } + + fReducedEvent->fNtracks[1] += 1; + } +} + + +//_________________________________________________________________________________ +void AliAnalysisTaskCorrelationTree::FillDielectronPairInfo(AliDielectron* die, Short_t iDie) +{ + // + // fill reduced pair information + // + Bool_t hasMC=AliDielectronMC::Instance()->HasMC(); + + for(Int_t iType=0; iType<3; ++iType) { + + const TObjArray* array = die->GetPairArray(iType); + if(!array || array->GetEntriesFast()==0) continue; + cout << "dielectron tracks (idie/type/n): " << iDie << " / 0 / " << die->GetTrackArray(0)->GetEntriesFast() << endl; + cout << "dielectron tracks (idie/type/n): " << iDie << " / 1 / " << die->GetTrackArray(1)->GetEntriesFast() << endl; + + cout << "diele idie/type/n-candidates: " << iDie << " / " << iType << " / " << array->GetEntriesFast() << endl; + for(Int_t iCandidate=0; iCandidateGetEntriesFast(); ++iCandidate) { + if(iCandidate%100==0) cout << "iCandidate = " << iCandidate << endl; + AliDielectronPair* pair = (AliDielectronPair*)array->At(iCandidate); + Double_t values[AliDielectronVarManager::kNMaxValues]; + AliDielectronVarManager::Fill(pair, values); + + TClonesArray& tracks = *(fReducedEvent->fCandidates); + AliCorrelationReducedPair *reducedParticle= + new (tracks[fReducedEvent->fNV0candidates[1]+fReducedEvent->fNDielectronCandidates]) AliCorrelationReducedPair(); + // !!! hardcoded flag for dielectron id + reducedParticle->fCandidateId = (iDie==0 ? AliCorrelationReducedPair::kJpsiToEE : AliCorrelationReducedPair::kPhiToKK); + reducedParticle->fPairType = values[AliDielectronVarManager::kPairType]; + reducedParticle->fLegIds[0] = (UShort_t)(static_cast(pair->GetFirstDaughter()))->GetID(); + reducedParticle->fLegIds[1] = (UShort_t)(static_cast(pair->GetSecondDaughter()))->GetID(); + reducedParticle->fMass[0] = values[AliDielectronVarManager::kM]; + reducedParticle->fMass[1] = -999.; + reducedParticle->fMass[2] = -999.; + reducedParticle->fPhi = values[AliDielectronVarManager::kPhi]; // in the [-pi,pi] interval + if(reducedParticle->fPhi<0.0) reducedParticle->fPhi = 2.0*TMath::Pi() + reducedParticle->fPhi; // converted to [0,2pi] + reducedParticle->fPt = values[AliDielectronVarManager::kPt]; + reducedParticle->fEta = values[AliDielectronVarManager::kEta]; + reducedParticle->fLxy = values[AliDielectronVarManager::kPseudoProperTime]; + reducedParticle->fOpeningAngle = values[AliDielectronVarManager::kOpeningAngle]; + + reducedParticle->fMCid = 0; + if(hasMC) { + AliDielectronMC::Instance()->ConnectMCEvent(); + const TObjArray* mcSignals = die->GetMCSignals(); + for(Int_t iSig=0; iSigGetEntries(); ++iSig) { + if(iSig>31) break; + AliDielectronMC *mc=AliDielectronMC::Instance(); + if(mc->IsMCTruth(pair, (AliDielectronSignalMC*)mcSignals->At(iSig))) { + reducedParticle->fMCid = reducedParticle->fMCid | (1<fNDielectronCandidates += 1; + } // end loop over candidates + } // end loop over pair type +} + + +//_________________________________________________________________________________ +void AliAnalysisTaskCorrelationTree::FillV0PairInfo() +{ + // + // fill reduced pair information + // + AliESDEvent* esd = (AliESDEvent*)InputEvent(); + const AliESDVertex *primaryVertex = esd->GetPrimaryVertex(); + AliKFVertex primaryVertexKF(*primaryVertex); + + fReducedEvent->fNV0candidates[0] = InputEvent()->GetNumberOfV0s(); + + Double_t valuesPos[AliDielectronVarManager::kNMaxValues]; + Double_t valuesNeg[AliDielectronVarManager::kNMaxValues]; + + cout << "n total V0s: " << InputEvent()->GetNumberOfV0s() << endl; + for(Int_t iV0=0; iV0GetNumberOfV0s(); ++iV0) { // loop over V0s + if(iV0%1000==0) cout << "iV0 = " << iV0 << endl; + AliESDv0 *v0 = esd->GetV0(iV0); + + AliESDtrack* legPos = esd->GetTrack(v0->GetPindex()); + AliESDtrack* legNeg = esd->GetTrack(v0->GetNindex()); + + if(legPos->GetSign() == legNeg->GetSign()) { + //cout << "V0 rejected: same sign legs" << endl; + continue; + } + + Bool_t v0ChargesAreCorrect = (legPos->GetSign()==+1 ? kTRUE : kFALSE); + legPos = (!v0ChargesAreCorrect ? esd->GetTrack(v0->GetNindex()) : legPos); + legNeg = (!v0ChargesAreCorrect ? esd->GetTrack(v0->GetPindex()) : legNeg); + + // apply the K0s filter + Bool_t goodK0s = kTRUE; + if(fK0sPionCuts && (!fK0sPionCuts->IsSelected(legPos) || !fK0sPionCuts->IsSelected(legNeg))) goodK0s = kFALSE; + if(fK0sCuts) { + TList k0sList; + k0sList.Add(v0); + k0sList.Add(legPos); k0sList.Add(legNeg); + k0sList.Add(const_cast(primaryVertex)); + if(!fK0sCuts->IsSelected(&k0sList)) goodK0s = kFALSE; + } + + // apply the lambda filter + Bool_t goodLambda = kTRUE; + if(fLambdaProtonCuts && !fLambdaProtonCuts->IsSelected(legPos)) goodLambda = kFALSE; + if(fLambdaPionCuts && !fLambdaPionCuts->IsSelected(legNeg)) goodLambda = kFALSE; + Bool_t goodALambda = kTRUE; + if(fLambdaProtonCuts && !fLambdaProtonCuts->IsSelected(legNeg)) goodALambda = kFALSE; + if(fLambdaPionCuts && !fLambdaPionCuts->IsSelected(legPos)) goodALambda = kFALSE; + if(fLambdaCuts) { + TList lambdaList; + lambdaList.Add(v0); + lambdaList.Add(legPos); lambdaList.Add(legNeg); + lambdaList.Add(const_cast(primaryVertex)); + if(!fLambdaCuts->IsSelected(&lambdaList)) {goodLambda = kFALSE; goodALambda = kFALSE;} + } + + if(!(goodK0s || goodLambda || goodALambda)) continue; + + // Fill the V0 information into the tree for 3 hypothesis: K0s, Lambda, Anti-Lambda + AliCorrelationReducedPair* k0sReducedPair = FillV0PairInfo(v0, AliCorrelationReducedPair::kK0sToPiPi, legPos, legNeg, &primaryVertexKF, v0ChargesAreCorrect); + AliCorrelationReducedPair* lambdaReducedPair = FillV0PairInfo(v0, AliCorrelationReducedPair::kLambda0ToPPi, legPos, legNeg, &primaryVertexKF, v0ChargesAreCorrect); + AliCorrelationReducedPair* alambdaReducedPair = FillV0PairInfo(v0, AliCorrelationReducedPair::kALambda0ToPPi, legPos, legNeg, &primaryVertexKF, v0ChargesAreCorrect); + + if(goodK0s && k0sReducedPair->fMass[0]>fK0sMassRange[0] && k0sReducedPair->fMass[0]fCandidates); + AliCorrelationReducedPair *goodK0sPair = new (tracks[fReducedEvent->fNV0candidates[1]]) AliCorrelationReducedPair(*k0sReducedPair); + goodK0sPair->fMass[0] = k0sReducedPair->fMass[0]; + goodK0sPair->fMass[1] = lambdaReducedPair->fMass[0]; + goodK0sPair->fMass[2] = alambdaReducedPair->fMass[0]; + fReducedEvent->fNV0candidates[1] += 1; + } else {goodK0s=kFALSE;} + if(goodLambda && lambdaReducedPair->fMass[0]>fLambdaMassRange[0] && lambdaReducedPair->fMass[0]fCandidates); + AliCorrelationReducedPair *goodLambdaPair = new (tracks[fReducedEvent->fNV0candidates[1]]) AliCorrelationReducedPair(*lambdaReducedPair); + fReducedEvent->fNV0candidates[1] += 1; + goodLambdaPair->fMass[0] = k0sReducedPair->fMass[0]; + goodLambdaPair->fMass[1] = lambdaReducedPair->fMass[0]; + goodLambdaPair->fMass[2] = alambdaReducedPair->fMass[0]; + } else {goodLambda=kFALSE;} + if(goodALambda && alambdaReducedPair->fMass[0]>fLambdaMassRange[0] && alambdaReducedPair->fMass[0]fCandidates); + AliCorrelationReducedPair *goodALambdaPair = new (tracks[fReducedEvent->fNV0candidates[1]]) AliCorrelationReducedPair(*alambdaReducedPair); + fReducedEvent->fNV0candidates[1] += 1; + goodALambdaPair->fMass[0] = k0sReducedPair->fMass[0]; + goodALambdaPair->fMass[1] = lambdaReducedPair->fMass[0]; + goodALambdaPair->fMass[2] = alambdaReducedPair->fMass[0]; + } else {goodALambda = kFALSE;} + delete k0sReducedPair; + delete lambdaReducedPair; + delete alambdaReducedPair; + + if(!(goodK0s || goodLambda || goodALambda)) continue; + + // Fill histograms and the CF container + AliDielectronVarManager::Fill(legPos, valuesPos); + AliDielectronVarManager::Fill(legNeg, valuesNeg); + + if(fV0Histos && fV0Histos->GetHistogramList()->FindObject("V0Track_Pos")) + fV0Histos->FillClass("V0Track_Pos", AliDielectronVarManager::kNMaxValues, valuesPos); + if(fV0Histos && fV0Histos->GetHistogramList()->FindObject("V0Track_Neg")) + fV0Histos->FillClass("V0Track_Neg", AliDielectronVarManager::kNMaxValues, valuesNeg); + } // end loop over V0s +} + + +//_________________________________________________________________________________ +AliCorrelationReducedPair* AliAnalysisTaskCorrelationTree::FillV0PairInfo(AliESDv0* v0, Int_t id, + AliESDtrack* legPos, AliESDtrack* legNeg, + AliKFVertex* vtxKF, Bool_t chargesAreCorrect) { + // + // Create a reduced V0 object and fill it + // + AliCorrelationReducedPair* reducedPair=new AliCorrelationReducedPair(); + reducedPair->fCandidateId = id; + reducedPair->fPairType = v0->GetOnFlyStatus(); // on the fly status + //reducedPair->fOnTheFly = v0->GetOnFlyStatus(); + reducedPair->fLegIds[0] = legPos->GetID(); + reducedPair->fLegIds[1] = legNeg->GetID(); + if(!reducedPair->fPairType) { // offline + UInt_t pidPos = AliPID::kPion; + if(id==AliCorrelationReducedPair::kLambda0ToPPi) pidPos = AliPID::kProton; + UInt_t pidNeg = AliPID::kPion; + if(id==AliCorrelationReducedPair::kALambda0ToPPi) pidNeg = AliPID::kProton; + reducedPair->fMass[0] = v0->GetEffMass(pidPos, pidNeg); + reducedPair->fPhi = v0->Phi(); + if(reducedPair->fPhi<0.0) reducedPair->fPhi = 2.0*TMath::Pi() + reducedPair->fPhi; // converted to [0,2pi] + reducedPair->fPt = v0->Pt(); + reducedPair->fEta = v0->Eta(); + reducedPair->fLxy = 0.0; //TODO + reducedPair->fOpeningAngle = 0.0; + } + else { + const AliExternalTrackParam *negHelix=v0->GetParamN(); + const AliExternalTrackParam *posHelix=v0->GetParamP(); + if(!chargesAreCorrect) { + negHelix = v0->GetParamP(); + posHelix = v0->GetParamN(); + } + AliKFParticle negKF(*(negHelix),(id==AliCorrelationReducedPair::kALambda0ToPPi ? -2212 : -211)); + AliKFParticle posKF(*(posHelix),(id==AliCorrelationReducedPair::kLambda0ToPPi ? +2212 : +211)); + AliKFParticle v0Refit; + v0Refit += negKF; + v0Refit += posKF; + Double_t massFit=0.0, massErrFit=0.0; + v0Refit.GetMass(massFit,massErrFit); + reducedPair->fMass[0] = massFit; + reducedPair->fPhi = v0Refit.GetPhi(); + if(reducedPair->fPhi<0.0) reducedPair->fPhi = 2.0*TMath::Pi() + reducedPair->fPhi; // converted to [0,2pi] + reducedPair->fPt = v0Refit.GetPt(); + reducedPair->fEta = v0Refit.GetEta(); + reducedPair->fLxy = v0Refit.GetPseudoProperDecayTime(*vtxKF, massFit); + reducedPair->fOpeningAngle = negKF.GetAngle(posKF); + } + return reducedPair; +} + + +//_________________________________________________________________________________ +UChar_t AliAnalysisTaskCorrelationTree::EncodeTPCClusterMap(AliESDtrack* track) { + // + // Encode the TPC cluster map into an UChar_t + // Divide the 159 bits from the bit map into 8 groups of adiacent clusters + // For each group enable its corresponding bit if in that group there are more clusters compared to + // a threshold. + // + const UChar_t threshold=5; + TBits tpcClusterMap = track->GetTPCClusterMap(); + UChar_t map=0; + UChar_t n=0; + UChar_t j=0; + for(UChar_t i=0; i<8; ++i) { + n=0; + for(j=i*20; j<(i+1)*20 && j<159; ++j) n+=tpcClusterMap.TestBitNumber(j); + if(n>=threshold) map |= (1<Write(); + fTreeFile->Close(); +} diff --git a/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliAnalysisTaskCorrelationTree.h b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliAnalysisTaskCorrelationTree.h new file mode 100644 index 00000000000..edcb2036e50 --- /dev/null +++ b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliAnalysisTaskCorrelationTree.h @@ -0,0 +1,104 @@ +#ifndef ALIANALYSISTASKCORRELATIONTREE_H +#define ALIANALYSISTASKCORRELATIONTREE_H + +// Analysis task for creating a reduced tree containing event, track and resonance candidate information +// Author: Ionut-Cristian Arsene (i.c.arsene@gsi.de) + +#include "TList.h" + +#include "AliAnalysisTaskSE.h" + +class AliAnalysisCuts; +class TTree; +class TFile; +class AliESDv0Cuts; +class AliKFVertex; +class AliCorrelationReducedEvent; +class AliCorrelationReducedEventFriend; +class AliCorrelationReducedPair; +class AliDielectron; +class AliFlowTrackCuts; + +//_________________________________________________________________________ +class AliAnalysisTaskCorrelationTree : public AliAnalysisTaskSE { + +public: + AliAnalysisTaskCorrelationTree(); + AliAnalysisTaskCorrelationTree(const char *name); + virtual ~AliAnalysisTaskCorrelationTree(){ } + + virtual void UserExec(Option_t *option); + virtual void UserCreateOutputObjects(); + virtual void FinishTaskOutput(); + + void UsePhysicsSelection(Bool_t phy=kTRUE) {fSelectPhysics=phy;} + void SetTriggerMask(UInt_t mask) {fTriggerMask=mask;} + UInt_t GetTriggerMask() const { return fTriggerMask; } + void SetRejectPileup(Bool_t pileup=kTRUE) { fRejectPileup=pileup; } + + // Cuts for selection of event to be written to tree + void SetEventFilter(AliAnalysisCuts * const filter) {fEventFilter=filter;} + // Cuts for selecting tracks included in the tree + void SetTrackFilter(AliAnalysisCuts * const filter) {fTrackFilter=filter;} + // Cuts for selecting tracks to be used for Q vector calculation + void SetFlowTrackFilter(AliAnalysisCuts * const filter) {fFlowTrackFilter = filter;} + + // Cuts for selecting V0s + void SetK0sPionCuts(AliAnalysisCuts * const filter) {fK0sPionCuts=filter;} + void SetLambdaProtonCuts(AliAnalysisCuts * const filter) {fLambdaProtonCuts=filter;} + void SetLambdaPionCuts(AliAnalysisCuts * const filter) {fLambdaPionCuts=filter;} + void SetK0sCuts(AliESDv0Cuts* const cuts) {fK0sCuts = cuts;} + void SetLambdaCuts(AliESDv0Cuts* const cuts) {fLambdaCuts = cuts;} + void SetK0sMassRange(Double_t min=0.4, Double_t max=0.6) {fK0sMassRange[0]=min; fK0sMassRange[1]=max;} + void SetLambdaMassRange(Double_t min=1.08, Double_t max=1.15) {fLambdaMassRange[0]=min; fLambdaMassRange[1]=max;} + void SetV0Histograms(AliDielectronHistos * const histos) {fV0Histos=histos;} + + // Add dielectron objects to the list. These contain cuts and histogram definitions + void AddDielectron(AliDielectron * const die) { fListDielectron.Add(die); } + + private: + + TList fListDielectron; // List of dielectron framework instances + TList fListHistos; //! List of histogram managers in the dielectron framework classes + + Bool_t fSelectPhysics; // Whether to use physics selection + UInt_t fTriggerMask; // Event trigger mask + Bool_t fRejectPileup; // pileup rejection wanted + + AliAnalysisCuts *fEventFilter; // event filter + AliAnalysisCuts *fTrackFilter; // filter for the hadrons to be correlated with the dielectrons + AliAnalysisCuts *fFlowTrackFilter; // filter for the barrel tracks to be used for the Q-vector + + AliESDv0Cuts *fK0sCuts; // v0 standard filter for K0s->pi+pi- + AliESDv0Cuts *fLambdaCuts; // v0 standard filter for Lambda0->p + pi + AliAnalysisCuts *fK0sPionCuts; // filter for pions from K0s + AliAnalysisCuts *fLambdaProtonCuts;// filter for protons from Lambda + AliAnalysisCuts *fLambdaPionCuts; // filter for pions from Lambda + Double_t fK0sMassRange[2]; // mass range for allowed K0s pairs + Double_t fLambdaMassRange[2]; // mass range for allowed Lambda pairs + AliDielectronHistos* fV0Histos; // histogram manager for V0s + + TFile *fTreeFile; // output file containing the tree + TTree *fTree; //! Reduced event tree + TTree *fFriendTreeFile; // output file containing the friend tree + TTree *fFriendTree; //! Reduced event tree with friend info (event plane, etc.) + AliCorrelationReducedEvent *fReducedEvent; // reduced event wise information + AliCorrelationReducedEventFriend *fReducedEventFriend; // friend reduced event wise information + + AliFlowTrackCuts* fFlowTrackCuts; // flow track cuts object + + void FillEventInfo(); // fill reduced event information + void FillFriendEventInfo(); // fill reduced event friend information + void FillTrackInfo(); // fill reduced track information + void FillDielectronPairInfo(AliDielectron* die, Short_t iDie); // fill dielectron reduced pair information + void FillV0PairInfo(); // fill V0 reduced pair information + AliCorrelationReducedPair* FillV0PairInfo(AliESDv0* v0, Int_t id, AliESDtrack* legPos, AliESDtrack* legNeg, AliKFVertex* vtxKF, Bool_t chargesAreCorrect); + UChar_t EncodeTPCClusterMap(AliESDtrack* track); + void FillCaloClusters(); + + AliAnalysisTaskCorrelationTree(const AliAnalysisTaskCorrelationTree &c); + AliAnalysisTaskCorrelationTree& operator= (const AliAnalysisTaskCorrelationTree &c); + + ClassDef(AliAnalysisTaskCorrelationTree, 1); //Analysis Task for creating a reduced event information tree +}; +#endif diff --git a/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliCorrelationReducedEvent.cxx b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliCorrelationReducedEvent.cxx new file mode 100644 index 00000000000..dc6295af07f --- /dev/null +++ b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliCorrelationReducedEvent.cxx @@ -0,0 +1,502 @@ +#ifndef ALICORRELATIONREDUCEDEVENT_H +#include "AliCorrelationReducedEvent.h" +#endif + +#include +#include + +ClassImp(AliCorrelationReducedTrack) +ClassImp(AliCorrelationReducedPair) +ClassImp(AliCorrelationReducedEvent) +ClassImp(AliCorrelationReducedEventFriend) +ClassImp(AliCorrelationReducedCaloCluster) + +TClonesArray* AliCorrelationReducedEvent::fgTracks = 0; +TClonesArray* AliCorrelationReducedEvent::fgCandidates = 0; +TClonesArray* AliCorrelationReducedEvent::fgCaloClusters = 0; + +//_______________________________________________________________________________ +AliCorrelationReducedTrack::AliCorrelationReducedTrack() : + fTrackId(-1), + fStatus(0), + fGlobalPhi(0.0), + fGlobalPt(0.0), + fGlobalEta(0.0), + fTPCPhi(0.0), + fTPCPt(0.0), + fTPCEta(0.0), + fMomentumInner(0.0), + fDCA(), + fITSclusterMap(0), + fITSsignal(0.0), + fTPCNcls(0), + fTPCCrossedRows(0), + fTPCNclsF(0), + fTPCNclsIter1(0), + fTPCClusterMap(0), + fTPCsignal(0), + fTPCnSig(), + fTOFbeta(0.0), + fTOFnSig(), + fTRDntracklets(), + fTRDpid(), + fCaloClusterId(-999), + fBayesPID(), + fFlags(0) +{ + // + // Constructor + // + fDCA[0] = 0.0; fDCA[1]=0.0; + for(Int_t i=0; i<4; ++i) {fTPCnSig[i]=-999.; fTOFnSig[i]=-999.;} + for(Int_t i=0; i<3; ++i) {fBayesPID[i]=-999.;} + fTRDpid[0]=-999.; fTRDpid[1]=-999.; +} + + +//_______________________________________________________________________________ +AliCorrelationReducedTrack::~AliCorrelationReducedTrack() +{ + // + // De-Constructor + // +} + + +//_______________________________________________________________________________ +AliCorrelationReducedPair::AliCorrelationReducedPair() : + fCandidateId(-1), + fPairType(-1), + fLegIds(), + fMass(), + fPhi(0.0), + fPt(0.0), + fEta(0.0), + fLxy(0.0), + fOpeningAngle(-1.0), + //fOnTheFly(kFALSE), + fMCid(0) +{ + // + // Constructor + // + fLegIds[0] = -1; fLegIds[1] = -1; + fMass[0]=-999.; fMass[1]=-999.; fMass[2]=-999.; +} + + +//_______________________________________________________________________________ +AliCorrelationReducedPair::AliCorrelationReducedPair(const AliCorrelationReducedPair &c) : + TObject(c), + fCandidateId(c.CandidateId()), + fPairType(c.PairType()), + fLegIds(), + fMass(), + fPhi(c.Phi()), + fPt(c.Pt()), + fEta(c.Eta()), + fLxy(c.Lxy()), + fOpeningAngle(c.OpeningAngle()), + //fOnTheFly(c.IsOnTheFly()), + fMCid(c.MCid()) +{ + // + // copy constructor + // + fLegIds[0] = c.LegId(0); + fLegIds[1] = c.LegId(1); + fMass[0] = c.Mass(0); fMass[1] = c.Mass(1); fMass[2] = c.Mass(2); +} + + +//_______________________________________________________________________________ +AliCorrelationReducedPair::~AliCorrelationReducedPair() { + // + // destructor + // +} + +//____________________________________________________________________________ +AliCorrelationReducedEvent::AliCorrelationReducedEvent() : + fRunNo(0), + fBC(0), + fTriggerMask(0), + fIsPhysicsSelection(kTRUE), + fVtx(), + fNVtxContributors(0), + fVtxTPC(), + fNVtxTPCContributors(0), + fCentrality(), + fCentQuality(0), + fNV0candidates(), + fNDielectronCandidates(0), + fNtracks(), + fSPDntracklets(0), + fVZEROMult(), + fZDCnEnergy(), + fNCaloClusters(0) +//fFMDMult() +{ + // + // Constructor + // + for(Int_t i=0; i<3; ++i) {fVtx[i]=-999.; fVtxTPC[i]=-999.;} + for(Int_t i=0; i<4; ++i) fCentrality[i]=-1.; + fNV0candidates[0]=0; fNV0candidates[1]=0; + fNtracks[0]=0; fNtracks[1]=0; + for(Int_t i=0; i<64; ++i) fVZEROMult[i] = 0.0; + for(Int_t i=0; i<8; ++i) fZDCnEnergy[i]=0.0; + + if(!fgTracks) fgTracks = new TClonesArray("AliCorrelationReducedTrack", 100000); + fTracks = fgTracks; + if(!fgCandidates) fgCandidates = new TClonesArray("AliCorrelationReducedPair", 100000); + fCandidates = fgCandidates; + if(!fgCaloClusters) fgCaloClusters = new TClonesArray("AliCorrelationReducedCaloCluster", 50000); + fCaloClusters = fgCaloClusters; +} + + +//____________________________________________________________________________ +AliCorrelationReducedEvent::~AliCorrelationReducedEvent() +{ + // + // De-Constructor + // + ClearEvent(); +} + + +//____________________________________________________________________________ +Float_t AliCorrelationReducedEvent::MultVZEROA() const +{ + // + // Total VZERO multiplicity in A side + // + Float_t mult=0.0; + for(Int_t i=32;i<64;++i) + mult += fVZEROMult[i]; + return mult; +} + + +//____________________________________________________________________________ +Float_t AliCorrelationReducedEvent::MultVZEROC() const +{ + // + // Total VZERO multiplicity in C side + // + Float_t mult=0.0; + for(Int_t i=0;i<32;++i) + mult += fVZEROMult[i]; + return mult; +} + + +//____________________________________________________________________________ +Float_t AliCorrelationReducedEvent::MultVZERO() const +{ + // + // Total VZERO multiplicity + // + return MultVZEROA()+MultVZEROC(); +} + + +//____________________________________________________________________________ +Float_t AliCorrelationReducedEvent::MultRingVZEROA(Int_t ring) const +{ + // + // VZERO multiplicity in a given ring on A side + // + if(ring<0 || ring>3) return -1.0; + + Float_t mult=0.0; + for(Int_t i=32+8*ring; i<32+8*(ring+1); ++i) + mult += fVZEROMult[i]; + return mult; +} + + +//____________________________________________________________________________ +Float_t AliCorrelationReducedEvent::MultRingVZEROC(Int_t ring) const +{ + // + // VZERO multiplicity in a given ring on C side + // + if(ring<0 || ring>3) return -1.0; + + Float_t mult=0.0; + for(Int_t i=8*ring; i<8*(ring+1); ++i) + mult += fVZEROMult[i]; + return mult; +} + +//_____________________________________________________________________________ +void AliCorrelationReducedEvent::ClearEvent() { + // + // clear the event + // + fTracks->Clear("C"); + fCandidates->Clear("C"); + fCaloClusters->Clear("C"); + fRunNo = 0; + fBC = 0; + fTriggerMask = 0; + fIsPhysicsSelection = kTRUE; + fCentQuality = 0; + fNV0candidates[0] = 0; fNV0candidates[1] = 0; + fNDielectronCandidates = 0; + fNtracks[0] = 0; fNtracks[1] = 0; + for(Int_t i=0; i<3; ++i) {fVtx[i]=-999.; fVtxTPC[i]=-999.; fCentrality[i]=-1.;} + for(Int_t i=0; i<64; ++i) fVZEROMult[i] = 0.0; + for(Int_t i=0; i<8; ++i) fZDCnEnergy[i]=0.0; +} + + +//_______________________________________________________________________________ +AliCorrelationReducedEventFriend::AliCorrelationReducedEventFriend() : + fQvector(), + fEventPlaneStatus() +{ + // + // default constructor + // + for(Int_t idet=0; idetQx(idet, ih+1); + fQvector[idet][ih][1] = event->Qy(idet, ih+1); + fEventPlaneStatus[idet][ih] = event->GetEventPlaneStatus(idet, ih+1); + } + } +} + + +//_____________________________________________________________________________ +AliCorrelationReducedCaloCluster::AliCorrelationReducedCaloCluster() : + fType(kUndefined), + fEnergy(-999.), + fTrackDx(-999.), + fTrackDz(-999.) +{ + // + // default constructor + // +} + + +//_____________________________________________________________________________ +AliCorrelationReducedCaloCluster::~AliCorrelationReducedCaloCluster() +{ + // + // destructor + // +} + + +//_______________________________________________________________________________ +void AliCorrelationReducedEvent::GetQvector(Double_t Qvec[][2], Int_t det) { + // + // Get the event plane for a specified detector + // + if(det==AliCorrelationReducedEventFriend::kTPC || + det==AliCorrelationReducedEventFriend::kTPCptWeights || + det==AliCorrelationReducedEventFriend::kTPCpos || + det==AliCorrelationReducedEventFriend::kTPCneg) { + GetTPCQvector(Qvec, det); + return; + } + if(det==AliCorrelationReducedEventFriend::kVZEROA || + det==AliCorrelationReducedEventFriend::kVZEROC) { + GetVZEROQvector(Qvec, det); + return; + } + if(det==AliCorrelationReducedEventFriend::kZDCA || + det==AliCorrelationReducedEventFriend::kZDCC) { + GetZDCQvector(Qvec, det); + return; + } + if(det==AliCorrelationReducedEventFriend::kFMD) { + //TODO implementation + return; + } + return; +} + + +//_________________________________________________________________ +void AliCorrelationReducedEvent::GetTPCQvector(Double_t Qvec[][2], Int_t det) { + // + // Construct the event plane using tracks in the barrel + // + if(!(det==AliCorrelationReducedEventFriend::kTPC || + det==AliCorrelationReducedEventFriend::kTPCpos || + det==AliCorrelationReducedEventFriend::kTPCneg)) + return; + AliCorrelationReducedTrack* track=0x0; + Double_t pt=0.0; Double_t x=0.0; Double_t y=0.0; + Short_t charge=0; + TIter nextTrack(fTracks); + while((track=static_cast(nextTrack()))) { + if(!track->UsedForQvector()) continue; + charge = track->Charge(); + if(det==AliCorrelationReducedEventFriend::kTPCpos && charge<0) continue; + if(det==AliCorrelationReducedEventFriend::kTPCneg && charge>0) continue; + pt = 1.0; + if(det==AliCorrelationReducedEventFriend::kTPCptWeights) { + pt = track->Pt(); + if(pt>2.0) pt = 2.0; // pt is the weight used for the event plane + } + x = TMath::Cos(track->Phi()); + y = TMath::Sin(track->Phi()); + // 1st harmonic + Qvec[0][0] += pt*x; + Qvec[0][1] += pt*y; + // 2nd harmonic + Qvec[1][0] += pt*(2.0*TMath::Power(x,2.0)-1); + Qvec[1][1] += pt*(2.0*x*y); + // 3rd harmonic + Qvec[2][0] += pt*(4.0*TMath::Power(x,3.0)-3.0*x); + Qvec[2][1] += pt*(3.0*y-4.0*TMath::Power(y,3.0)); + // 4th harmonic + Qvec[3][0] += pt*(1.0-8.0*TMath::Power(x*y,2.0)); + Qvec[3][1] += pt*(4.0*x*y-8.0*x*TMath::Power(y,3.0)); + // 5th harmonic + Qvec[4][0] += pt*(16.0*TMath::Power(x,5.0)-20.0*TMath::Power(x, 3.0)+5.0*x); + Qvec[4][1] += pt*(16.0*TMath::Power(y,5.0)-20.0*TMath::Power(y, 3.0)+5.0*y); + // 6th harmonic + Qvec[5][0] += pt*(32.0*TMath::Power(x,6.0)-48.0*TMath::Power(x, 4.0)+18.0*TMath::Power(x, 2.0)-1.0); + Qvec[5][1] += pt*(x*y*(32.0*TMath::Power(y,4.0)-32.0*TMath::Power(y, 2.0)+6.0)); + } +} + + +//_________________________________________________________________ +void AliCorrelationReducedEvent::GetVZEROQvector(Double_t Qvec[][2], Int_t det) { + // + // Get the reaction plane from the VZERO detector for a given harmonic + // + GetVZEROQvector(Qvec, det, fVZEROMult); +} + + +//_________________________________________________________________ +void AliCorrelationReducedEvent::GetVZEROQvector(Double_t Qvec[][2], Int_t det, Float_t* vzeroMult) { + // + // Get the reaction plane from the VZERO detector for a given harmonic + // + // Q{x,y} = SUM_i mult(i) * {cos(n*phi_i), sin(n*phi_i)} + // phi_i - phi angle of the VZERO sector i + // Each sector covers 45 degrees(8 sectors per ring). Middle of sector 0 is at 45/2 + // channel 0: 22.5 + // 1: 22.5+45 + // 2: 22.5+45*2 + // ... + // at the next ring continues the same + // channel 8: 22.5 + // channel 9: 22.5 + 45 + // + if(!(det==AliCorrelationReducedEventFriend::kVZEROA || + det==AliCorrelationReducedEventFriend::kVZEROC)) + return; + + const Double_t kX[8] = {0.92388, 0.38268, -0.38268, -0.92388, -0.92388, -0.38268, 0.38268, 0.92388}; // cosines of the angles of the VZERO sectors (8 per ring) + const Double_t kY[8] = {0.38268, 0.92388, 0.92388, 0.38268, -0.38268, -0.92388, -0.92388, -0.38268}; // sines -- " -- + Int_t phi; + + for(Int_t iChannel=0; iChannel<64; ++iChannel) { + if(iChannel<32 && det==AliCorrelationReducedEventFriend::kVZEROA) continue; + if(iChannel>=32 && det==AliCorrelationReducedEventFriend::kVZEROC) continue; + phi=iChannel%8; + // 1st harmonic + Qvec[0][0] += vzeroMult[iChannel]*kX[phi]; + Qvec[0][1] += vzeroMult[iChannel]*kY[phi]; + // 2nd harmonic + Qvec[1][0] += vzeroMult[iChannel]*(2.0*TMath::Power(kX[phi],2.0)-1); + Qvec[1][1] += vzeroMult[iChannel]*(2.0*kX[phi]*kY[phi]); + // 3rd harmonic + Qvec[2][0] += vzeroMult[iChannel]*(4.0*TMath::Power(kX[phi],3.0)-3.0*kX[phi]); + Qvec[2][1] += vzeroMult[iChannel]*(3.0*kY[phi]-4.0*TMath::Power(kY[phi],3.0)); + // 4th harmonic + Qvec[3][0] += vzeroMult[iChannel]*(1.0-8.0*TMath::Power(kX[phi]*kY[phi],2.0)); + Qvec[3][1] += vzeroMult[iChannel]*(4.0*kX[phi]*kY[phi]-8.0*kX[phi]*TMath::Power(kY[phi],3.0)); + // 5th harmonic + Qvec[4][0] += vzeroMult[iChannel]*(16.0*TMath::Power(kX[phi],5.0)-20.0*TMath::Power(kX[phi], 3.0)+5.0*kX[phi]); + Qvec[4][1] += vzeroMult[iChannel]*(16.0*TMath::Power(kY[phi],5.0)-20.0*TMath::Power(kY[phi], 3.0)+5.0*kY[phi]); + // 6th harmonic + Qvec[5][0] += vzeroMult[iChannel]*(32.0*TMath::Power(kX[phi],6.0)-48.0*TMath::Power(kX[phi], 4.0)+18.0*TMath::Power(kX[phi], 2.0)-1.0); + Qvec[5][1] += vzeroMult[iChannel]*(kX[phi]*kY[phi]*(32.0*TMath::Power(kY[phi],4.0)-32.0*TMath::Power(kY[phi], 2.0)+6.0)); + } // end loop over channels +} + + +//_________________________________________________________________ +void AliCorrelationReducedEvent::GetZDCQvector(Double_t Qvec[][2], Int_t det) { + // + // Construct the event plane using the ZDC + // ZDC has 2 side (A and C) with 4 calorimeters on each side + // The XY position of each calorimeter is specified by the + // zdcNTowerCenters_x and zdcNTowerCenters_y arrays + if(!(det==AliCorrelationReducedEventFriend::kZDCA || + det==AliCorrelationReducedEventFriend::kZDCC)) return; // bad detector option + const Float_t zdcTowerCenter = 1.75; + const Float_t zdcNTowerCenters_x[4] = {-zdcTowerCenter, zdcTowerCenter, -zdcTowerCenter, zdcTowerCenter}; + const Float_t zdcNTowerCenters_y[4] = {-zdcTowerCenter, -zdcTowerCenter, zdcTowerCenter, zdcTowerCenter}; + + Qvec[0][0] = 0.0; Qvec[0][1] = 0.0; // first harmonic Q-vector + Float_t zdcNCentroidSum = 0; + Float_t zdcN_alpha = 0.5; + + for(Int_t i=0; i<4; ++i) { + if(fZDCnEnergy[i+(det==AliCorrelationReducedEventFriend::kZDCA ? 4 : 0)]>0.0) { + Float_t zdcNenergy_alpha = TMath::Power(fZDCnEnergy[i+(det==AliCorrelationReducedEventFriend::kZDCA ? 4 : 0)], zdcN_alpha); + Qvec[0][0] += zdcNTowerCenters_x[i-1]*zdcNenergy_alpha; + Qvec[0][1] += zdcNTowerCenters_y[i-1]*zdcNenergy_alpha; + zdcNCentroidSum += zdcNenergy_alpha; + } + } // end loop over channels + + if(zdcNCentroidSum>0.0) { + Qvec[0][0] /= zdcNCentroidSum; + Qvec[0][1] /= zdcNCentroidSum; + } +} diff --git a/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliCorrelationReducedEvent.h b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliCorrelationReducedEvent.h new file mode 100644 index 00000000000..bae25178eeb --- /dev/null +++ b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/AliCorrelationReducedEvent.h @@ -0,0 +1,458 @@ +// Classes used for creating a reduced information tree +// Author: Ionut-Cristian Arsene (i.c.arsene@gsi.de) +// +// Basic structure: +// 1. Event wise information +// 2. List of tracks in the event +// 3. List of resonance candidates + +#ifndef ALICORRELATIONREDUCEDEVENT_H +#define ALICORRELATIONREDUCEDEVENT_H + +#include +#include +#include + + +const Int_t fgkNMaxHarmonics = 10; +/*class AliCorrelationReducedTrack; +class AliCorrelationReducedPair; +class AliCorrelationReducedEventFriend; +class AliCorrelationReducedEvent; +class AliCorrelationReducedCaloCluster;*/ + +//_____________________________________________________________________ +class AliCorrelationReducedTrack : public TObject { + + friend class AliAnalysisTaskCorrelationTree; + + public: + AliCorrelationReducedTrack(); + ~AliCorrelationReducedTrack(); + + // getters + UShort_t TrackId() const {return fTrackId;} + ULong_t Status() const {return fStatus;} + Bool_t CheckTrackStatus(UInt_t flag) const {return (flag<8*sizeof(ULong_t) ? (fStatus&(1<0.0 ? +1 : -1);} + Float_t Px() const {return TMath::Abs(fGlobalPt)*TMath::Cos(fGlobalPhi);} + Float_t Py() const {return TMath::Abs(fGlobalPt)*TMath::Sin(fGlobalPhi);} + Float_t Pz() const {return TMath::Abs(fGlobalPt)*TMath::SinH(fGlobalEta);} + Float_t P() const {return TMath::Abs(fGlobalPt)*TMath::CosH(fGlobalEta);}; + Float_t Phi() const {return fGlobalPhi;} + Float_t Pt() const {return TMath::Abs(fGlobalPt);} + Float_t Eta() const {return fGlobalEta;} + Float_t Theta() const {return TMath::ACos(TMath::TanH(fGlobalEta));} + Float_t PxTPC() const {return fTPCPt*TMath::Cos(fTPCPhi);} + Float_t PyTPC() const {return fTPCPt*TMath::Sin(fTPCPhi);} + Float_t PzTPC() const {return fTPCPt*TMath::SinH(fTPCEta);} + Float_t PTPC() const {return fTPCPt*TMath::CosH(fTPCEta);}; + Float_t PhiTPC() const {return fTPCPhi;} + Float_t PtTPC() const {return fTPCPt;} + Float_t EtaTPC() const {return fTPCEta;} + Float_t ThetaTPC() const {return TMath::ACos(TMath::TanH(fTPCEta));} + Float_t Pin() const {return fMomentumInner;} + Float_t DCAxy() const {return fDCA[0];} + Float_t DCAz() const {return fDCA[1];} + + UShort_t ITSncls() const; + UChar_t ITSclusterMap() const {return fITSclusterMap;} + Bool_t ITSLayerHit(Int_t layer) const {return (layer>=0 && layer<6 ? (fITSclusterMap&(1<=0 && bit<8 ? (fTPCClusterMap&(1<=0 && specie<=3 ? fTPCnSig[specie] : -999.);} + + Float_t TOFbeta() const {return fTOFbeta;} + Float_t TOFnSig(Int_t specie) const {return (specie>=0 && specie<=3 ? fTOFnSig[specie] : -999.);} + + Int_t TRDntracklets(Int_t type) const {return (type==0 || type==1 ? fTRDntracklets[type] : -1);} + Float_t TRDpid(Int_t specie) const {return (specie>=0 && specie<=1 ? fTRDpid[specie] : -999.);} + + Int_t CaloClusterId() const {return fCaloClusterId;} + //Float_t CaloEnergy(AliCorrelationReducedEvent* event) const {if(fCaloClusterId>0) return event->GetCaloCluster(fCaloClusterId)->Energy();} + //Float_t CaloDx(AliCorrelationReducedEvent* event) const {if(fCaloClusterId>0) return event->GetCaloCluster(fCaloClusterId)->Dx();} + //Float_t CaloDz(AliCorrelationReducedEvent* event) const {if(fCaloClusterId>0) return event->GetCaloCluster(fCaloClusterId)->Dz();} + + Float_t BayesPID(Int_t specie) const {return (specie>=0 && specie<=2 ? fBayesPID[specie] : -999.);} + Bool_t UsedForQvector() const {return fFlags&(1<<0);} + + private: + UShort_t fTrackId; // track id + ULong_t fStatus; // tracking status + Float_t fGlobalPhi; // phi at the vertex from global track, in the [0,2pi) interval + Float_t fGlobalPt; // pt*charge at the vertex from global track + Float_t fGlobalEta; // eta at the vertex from global track + Float_t fTPCPhi; // phi at the vertex from TPC alone tracking , in the [0,2pi) interval + Float_t fTPCPt; // pt at the vertex from TPC alone tracking + Float_t fTPCEta; // eta at the vertex from TPC alone tracking + Float_t fMomentumInner; // inner param momentum (only the magnitude) + Float_t fDCA[2]; // DCA xy,z + + // ITS + UChar_t fITSclusterMap; // ITS cluster map + Float_t fITSsignal; // ITS signal + + // TPC + UChar_t fTPCNcls; // TPC ncls + UChar_t fTPCCrossedRows; // TPC crossed rows + UChar_t fTPCNclsF; // TPC findable ncls + UChar_t fTPCNclsIter1; // TPC no clusters after first iteration + UChar_t fTPCClusterMap; // TPC cluster distribution map + Float_t fTPCsignal; // TPC de/dx + Float_t fTPCnSig[4]; // 0-electron; 1-pion; 2-kaon; 3-proton + + // TOF + Float_t fTOFbeta; // TOF pid info + Float_t fTOFnSig[4]; // TOF n-sigma deviation from expected signal + + // TRD + UChar_t fTRDntracklets[2]; // 0 - AliESDtrack::GetTRDntracklets(); 1 - AliESDtrack::GetTRDntrackletsPID() TODO: use only 1 char + Float_t fTRDpid[2]; // TRD pid probabilities, [0]-electron, [1]-pion + + // EMCAL/PHOS + Int_t fCaloClusterId; // ID for the calorimeter cluster (if any) + + // Bayesian PID + Float_t fBayesPID[3]; // Combined Bayesian PID pi/K/p + + UShort_t fFlags; // BIT0 toggled if track used for TPC event plane TODO combine with other posible flags, use for MC pid? + // TODO flag for which TPC part used for pid --> Char_t used in 2011 data + + AliCorrelationReducedTrack(const AliCorrelationReducedTrack &c); + AliCorrelationReducedTrack& operator= (const AliCorrelationReducedTrack &c); + + ClassDef(AliCorrelationReducedTrack, 2); +}; + + +//_____________________________________________________________________ +class AliCorrelationReducedPair : public TObject { + + friend class AliAnalysisTaskCorrelationTree; + + public: + enum CandidateType { + kK0sToPiPi=0, + kPhiToKK, + kLambda0ToPPi, + kALambda0ToPPi, + kJpsiToEE, + kNMaxCandidateTypes + }; + AliCorrelationReducedPair(); + AliCorrelationReducedPair(const AliCorrelationReducedPair &c); + ~AliCorrelationReducedPair(); + + // getters + Char_t CandidateId() const {return fCandidateId;} + Char_t PairType() const {return fPairType;} + Int_t LegId(Int_t leg) const {return (leg==0 || leg==1 ? fLegIds[leg] : -1);} + Float_t Mass(Int_t idx=0) const {return (idx>=0 && idx<3 ? fMass[idx] : -999.);} + Float_t Px() const {return fPt*TMath::Cos(fPhi);} + Float_t Py() const {return fPt*TMath::Sin(fPhi);} + Float_t Pz() const {return fPt*TMath::SinH(fEta);} + Float_t P() const {return fPt*TMath::CosH(fEta);} + Float_t Phi() const {return fPhi;} + Float_t Pt() const {return fPt;} + Float_t Eta() const {return fEta;} + Float_t Energy() const; + Float_t Rapidity() const; + Float_t Theta() const {return TMath::ACos(TMath::TanH(fEta));} + Float_t Lxy() const {return fLxy;} + Float_t OpeningAngle() const {return fOpeningAngle;} + Bool_t IsOnTheFly() const {return fPairType;} + UInt_t MCid() const {return fMCid;} + Bool_t CheckMC(const Int_t flag) const {return (flag<32 ? (fMCid&(1< K0s assumption; idx=1 -> Lambda; idx=2 -> anti-Lambda + Float_t fPhi; // pair phi in the [0,2*pi) interval + Float_t fPt; // pair pt + Float_t fEta; // pair eta + Float_t fLxy; // pseudo-proper decay length + Float_t fOpeningAngle; // opening angle TODO remove ??? + UInt_t fMCid; // Bit map with Monte Carlo info about the pair + + AliCorrelationReducedPair& operator= (const AliCorrelationReducedPair &c); + + ClassDef(AliCorrelationReducedPair, 1); +}; + + +//_________________________________________________________________________ +class AliCorrelationReducedEventFriend : public TObject { + + friend class AliAnalysisTaskCorrelationTree; + + public: + enum EventPlaneStatus { + kRaw=0, + kCalibrated, + kRecentered, + kShifted, + kNMaxFlowFlags + }; + enum EventPlaneDetector { + kTPC=0, + kTPCptWeights, + kTPCpos, + kTPCneg, + kVZEROA, + kVZEROC, + kFMD, + kZDCA, + kZDCC, + kNdetectors + }; + + AliCorrelationReducedEventFriend(); + ~AliCorrelationReducedEventFriend(); + + Double_t Qx(Int_t det, Int_t harmonic) const {return (det>=0 && det0 && harmonic<=fgkNMaxHarmonics ? fQvector[det][harmonic-1][0] : -999.);} + Double_t Qy(Int_t det, Int_t harmonic) const {return (det>=0 && det0 && harmonic<=fgkNMaxHarmonics ? fQvector[det][harmonic-1][1] : -999.);} + Double_t EventPlane(Int_t det, Int_t h) const; + UChar_t GetEventPlaneStatus(Int_t det, Int_t h) const {return (det>=0 && det0 && h<=fgkNMaxHarmonics ? fEventPlaneStatus[det][h] : -999.);} + Bool_t CheckEventPlaneStatus(Int_t det, Int_t h, EventPlaneStatus flag) const; + void CopyEvent(AliCorrelationReducedEventFriend* event); + + void SetQx(Int_t det, Int_t harmonic, Float_t qx) { if(det>=0 && det0 && harmonic<=fgkNMaxHarmonics) fQvector[det][harmonic-1][0]=qx;} + void SetQy(Int_t det, Int_t harmonic, Float_t qy) { if(det>=0 && det0 && harmonic<=fgkNMaxHarmonics) fQvector[det][harmonic-1][1]=qy;} + void SetEventPlaneStatus(Int_t det, Int_t harmonic, EventPlaneStatus status) { + if(det>=0 && det0 && harmonic<=fgkNMaxHarmonics) + fEventPlaneStatus[det][harmonic-1] |= (1<=0 && axis<=2 ? fVtx[axis] : 0);} + Int_t VertexNContributors() const {return fNVtxContributors;} + Float_t VertexTPC(Int_t axis) const {return (axis>=0 && axis<=2 ? fVtxTPC[axis] : 0);} + Int_t VertexTPCContributors() const {return fNVtxTPCContributors;} + Float_t CentralityVZERO() const {return fCentrality[0];} + Float_t CentralitySPD() const {return fCentrality[1];} + Float_t CentralityTPC() const {return fCentrality[2];} + Float_t CentralityZEMvsZDC() const {return fCentrality[3];} + Int_t CentralityQuality() const {return fCentQuality;} + Int_t NV0CandidatesTotal() const {return fNV0candidates[0];} + Int_t NV0Candidates() const {return fNV0candidates[1];} + Int_t NDielectrons() const {return fNDielectronCandidates;} + Int_t NTracksTotal() const {return fNtracks[0];} + Int_t NTracks() const {return fNtracks[1];} + Int_t SPDntracklets() const {return fSPDntracklets;} + + Float_t MultChannelVZERO(Int_t channel) const {return (channel>=0 && channel<=63 ? fVZEROMult[channel] : -999.);} + Float_t MultVZEROA() const; + Float_t MultVZEROC() const; + Float_t MultVZERO() const; + Float_t MultRingVZEROA(Int_t ring) const; + Float_t MultRingVZEROC(Int_t ring) const; + + Float_t EnergyZDC(Int_t channel) const {return (channel>=0 && channel<8 ? fZDCnEnergy[channel] : -999.);} + Float_t EnergyZDCnA(Int_t channel) const {return (channel>=0 && channel<4 ? fZDCnEnergy[channel+4] : -999.);} + Float_t EnergyZDCnC(Int_t channel) const {return (channel>=0 && channel<4 ? fZDCnEnergy[channel] : -999.);} + + AliCorrelationReducedTrack* GetTrack(Int_t i) const {return (iAt(i) : 0x0);} + AliCorrelationReducedPair* GetV0Pair(Int_t i) const {return (i>=0 && iAt(i) : 0x0);} + AliCorrelationReducedPair* GetDielectronPair(Int_t i) const {return (i>=0 && iAt(i+fNV0candidates[1]) : 0x0);} + TClonesArray* GetPairs() const {return fCandidates;} + TClonesArray* GetTracks() const {return fTracks;} + + Int_t GetNCaloClusters() const {return fNCaloClusters;} + AliCorrelationReducedCaloCluster* GetCaloCluster(Int_t i) const {return (i>=0 && iAt(i) : 0x0);} + + void GetQvector(Double_t Qvec[][2], Int_t det); + void GetTPCQvector(Double_t Qvec[][2], Int_t det); + void GetVZEROQvector(Double_t Qvec[][2], Int_t det); + void GetVZEROQvector(Double_t Qvec[][2], Int_t det, Float_t* vzeroMult); + void GetZDCQvector(Double_t Qvec[][2], Int_t det); + + private: + Int_t fRunNo; // run number + UShort_t fBC; // bunch crossing + ULong64_t fTriggerMask; // trigger mask + Bool_t fIsPhysicsSelection; // PhysicsSelection passed event + Float_t fVtx[3]; // global event vertex vector in cm + Int_t fNVtxContributors; // global event vertex contributors + Float_t fVtxTPC[3]; // TPC only event vertex + Int_t fNVtxTPCContributors; // TPC only event vertex contributors + Float_t fCentrality[4]; // centrality; 0-VZERO, 1-SPD, 2-TPC, 3-ZEMvsZDC + Int_t fCentQuality; // quality flag for the centrality + Int_t fNV0candidates[2]; // number of V0 candidates, [0]-total, [1]-selected for the tree + Int_t fNDielectronCandidates; // number of pairs selected as dielectrons + Int_t fNtracks[2]; // number of tracks, [0]-total, [1]-selected for the tree + Int_t fSPDntracklets; // number of SPD tracklets in |eta|<1.0 + + Float_t fVZEROMult[64]; // VZERO multiplicity in all 64 channels + Float_t fZDCnEnergy[8]; // neutron ZDC energy in all 8 channels + + TClonesArray* fTracks; //-> array containing global tracks + static TClonesArray* fgTracks; + + TClonesArray* fCandidates; //-> array containing pair candidates + static TClonesArray* fgCandidates; + + Int_t fNCaloClusters; // number of calorimeter clusters + TClonesArray* fCaloClusters; //-> array containing calorimeter clusters + static TClonesArray* fgCaloClusters; + + void ClearEvent(); + AliCorrelationReducedEvent(const AliCorrelationReducedEvent &c); + AliCorrelationReducedEvent& operator= (const AliCorrelationReducedEvent &c); + + ClassDef(AliCorrelationReducedEvent, 2); +}; + +//_______________________________________________________________________________ +inline UShort_t AliCorrelationReducedTrack::ITSncls() const +{ + // + // ITS number of clusters from the cluster map + // + UShort_t ncls=0; + for(Int_t i=0; i<6; ++i) ncls += (ITSLayerHit(i) ? 1 : 0); + return ncls; +} + + +//_______________________________________________________________________________ +inline Int_t AliCorrelationReducedTrack::TPCClusterMapBitsFired() const +{ + // + // Count the number of bits fired in the TPC cluster map + // + Int_t nbits=0; + for(Int_t i=0; i<8; ++i) nbits += (TPCClusterMapBitFired(i) ? 1 : 0); + return nbits; +} + + +//_______________________________________________________________________________ +inline Float_t AliCorrelationReducedPair::Energy() const +{ + // + // Return the energy + // + Float_t mass=fMass[0]; + switch (fCandidateId) { + case kK0sToPiPi: + mass = fMass[0]; + break; + case kLambda0ToPPi: + mass = fMass[1]; + break; + case kALambda0ToPPi: + mass = fMass[2]; + break; + default: + mass = fMass[0]; + break; + } + Float_t p = P(); + return TMath::Sqrt(mass*mass+p*p); +} + + +//_______________________________________________________________________________ +inline Float_t AliCorrelationReducedPair::Rapidity() const +{ + // + // return rapidity + // + Float_t e = Energy(); + Float_t pz = Pz(); + if(e-TMath::Abs(pz)>1.0e-10) + return 0.5*TMath::Log((e+pz)/(e-pz)); + else + return -999.; +} + + +//_______________________________________________________________________________ +inline Double_t AliCorrelationReducedEventFriend::EventPlane(Int_t det, Int_t harmonic) const +{ + // + // Event plane from detector "det" and harmonic "harmonic" + // + if(det<0 || det>=kNdetectors || harmonic<1 || harmonic>fgkNMaxHarmonics) return -999.; + return TMath::ATan2(fQvector[det][harmonic-1][1], fQvector[det][harmonic-1][0])/Double_t(harmonic); +} + +//_______________________________________________________________________________ +inline Bool_t AliCorrelationReducedEventFriend::CheckEventPlaneStatus(Int_t det, Int_t h, EventPlaneStatus flag) const { + // + // Check the status of the event plane for a given detector and harmonic + // + if(det<0 || det>=kNdetectors || h<1 || h>fgkNMaxHarmonics) return kFALSE; + return (flag +#include +using namespace std; + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#ifndef ALICORRELATIONREDUCEDEVENT_H +#include "AliCorrelationReducedEvent.h" +#endif + +namespace DstCommonMacros { + + enum ParticleId { + kUnknown = -1, + kElectron = 0, + kPion, + kKaon, + kProton + }; + + enum Variables { + kNothing = -1, + // Event wise variables + kRunNo = 0, // run number + kBC, // bunch crossing + kTriggerMask, // trigger mask + kOfflineTrigger, // offline trigger + kOfflineTriggerFired, // offline trigger fired + kOfflineTriggerFired2, // offline trigger if fired, -1 if not fired + kIsPhysicsSelection, // physics selection + kVtxX, // vtx X + kVtxY, // vtx Y + kVtxZ, // vtx Z + kVtxXtpc, // vtx X from tpc + kVtxYtpc, // vtx Y from tpc + kVtxZtpc, // vtx Z from tpc + kCentVZERO, // centrality from VZERO + kCentSPD, // centrality from SPD + kCentTPC, // centrality from TPC + kCentZDC, // centrality from ZDC + kCentQuality, // centrality quality + kNV0total, // total number of V0s in the esd + kNV0selected, // number of V0s selected + kNdielectrons, // number of dielectron pairs + kNpairsSelected, // number of selected dielectron pairs per event + kNtracksTotal, // total number of tracks + kNtracksSelected, // number of selected tracks + kSPDntracklets, // SPD number of tracklets in |eta|<1.0 + kEventMixingId, // Id of the event mixing category + // VZERO event plane related variables + kVZEROAemptyChannels, // Number of empty VZERO channels in A side + kVZEROCemptyChannels, // Number of empty VZERO channels in C side + kVZEROChannelMult, // VZERO multiplicity per channel + kVZEROChannelEta = kVZEROChannelMult+64, // pseudo-rapidity of a VZERO channel + kVZEROQvecX = kVZEROChannelEta+64, // Q-vector components for harmonics 1-6 and + kVZEROQvecY = kVZEROQvecX+6*3, // 6- n-harmonics; 3- A,C and A&C options + kVZERORP = kVZEROQvecY+6*3, // VZERO reaction plane from A,C and A&C sides (harmonics 1-6) + kVZERORPres = kVZERORP+6*3, // VZERO reaction plane resolution (sqrt(n*(RPa-RPc))) + kVZEROXaXc = kVZERORPres+6, // correlations for the components of the Q vector + kVZEROXaYa = kVZEROXaXc+6, + kVZEROXaYc = kVZEROXaYa+6, + kVZEROYaXc = kVZEROXaYc+6, + kVZEROYaYc = kVZEROYaXc+6, + kVZEROXcYc = kVZEROYaYc+6, + kVZEROdeltaRPac = kVZEROXcYc+6, // Psi_VZEROA-Psi_VZEROC + kVZEROdeltaRPa = kVZEROdeltaRPac+6, // Psi_VZEROA6-Psi_VZEROA5, Psi_VZEROA7-Psi_VZEROA5, Psi_VZEROA6-Psi_VZEROA5 + kVZEROdeltaRPc = kVZEROdeltaRPa+6*2, // Psi_VZEROA2-Psi_VZEROA1, Psi_VZEROA3-Psi_VZEROA1, Psi_VZEROA4-Psi_VZEROA1 + kVZEROflowV2TPC = kVZEROdeltaRPc+6*2, // vzero v2 using TPC event plane + // TPC event plane variables + kTPCQvecX = kVZEROflowV2TPC+64, // TPC Q-vector components for harmonics 1-6 + kTPCQvecY = kTPCQvecX+6, + kTPCRP = kTPCQvecY+6, // Event plane using TPC + kTPCRPres = kTPCRP+6, // Event plane resolution variables sqrt(n*(RPtpc-RPvzeroa)),sqrt(n*(RPtpc-RPvzeroc)) + // Correlations between TPC and VZERO event planes + kRPXtpcXvzeroa = kTPCRPres+6*2, + kRPXtpcXvzeroc = kRPXtpcXvzeroa+6, + kRPYtpcYvzeroa = kRPXtpcXvzeroc+6, + kRPYtpcYvzeroc = kRPYtpcYvzeroa+6, + kRPXtpcYvzeroa = kRPYtpcYvzeroc+6, + kRPXtpcYvzeroc = kRPXtpcYvzeroa+6, + kRPYtpcXvzeroa = kRPXtpcYvzeroc+6, + kRPYtpcXvzeroc = kRPYtpcXvzeroa+6, + kRPdeltaVZEROAtpc = kRPYtpcXvzeroc+6, + kRPdeltaVZEROCtpc = kRPdeltaVZEROAtpc+6, + kNEventVars = kRPdeltaVZEROCtpc+6, // number of event variables + // Pair variables -------------------------------------- + kMass=kNEventVars, + kCandidateId, + kPairType, // 0 ++; 1 +-; 2 -- + kPairPt, + kPairPx, + kPairPy, + kPairPz, + kPairP, + kPairRap, + kPairEta, + kPairTheta, + kPairPhi, + kPairLxy, + kPairOpeningAngle, + kPairCosNPhi, // cos (n*phi) + kPairSinNPhi = kPairCosNPhi+6, // sin (n*phi) + kPairDeltaPhiVZEROFlowVn = kPairSinNPhi+6, // phi - Psi_{VZERO} + kPairDeltaPhiTPCFlowVn = kPairDeltaPhiVZEROFlowVn+6*3, // phi - Psi_{TPC} + kPairVZEROFlowVn = kPairDeltaPhiTPCFlowVn+6, // vn{Psi_{n,VZERO}} + kPairTPCFlowVn = kPairVZEROFlowVn+6*3, // vn{Psi_{n,TPC}} + kPairVZEROFlowV3Psi2 = kPairTPCFlowVn+6, // v3{Psi_{2,VZERO}} + kPairTPCFlowV3Psi2 = kPairVZEROFlowV3Psi2+3, // v3{Psi_{2,TPC}} + // Track variables ------------------------------------- + kPt=kPairTPCFlowV3Psi2+1, + kP, + kTheta, + kPhi, + kEta, + kRap, + kPtTPC, + kPhiTPC, + kEtaTPC, + kPin, + kTrackDeltaPhiVZEROFlowVn, + kTrackVZEROFlowVn = kTrackDeltaPhiVZEROFlowVn+6*2, + kDcaXY = kTrackVZEROFlowVn+6*2, + kDcaZ, + kITSncls, + kITSsignal, + kTPCncls, + kTPCcrossedRows, + kTPCnclsIter1, + kTPCnclsF, + kTPCnclsRatio, + kTPCsignal, + kTPCnSig, + kTOFbeta=kTPCnSig+4, + kTOFnSig, + kTRDntracklets=kTOFnSig+4, + kTRDntrackletsPID, + kTRDpidProbabilities, + kEMCALmatchedEnergy=kTRDpidProbabilities+2, + kEMCALmatchedEOverP, + // Calorimeter cluster variables -------------------------------------- + kEMCALclusterEnergy, + kEMCALclusterDx, + kEMCALclusterDz, + kEMCALdetector, // 0 - EMCAL; 1 - PHOS + // Tracking flags ----------------------------------------------------- + kTrackingFlag, + // Correlation variables ---------------------------------------------- + kDeltaPhi, + kDeltaTheta, + kDeltaEta, + kNVars + }; + + + Bool_t gUsedVars[kNVars] = {kFALSE}; + + // tracking flags as in AliESDtrack.h + // NOTE: check consistency with aliroot + enum TrackingFlags { + kITSin=0, + kITSout, + kITSrefit, + kITSpid, + kTPCin, + kTPCout, + kTPCrefit, + kTPCpid, + kTRDin, + kTRDout, + kTRDrefit, + kTRDpid, + kTOFin, + kTOFout, + kTOFrefit, + kTOFpid, + kTOFmismatch, + kHMPIDout, + kHMPIDpid, + kEMCALmatch, + kPHOSmatch, + kTRDbackup, + kTRDStop, + kESDpid, + kTIME, + kGlobalMerge, + kITSpureSA, + kMultInV0, + kMultSec, + kTRDnPlanes, + kEMCALNoMatch, + kNTrackingFlags + }; + + const Char_t* gkTrackingFlagNames[kNTrackingFlags] = { + "kITSin", "kITSout", "kITSrefit", "kITSpid", + "kTPCin", "kTPCout", "kTPCrefit", "kTPCpid", + "kTRDin", "kTRDout", "kTRDrefit", "kTRDpid", + "kTOFin", "kTOFout", "kTOFrefit", "kTOFpid", "kTOFmismatch", + "kHMPIDout", "kHMPIDpid", + "kEMCALmatch", "kPHOSmatch", + "kTRDbackup", "kTRDStop", + "kESDpid", "kTIME", "kGlobalMerge", + "kITSpureSA", + "kMultInV0", + "kMultSec", + "kTRDnPlanes", + "kEMCALNoMatch" + }; + + // offline triggers as defined in AliVEvent.h + // NOTE: Check consistency with updates in aliroot!!! + enum EOfflineTriggerTypes { + kMB = BIT(0), // Minimum bias trigger, i.e. interaction trigger, offline SPD or V0 selection + kINT7 = BIT(1), // V0AND trigger, offline V0 selection + kMUON = BIT(2), // Muon trigger, offline SPD or V0 selection + kHighMult = BIT(3), // High-multiplicity trigger (threshold defined online), offline SPD or V0 selection + kEMC1 = BIT(4), // EMCAL trigger + kCINT5 = BIT(5), // Minimum bias trigger without SPD. i.e. interaction trigger, offline V0 selection + kCMUS5 = BIT(6), // Muon trigger, offline V0 selection + kMUSPB = BIT(6), // idem for PbPb + kMUSH7 = BIT(7), // Muon trigger: high pt, single muon, offline V0 selection, CINT7 suite + kMUSHPB = BIT(7), // idem for PbPb + kMUL7 = BIT(8), // Muon trigger: like sign dimuon, offline V0 selection, CINT7 suite + kMuonLikePB = BIT(8), // idem for PbPb + kMUU7 = BIT(9), // Muon trigger, unlike sign dimuon, offline V0 selection, CINT7 suite + kMuonUnlikePB = BIT(9), // idem for PbPb + kEMC7 = BIT(10), // EMCAL trigger, CINT7 suite + kMUS7 = BIT(11), // Muon trigger: low pt, single muon, offline V0 selection, CINT7 suite + kPHI1 = BIT(12), // PHOS trigger, CINT1 suite + kPHI7 = BIT(13), // PHOS trigger, CINT7 suite + kPHOSPb = BIT(13), // idem for PbPb + kEMCEJE = BIT(14), // EMCAL jet patch trigger + kEMCEGA = BIT(15), // EMCAL gamma trigger + kCentral = BIT(16), // PbPb central collision trigger + kSemiCentral = BIT(17), // PbPb semicentral collision trigger + kDG5 = BIT(18), // Double gap diffractive + kZED = BIT(19), // ZDC electromagnetic dissociation + kUserDefined = BIT(27), // Set when custom trigger classes are set in AliPhysicsSelection, offline SPD or V0 selection + // Bits 28 and above are reserved for FLAGS + kFastOnly = BIT(30), // The fast cluster fired. This bit is set in to addition another trigger bit, e.g. kMB + kAny = 0xffffffff, // to accept any trigger + kAnyINT = kMB | kINT7 | kCINT5 // to accept any interaction (aka minimum bias) trigger + }; + + const Char_t* gkOfflineTriggerNames[64] = { + "MB", "INT7", "MUON", "HighMult", "EMC1", "CINT5", "CMUS5/MUSPB", "MUSH7/MUSHPB", + "MUL7/MuonLikePB", "MUU7/MuonUnlikePB", "EMC7", "MUS7", "PHI1", "PHI7/PHOSPb", "EMCEJE", "EMCEGA", + "Central", "SemiCentral", "DG5", "ZED", "N/A", "N/A", "N/A", "N/A", + "N/A", "N/A", "N/A", "UserDefined", "N/A", "N/A", "FastOnly", "N/A", + "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", + "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", + "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", + "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A", "N/A" + }; + + enum ITSLayerMap { + kITSfirst = 1, + kITSsecond = 2, + kITSthird = 4, + kITSfourth = 8, + kITSfifth = 16, + kITSsixth = 32 + }; + + // radii of VZERO channels centers (in cm) + const Double_t gkVZEROChannelRadii[64] = {6.0567, 6.0567, 6.0567, 6.0567, 6.0567, 6.0567, 6.0567, 6.0567, + 9.6977, 9.6977, 9.6977, 9.6977, 9.6977, 9.6977, 9.6977, 9.6977, + 15.9504, 15.9504, 15.9504, 15.9504, 15.9504, 15.9504, 15.9504, 15.9504, + 26.4031, 26.4031, 26.4031, 26.4031, 26.4031, 26.4031, 26.4031, 26.4031, + 5.9347, 5.9347, 5.9347, 5.9347, 5.9347, 5.9347, 5.9347, 5.9347, + 10.685, 10.685, 10.685, 10.685, 10.685, 10.685, 10.685, 10.685, + 18.116, 18.116, 18.116, 18.116, 18.116, 18.116, 18.116, 18.116, + 31.84, 31.84, 31.84, 31.84, 31.84, 31.84, 31.84, 31.84}; + const Double_t gkVZEROAz = 340.0; // cm + const Double_t gkVZEROCz = 90.0; // cm + const Double_t gkVZEROminMult = 0.5; // minimum VZERO channel multiplicity + + // Pointer to the current event + AliCorrelationReducedEvent* gCurrentEvent = 0x0; + AliCorrelationReducedEventFriend* gCurrentEventFriend = 0x0; + + // Event mixing variables + Int_t gEMCategories=0; + TString* gEMCategoryNames; + + TObjArray* gHistLists=0x0; // main histogram list for the current running process + TDirectoryFile* gHistListsOld=0x0; // main directory for a standard tree analysis output (used for calibration, plotting etc.) + TFile* gHistFile = 0x0; // pointer to a TFile opened for reading + + TProfile2D* gVzeroAvMult[64] = {0x0}; // pointer to the array of average VZERO multiplicity per channel + TProfile2D* gQvecCentering[AliCorrelationReducedEventFriend::kNdetectors][fgkNMaxHarmonics][2] = {{{0x0}}}; // pointer to the array of Qvec centering histograms + + // pt range for J/psi's + const Float_t gkJpsiPtCut[2] = {0.0, 20.0}; + + // Function prototypes + void WriteOutput(TFile* saveFile); + //void DefineHistograms(const Char_t* histClasses); + TChain* GetChain(const Char_t* filename, Int_t howMany, Int_t offset, Long64_t& entries, TChain* friendChain=0x0, const Char_t* friendChainFile=0x0); + void FillEventInfo(AliCorrelationReducedEvent* event, Float_t* values, AliCorrelationReducedEventFriend* eventF=0x0); + void FillEventOfflineTriggers(UShort_t triggerBit, Float_t* values); + void FillTrackingFlag(AliCorrelationReducedTrack* track, UShort_t flag, Float_t* values); + void FillTrackInfo(AliCorrelationReducedTrack* p, Float_t* values); + void FillPairInfo(AliCorrelationReducedPair* p, Float_t* values); + void FillPairInfo(AliCorrelationReducedTrack* t1, AliCorrelationReducedTrack* t2, Int_t type, Float_t* values); + void FillCorrelationInfo(AliCorrelationReducedPair* p, AliCorrelationReducedTrack* t, Float_t* values); + void FillCaloClusterInfo(AliCorrelationReducedCaloCluster* cl, Float_t* values); + void FillHistClass(const Char_t* className, Float_t* values); + void DoEventMixing(TList* list1, TList* list2, Float_t* values, Int_t mixingType, const Char_t* histClass); + void EventMixingPairTracks(TList* pairs, TList* tracks, Float_t* values); + void EventMixingResonanceLegs(TList* posLegs, TList* negLegs, Float_t* values, Int_t type, const Char_t* histClass); + Double_t DeltaPhi(Double_t phi1, Double_t phi2); // calculate delta phi in the (-pi,+pi) interval + Bool_t IsPairSelectedEM(Float_t* values); // pair selection used in the mixed event + void AddHistClass(const Char_t* histClass); + Int_t ValidateHistogramName(THashList* hList, const Char_t* name); + void AddHistogram(const Char_t* histClass, + const Char_t* name, const Char_t* title, Bool_t isProfile, + Int_t nXbins, Double_t xmin, Double_t xmax, Int_t varX, + Int_t nYbins=0, Double_t ymin=0, Double_t ymax=0, Int_t varY=kNothing, + Int_t nZbins=0, Double_t zmin=0, Double_t zmax=0, Int_t varZ=kNothing, + const Char_t* xLabels="", const Char_t* yLabels="", const Char_t* zLabels=""); + void AddHistogram(const Char_t* histClass, + const Char_t* name, const Char_t* title, Bool_t isProfile, + Int_t nXbins, Double_t* xbins, Int_t varX, + Int_t nYbins=0, Double_t* ybins=0x0, Int_t varY=kNothing, + Int_t nZbins=0, Double_t* zbins=0x0, Int_t varZ=kNothing, + const Char_t* xLabels="", const Char_t* yLabels="", const Char_t* zLabels=""); + void MakeAxisLabels(TAxis* ax, const Char_t* labels); + void InitFile(const Char_t* filename); // open an old output filename + void CloseFile(); + TObject* GetHistogram(const Char_t* listname, const Char_t* hname); // get a histogram from an old output + +} // end declarations for namespace DstCommonMacros + +//__________________________________________________________________ +void DstCommonMacros::FillEventInfo(AliCorrelationReducedEvent* event, Float_t* values, AliCorrelationReducedEventFriend* eventF/*=0x0*/) { + // + // fill event wise info + // + values[kRunNo] = event->RunNo(); + values[kBC] = event->BC(); + values[kTriggerMask] = event->TriggerMask(); + values[kIsPhysicsSelection] = (event->IsPhysicsSelection() ? 1.0 : 0.0); + values[kVtxX] = event->Vertex(0); + values[kVtxY] = event->Vertex(1); + values[kVtxZ] = event->Vertex(2); + values[kVtxXtpc] = event->VertexTPC(0); + values[kVtxYtpc] = event->VertexTPC(1); + values[kVtxZtpc] = event->VertexTPC(2); + values[kCentVZERO] = event->CentralityVZERO(); + values[kCentSPD] = event->CentralitySPD(); + values[kCentTPC] = event->CentralityTPC(); + values[kCentZDC] = event->CentralityZEMvsZDC(); + values[kCentQuality] = event->CentralityQuality(); + values[kNV0total] = event->NV0CandidatesTotal(); + values[kNV0selected] = event->NV0Candidates(); + values[kNdielectrons] = event->NDielectrons(); + values[kNtracksTotal] = event->NTracksTotal(); + values[kNtracksSelected] = event->NTracks(); + values[kSPDntracklets] = event->SPDntracklets(); + values[kVZEROAemptyChannels] = 0; + values[kVZEROCemptyChannels] = 0; + for(Int_t ich=0;ich<64;++ich) gUsedVars[kVZEROChannelMult+ich] = kTRUE; + Float_t theta=0.0; + for(Int_t ich=0;ich<64;++ich) { + if(gUsedVars[kVZEROChannelMult+ich]) { + values[kVZEROChannelMult+ich] = event->MultChannelVZERO(ich); + if(!gVzeroAvMult[0] && values[kVZEROChannelMult+ich]Qx(AliCorrelationReducedEventFriend::kVZEROA+iVZEROside, ih+1); + values[kVZEROQvecY+iVZEROside*6+ih] = eventF->Qy(AliCorrelationReducedEventFriend::kVZEROA+iVZEROside, ih+1); + values[kVZERORP +iVZEROside*6+ih] = eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROA+iVZEROside, ih+1); + values[kVZEROQvecX+2*6 +ih] += values[kVZEROQvecX+iVZEROside*6+ih]; + values[kVZEROQvecY+2*6 +ih] += values[kVZEROQvecY+iVZEROside*6+ih]; + // cos(n(EPtpc-EPvzero A/C)) + values[kTPCRPres+iVZEROside*6+ih] = DeltaPhi(eventF->EventPlane(AliCorrelationReducedEventFriend::kTPC, ih+1), eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROA+iVZEROside, ih+1)); + values[kTPCRPres+iVZEROside*6+ih] = TMath::Cos(values[kTPCRPres+iVZEROside*6+ih]*(ih+1)); + } + values[kVZERORP +2*6+ih] = TMath::ATan2(values[kVZEROQvecY+2*6+ih],values[kVZEROQvecX+2*6+ih]); + // cos (n*(psi_A-psi_C)) + values[kVZERORPres + ih] = DeltaPhi(eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROA, ih+1), + eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROC, ih+1)); + values[kVZERORPres + ih] = TMath::Cos(values[kVZERORPres + ih]*(ih+1)); + // Qx,Qy correlations for VZERO + values[kVZEROXaXc+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kVZEROA, ih+1)*eventF->Qx(AliCorrelationReducedEventFriend::kVZEROC, ih+1); + values[kVZEROXaYa+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kVZEROA, ih+1)*eventF->Qy(AliCorrelationReducedEventFriend::kVZEROA, ih+1); + values[kVZEROXaYc+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kVZEROA, ih+1)*eventF->Qy(AliCorrelationReducedEventFriend::kVZEROC, ih+1); + values[kVZEROYaXc+ih] = eventF->Qy(AliCorrelationReducedEventFriend::kVZEROA, ih+1)*eventF->Qx(AliCorrelationReducedEventFriend::kVZEROC, ih+1); + values[kVZEROYaYc+ih] = eventF->Qy(AliCorrelationReducedEventFriend::kVZEROA, ih+1)*eventF->Qy(AliCorrelationReducedEventFriend::kVZEROC, ih+1); + values[kVZEROXcYc+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kVZEROC, ih+1)*eventF->Qy(AliCorrelationReducedEventFriend::kVZEROC, ih+1); + // Psi_A - Psi_C + values[kVZEROdeltaRPac+ih] = DeltaPhi(eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROA, ih+1), + eventF->EventPlane(AliCorrelationReducedEventFriend::kVZEROC, ih+1)); + + // TPC event plane + values[kTPCQvecX+ih] = eventF->Qx(AliCorrelationReducedEventFriend::kTPC, ih+1); + values[kTPCQvecY+ih] = eventF->Qy(AliCorrelationReducedEventFriend::kTPC, ih+1); + values[kTPCRP +ih] = eventF->EventPlane(AliCorrelationReducedEventFriend::kTPC, ih+1); + // TPC VZERO Q-vector correlations + values[kRPXtpcXvzeroa+ih] = values[kTPCQvecX+ih]*values[kVZEROQvecX+ih]; + values[kRPXtpcXvzeroc+ih] = values[kTPCQvecX+ih]*values[kVZEROQvecX+6+ih]; + values[kRPYtpcYvzeroa+ih] = values[kTPCQvecY+ih]*values[kVZEROQvecY+ih]; + values[kRPYtpcYvzeroc+ih] = values[kTPCQvecY+ih]*values[kVZEROQvecY+6+ih]; + values[kRPXtpcYvzeroa+ih] = values[kTPCQvecX+ih]*values[kVZEROQvecY+ih]; + values[kRPXtpcYvzeroc+ih] = values[kTPCQvecX+ih]*values[kVZEROQvecY+6+ih]; + values[kRPYtpcXvzeroa+ih] = values[kTPCQvecY+ih]*values[kVZEROQvecX+ih]; + values[kRPYtpcXvzeroc+ih] = values[kTPCQvecY+ih]*values[kVZEROQvecX+6+ih]; + // Psi_TPC - Psi_VZERO A/C + values[kRPdeltaVZEROAtpc+ih] = DeltaPhi(values[kVZERORP+0*6+ih], values[kTPCRP+ih]); + values[kRPdeltaVZEROCtpc+ih] = DeltaPhi(values[kVZERORP+1*6+ih], values[kTPCRP+ih]); + } // end loop over harmonics + + // VZERO v2 using TPC event plane + Double_t vzeroChannelPhi[8] = {0.3927, 1.1781, 1.9635, 2.7489, -2.7489, -1.9635, -1.1781, -0.3927}; + for(Int_t ich=0; ich<64; ++ich) { + if(gUsedVars[kVZEROflowV2TPC]) { + values[kVZEROflowV2TPC+ich] = values[kVZEROChannelMult+ich]*TMath::Cos(2.0*DeltaPhi(vzeroChannelPhi[ich%8],values[kTPCRP+1])); + } + } + } // end if (eventF) +} + + +//_________________________________________________________________ +void DstCommonMacros::FillTrackingFlag(AliCorrelationReducedTrack* track, UShort_t flag, Float_t* values) { + // + // fill the tracking flag + // + values[kTrackingFlag] = -1; + if(track->CheckTrackStatus(flag)) values[kTrackingFlag] = flag; +} + + +//_________________________________________________________________ +void DstCommonMacros::FillEventOfflineTriggers(UShort_t triggerBit, Float_t* values) { + // + // fill the trigger bit input + // + if(triggerBit>=64) return; + if(!gCurrentEvent) return; + ULong64_t trigger = BIT(0); + values[kOfflineTrigger] = triggerBit; + values[kOfflineTriggerFired] = (gCurrentEvent->TriggerMask()&(trigger<0.01 ? triggerBit : -1.0); +} + + +//_________________________________________________________________ +void DstCommonMacros::FillTrackInfo(AliCorrelationReducedTrack* p, Float_t* values) { + // + // fill track information + // + values[kPt] = p->Pt(); + values[kPtTPC] = p->PtTPC(); + if(gUsedVars[kP]) values[kP] = p->P(); + if(gUsedVars[kTheta]) values[kTheta] = p->Theta(); + values[kPhi] = p->Phi(); + values[kPhiTPC] = p->PhiTPC(); + values[kEta] = p->Eta(); + values[kEtaTPC] = p->EtaTPC(); + values[kPin] = p->Pin(); + values[kDcaXY] = p->DCAxy(); + values[kDcaZ] = p->DCAz(); + + if(gUsedVars[kITSncls]) values[kITSncls] = p->ITSncls(); + values[kITSsignal] = p->ITSsignal(); + + values[kTPCncls] = p->TPCncls(); + if(gUsedVars[kTPCnclsRatio]) + values[kTPCnclsRatio] = (p->TPCFindableNcls()>0 ? Float_t(p->TPCncls())/Float_t(p->TPCFindableNcls()) : 0.0); + values[kTPCnclsIter1] = p->TPCnclsIter1(); + values[kTPCnclsF] = p->TPCFindableNcls(); + values[kTPCsignal] = p->TPCsignal(); + + values[kTOFbeta] = p->TOFbeta(); + for(Int_t specie=kElectron; specie<=kProton; ++specie) { + values[kTPCnSig+specie] = p->TPCnSig(specie); + values[kTOFnSig+specie] = p->TOFnSig(specie); + } + values[kTRDpidProbabilities] = p->TRDpid(0); + values[kTRDpidProbabilities+1] = p->TRDpid(1); + + values[kTRDntracklets] = p->TRDntracklets(0); + values[kTRDntrackletsPID] = p->TRDntracklets(1); + + if(gUsedVars[kEMCALmatchedEnergy] || gUsedVars[kEMCALmatchedEOverP]) { + AliCorrelationReducedCaloCluster* cluster = gCurrentEvent->GetCaloCluster(p->CaloClusterId()); + values[kEMCALmatchedEnergy] = (cluster ? cluster->Energy() : -999.0); + Float_t mom = p->P(); + values[kEMCALmatchedEnergy] = (TMath::Abs(mom)>1.e-8 && cluster ? values[kEMCALmatchedEOverP]/mom : -999.0); + } + + // Fill track flow variables + for(Int_t iVZEROside=0; iVZEROside<2; ++iVZEROside) { + for(Int_t ih=0; ih<6; ++ih) { + if(gUsedVars[kTrackVZEROFlowVn+iVZEROside*6+ih] || gUsedVars[kTrackDeltaPhiVZEROFlowVn+iVZEROside*6+ih]) { + values[kTrackVZEROFlowVn+iVZEROside*6+ih] = TMath::Cos(DeltaPhi(values[kPhi],values[kVZERORP+iVZEROside*6+ih])*(ih+1)); + values[kTrackDeltaPhiVZEROFlowVn+iVZEROside*6+ih] = DeltaPhi(values[kPhi],values[kVZERORP+iVZEROside*6+ih]); + } + } + } +} + + +//_________________________________________________________________ +void DstCommonMacros::FillCaloClusterInfo(AliCorrelationReducedCaloCluster* cl, Float_t* values) { + // + // Fill calorimeter cluster information + // + values[kEMCALclusterEnergy] = cl->Energy(); + values[kEMCALclusterDx] = cl->Dx(); + values[kEMCALclusterDz] = cl->Dz(); + values[kEMCALdetector] = (cl->IsEMCAL() ? AliCorrelationReducedCaloCluster::kEMCAL : AliCorrelationReducedCaloCluster::kPHOS); +} + + +//_________________________________________________________________ +void DstCommonMacros::FillPairInfo(AliCorrelationReducedPair* p, Float_t* values) { + // + // fill pair information + // + + values[kCandidateId] = p->CandidateId(); + values[kPairType] = p->PairType(); + if(gUsedVars[kMass]) { + values[kMass] = p->Mass(); + if(p->CandidateId()==AliCorrelationReducedPair::kLambda0ToPPi) values[kMass] = p->Mass(1); + if(p->CandidateId()==AliCorrelationReducedPair::kALambda0ToPPi) values[kMass] = p->Mass(2); + } + values[kPairPt] = p->Pt(); + if(gUsedVars[kPairP]) values[kPairP] = p->P(); + if(gUsedVars[kPairPx]) values[kPairPx] = p->Px(); + if(gUsedVars[kPairPy]) values[kPairPy] = p->Py(); + if(gUsedVars[kPairPz]) values[kPairPz] = p->Pz(); + values[kPairEta] = p->Eta(); + if(gUsedVars[kPairRap]) values[kPairRap] = p->Rapidity(); + values[kPairPhi] = p->Phi(); + values[kPairLxy] = p->Lxy(); + values[kPairOpeningAngle] = p->OpeningAngle(); + if(gUsedVars[kPairTheta]) values[kPairTheta] = p->Theta(); + + // Flow variables + // cos(n*phi), sin(n*phi) + for(Int_t ih=0; ih<6; ++ih) { + if(gUsedVars[kPairCosNPhi+ih]) values[kPairCosNPhi+ih] = TMath::Cos(values[kPairPhi]*(1.0+ih)); + if(gUsedVars[kPairSinNPhi+ih]) values[kPairSinNPhi+ih] = TMath::Sin(values[kPairPhi]*(1.0+ih)); + } + // VZERO + for(Int_t iVZEROside=0; iVZEROside<3; ++iVZEROside) { + // v3 using VZERO Psi_2 + if(gUsedVars[kPairVZEROFlowV3Psi2+iVZEROside]) + values[kPairVZEROFlowV3Psi2+iVZEROside] = TMath::Cos(3.0*DeltaPhi(values[kPairPhi],values[kVZERORP+iVZEROside*6+1])); + for(Int_t ih=0; ih<6; ++ih) { + // vn using VZERO Psi_n + if(gUsedVars[kPairVZEROFlowVn+iVZEROside*6+ih]) + values[kPairVZEROFlowVn+iVZEROside*6+ih] = TMath::Cos(DeltaPhi(values[kPairPhi], values[kVZERORP+iVZEROside*6+ih])*(ih+1)); + // phi - Psi_n + if(gUsedVars[kPairDeltaPhiVZEROFlowVn+iVZEROside*6+ih]) + values[kPairDeltaPhiVZEROFlowVn+iVZEROside*6+ih] = DeltaPhi(values[kPairPhi], values[kVZERORP+iVZEROside*6+ih]); + } + } + // TPC + // v3 using Psi_2 + if(gUsedVars[kPairTPCFlowV3Psi2]) + values[kPairTPCFlowV3Psi2] = TMath::Cos(3.0*DeltaPhi(values[kPairPhi],values[kTPCRP+1])); + for(Int_t ih=0; ih<6; ++ih) { + // vn using Psi_n + if(gUsedVars[kPairTPCFlowVn+ih]) + values[kPairTPCFlowVn+ih] = TMath::Cos(DeltaPhi(values[kPairPhi],values[kTPCRP+ih])*(ih+1)); + if(gUsedVars[kPairDeltaPhiTPCFlowVn+ih]) + values[kPairDeltaPhiTPCFlowVn+ih] = DeltaPhi(values[kPairPhi], values[kTPCRP+ih]); + } +} + + +//_________________________________________________________________ +void DstCommonMacros::FillPairInfo(AliCorrelationReducedTrack* t1, AliCorrelationReducedTrack* t2, Int_t type, Float_t* values) { + // + // fill pair information from 2 tracks + // + // type - Parameter encoding the resonance type + // This is needed for making a mass assumption on the legs + // + if(gUsedVars[kPairType]) { + if(t1->Charge()*t2->Charge()<0) values[kPairType] = 1; + else if(t1->Charge()>0) values[kPairType] = 0; + else values[kPairType] = 2; + } + Float_t kMass1 = 0.0; + Float_t kMass2 = 0.0; + switch (type) { + case AliCorrelationReducedPair::kK0sToPiPi : + kMass1 = 0.13957; kMass2 = 0.13957; + break; + case AliCorrelationReducedPair::kPhiToKK : + kMass1 = 0.493677; kMass2 = 0.493677; + break; + case AliCorrelationReducedPair::kLambda0ToPPi : + kMass1 = 0.938272; kMass2 = 0.13957; + break; + case AliCorrelationReducedPair::kALambda0ToPPi : + kMass1 = 0.13957; kMass2 = 0.938272; + break; + case AliCorrelationReducedPair::kJpsiToEE : + kMass1 = 0.000511; kMass2 = 0.000511; + break; + default : + break; + } + + if(gUsedVars[kMass]) { + values[kMass] = kMass1*kMass1+kMass2*kMass2 + + 2.0*(TMath::Sqrt(kMass1*kMass1+t1->P()*t1->P())*TMath::Sqrt(kMass2*kMass2+t2->P()*t2->P()) - + t1->Px()*t2->Px() - t1->Py()*t2->Py() - t1->Pz()*t2->Pz()); + if(values[kMass]<0.0) { + cout << "FillPairInfo(track, track, type, values): Warning: Very small squared mass found. " + << " Could be negative due to resolution of Float_t so it will be set to a small positive value." << endl; + cout << " mass2: " << values[kMass] << endl; + cout << "p1(p,x,y,z): " << t1->P() << ", " << t1->Px() << ", " << t1->Py() << ", " << t1->Pz() << endl; + cout << "p2(p,x,y,z): " << t2->P() << ", " << t2->Px() << ", " << t2->Py() << ", " << t2->Pz() << endl; + values[kMass] = 0.0; + } + else + values[kMass] = TMath::Sqrt(values[kMass]); + } + + if(gUsedVars[kPairPt]) values[kPairPt] = TMath::Sqrt((t1->Px()+t2->Px())*(t1->Px()+t2->Px()) + + (t1->Py()+t2->Py())*(t1->Py()+t2->Py())); + if(gUsedVars[kPairP]) values[kPairP] = TMath::Sqrt((t1->Px()+t2->Px())*(t1->Px()+t2->Px()) + + (t1->Py()+t2->Py())*(t1->Py()+t2->Py()) + + (t1->Pz()+t2->Pz())*(t1->Pz()+t2->Pz())); + if(gUsedVars[kPairEta]) { + Float_t p = TMath::Sqrt((t1->Px()+t2->Px())*(t1->Px()+t2->Px()) + + (t1->Py()+t2->Py())*(t1->Py()+t2->Py()) + + (t1->Pz()+t2->Pz())*(t1->Pz()+t2->Pz())); + values[kPairEta] = p-t1->Pz()-t2->Pz(); + values[kPairEta] = (TMath::Abs(values[kPairEta])>1.0e-8 ? (p+t1->Pz()+t2->Pz())/values[kPairEta] : 0.0); + values[kPairEta] = (values[kPairEta]>1.0e-8 ? 0.5*TMath::Log(values[kPairEta]) : -999.); + } + if(gUsedVars[kPairRap]) { + Float_t mass = kMass1*kMass1+kMass2*kMass2 + + 2.0*(TMath::Sqrt(kMass1*kMass1+t1->P()*t1->P())*TMath::Sqrt(kMass2*kMass2+t2->P()*t2->P()) - + t1->Px()*t2->Px() - t1->Py()*t2->Py() - t1->Pz()*t2->Pz()); + if(mass<0.0) { + cout << "Negative squared mass (Float_t resolution). Setting mass to zero" << endl; + mass = 0.0; + } + else mass = TMath::Sqrt(mass); + Float_t e = TMath::Sqrt(mass*mass+ + (t1->Px()+t2->Px())*(t1->Px()+t2->Px()) + + (t1->Py()+t2->Py())*(t1->Py()+t2->Py()) + + (t1->Pz()+t2->Pz())*(t1->Pz()+t2->Pz())); + values[kPairRap] = 0.5*TMath::Log((e+t1->Pz()+t2->Pz())/(e-t1->Pz()-t2->Pz())); + } + if(gUsedVars[kPairPhi]) { + values[kPairPhi] = TMath::ATan2(t1->Py()+t2->Py(),t1->Px()+t2->Px()); + if(values[kPairPhi]<0.0) values[kPairPhi] = 2.0*TMath::Pi() + values[kPairPhi]; + } + if(gUsedVars[kPairTheta]) values[kPairTheta] = TMath::ACos((t1->Pz()+t2->Pz())/ + TMath::Sqrt((t1->Px()+t2->Px())*(t1->Px()+t2->Px()) + + (t1->Py()+t2->Py())*(t1->Py()+t2->Py()) + + (t1->Pz()+t2->Pz())*(t1->Pz()+t2->Pz()))); + // Flow variables + // cos(n*phi), sin(n*phi) + for(Int_t ih=0; ih<6; ++ih) { + if(gUsedVars[kPairCosNPhi+ih]) values[kPairCosNPhi+ih] = TMath::Cos(values[kPairPhi]*(1.0+ih)); + if(gUsedVars[kPairSinNPhi+ih]) values[kPairSinNPhi+ih] = TMath::Sin(values[kPairPhi]*(1.0+ih)); + } + // VZERO + for(Int_t iVZEROside=0; iVZEROside<3; ++iVZEROside) { + // v3 using VZERO Psi_2 + if(gUsedVars[kPairVZEROFlowV3Psi2+iVZEROside]) + values[kPairVZEROFlowV3Psi2+iVZEROside] = TMath::Cos(3.0*DeltaPhi(values[kPairPhi],values[kVZERORP+iVZEROside*6+1])); + for(Int_t ih=0; ih<6; ++ih) { + // vn using VZERO Psi_n + if(gUsedVars[kPairVZEROFlowVn+iVZEROside*6+ih]) + values[kPairVZEROFlowVn+iVZEROside*6+ih] = TMath::Cos(DeltaPhi(values[kPairPhi], values[kVZERORP+iVZEROside*6+ih])*(ih+1)); + // phi - Psi_n + if(gUsedVars[kPairDeltaPhiVZEROFlowVn+iVZEROside*6+ih]) + values[kPairDeltaPhiVZEROFlowVn+iVZEROside*6+ih] = DeltaPhi(values[kPairPhi], values[kVZERORP+iVZEROside*6+ih]); + } + } + // TPC + // v3 using Psi_2 + if(gUsedVars[kPairTPCFlowV3Psi2]) + values[kPairTPCFlowV3Psi2] = TMath::Cos(3.0*DeltaPhi(values[kPairPhi],values[kTPCRP+1])); + for(Int_t ih=0; ih<6; ++ih) { + // vn using Psi_n + if(gUsedVars[kPairTPCFlowVn+ih]) + values[kPairTPCFlowVn+ih] = TMath::Cos(DeltaPhi(values[kPairPhi],values[kTPCRP+ih])*(ih+1)); + if(gUsedVars[kPairDeltaPhiTPCFlowVn+ih]) + values[kPairDeltaPhiTPCFlowVn+ih] = DeltaPhi(values[kPairPhi], values[kTPCRP+ih]); + } +} + + +//__________________________________________________________________ +void DstCommonMacros::FillCorrelationInfo(AliCorrelationReducedPair* p, AliCorrelationReducedTrack* t, Float_t* values) { + // + // fill pair-track correlation information + // + if(gUsedVars[kDeltaPhi]) + values[kDeltaPhi] = DeltaPhi(p->Phi(), t->Phi()); + + if(gUsedVars[kDeltaTheta]) + values[kDeltaTheta] = p->Theta() - t->Theta(); + + if(gUsedVars[kDeltaEta]) values[kDeltaEta] = p->Eta() - t->Eta(); +} + + +//__________________________________________________________________ +void DstCommonMacros::DoEventMixing(TList* list1, TList* list2, Float_t* values, Int_t type, const Char_t* histClass) { + // + // Do the event mixing + // + if(type==0) return; + + cout << "mixing ..." << endl; + switch(type) { + case -1: + EventMixingPairTracks(list1, list2, values); // list1 contains pairs; list2 contains tracks + break; + case AliCorrelationReducedPair::kJpsiToEE: + // type is needed to encode the type of resonance which is needed for the mass assumption of legs when calculating the invariant mass + EventMixingResonanceLegs(list1, list2, values, type, histClass); // list1 contains positive legs; list2 contains negative legs + break; + default: + break; + } + + for(Int_t i=list1->GetEntries()-1; i>=0; --i) ((TList*)list1->At(i))->RemoveAll(); + for(Int_t j=list2->GetEntries()-1; j>=0; --j) ((TList*)list2->At(j))->RemoveAll(); + list1->Clear(); + list2->Clear(); +} + + +//__________________________________________________________________ +void DstCommonMacros::EventMixingPairTracks(TList* pairLists, TList* trackLists, Float_t* values) { + // + // Make event mixing for pair - track correlations + // + cout << "Mixing category " << gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data() << endl; + Int_t entries = pairLists->GetEntries(); // we should have the same number of entries in both pair and track lists + if(entries<2) return; + + TList* pairs = 0x0; + TList* tracks = 0x0; + AliCorrelationReducedPair* pair = 0x0; + AliCorrelationReducedTrack* track = 0x0; + + TIter nextPairList(pairLists); + for(Int_t ie1=0; ie1GetEntries(); ++ip) { + pair = (AliCorrelationReducedPair*)nextPair(); + FillPairInfo(pair, values); + + TIter nextTrack(tracks); + Int_t leadingIdx = -1; + Float_t leadingPt = 0.0; + for(Int_t it=0; itGetEntries(); ++it) { + track = (AliCorrelationReducedTrack*)nextTrack(); + FillTrackInfo(track, values); + + FillCorrelationInfo(pair, track, values); + if(pair->PairType()==1) + FillHistClass(Form("Correlation_ME_US%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values); + if(pair->PairType()==0) { + FillHistClass(Form("Correlation_ME_LS%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()), values); + FillHistClass(Form("Correlation_ME_LSpp%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()), values); + } + if(pair->PairType()==2) { + FillHistClass(Form("Correlation_ME_LS%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()), values); + FillHistClass(Form("Correlation_ME_LSmm%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()), values); + } + if(track->Pt()>leadingPt) { + leadingPt = track->Pt(); + leadingIdx = it; + } + } // end loop over tracks + if(leadingIdx!=-1) { + track = (AliCorrelationReducedTrack*)tracks->At(leadingIdx); + FillTrackInfo(track, values); + FillCorrelationInfo(pair, track, values); + if(pair->PairType()==1) + FillHistClass(Form("Correlation_LeadingPt_ME_US%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values); + if(pair->PairType()==0) { + FillHistClass(Form("Correlation_LeadingPt_ME_LS%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values); + FillHistClass(Form("Correlation_LeadingPt_ME_LSpp%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values); + } + if(pair->PairType()==2) { + FillHistClass(Form("Correlation_LeadingPt_ME_LS%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values); + FillHistClass(Form("Correlation_LeadingPt_ME_LSmm%s", gEMCategoryNames[TMath::Nint(values[kEventMixingId])].Data()),values); + } + } + } // end loop over pairs + } // end loop over second event + } // end loop over first event +} + + +//__________________________________________________________________ +void DstCommonMacros::EventMixingResonanceLegs(TList* posLists, TList* negLists, Float_t* values, Int_t type, const Char_t* histClass) { + // + // Do event mixing with the legs of a resonance to obtain the background + // + cout << "Mixing event class " << histClass << ", (cent/vtx): " << values[kCentVZERO] << "/" << values[kVtxZ] << endl; + Int_t entries = posLists->GetEntries(); + cout << "mixing positives: " << entries << endl; + cout << "mixing negatives: " << negLists->GetEntries() << endl; + if(entries<2) return; + + TList* positives1 = 0x0; // list of tracks in the first event + //TList* negatives1 = 0x0; + //TList* positives2 = 0x0; // list of tracks in the second event + TList* negatives2 = 0x0; + AliCorrelationReducedTrack* posTrack1 = 0x0; + //AliCorrelationReducedTrack* negTrack1 = 0x0; + //AliCorrelationReducedTrack* posTrack2 = 0x0; + AliCorrelationReducedTrack* negTrack2 = 0x0; + + TIter nextPosList1(posLists); + TIter nextNegList1(negLists); + for(Int_t ie1=0; ie12.0*TMath::Pi()) delta -= 2.0*TMath::Pi(); + if(delta<0.0) delta += 2.0*TMath::Pi(); + /*Double_t delta = phi2; + if(phi2<0.0) delta += 2.0*TMath::Pi(); + delta = phi1-delta; + if(delta>TMath::Pi()) delta = delta - 2.0*TMath::Pi(); + if(delta<-1.*TMath::Pi()) delta = 2.0*TMath::Pi() + delta; + */ + return delta; +} + + +//__________________________________________________________________ +void DstCommonMacros::FillHistClass(const Char_t* className, Float_t* values) { + // + // fill a class of histograms + // + THashList* hList = (THashList*)gHistLists->FindObject(className); + if(!hList) { + //cout << "Warning in FillHistClass(): Histogram class " << className << " not found" << endl; + return; + } + + TIter next(hList); + TObject* h=0x0; + while((h=next())) { + UInt_t id = h->GetUniqueID(); + UInt_t varX = id%kNVars; + UInt_t varY = (id/kNVars)%kNVars; + UInt_t varZ = (id/kNVars/kNVars)%kNVars; + if(idFill(values[varX]); + } else if(idGetDimension()==1) ((TProfile*)h)->Fill(values[varX],values[varY]); + if(((TH1*)h)->GetDimension()==2) ((TH2F*)h)->Fill(values[varX],values[varY]); + } else { + if(((TH1*)h)->GetDimension()==2) ((TProfile2D*)h)->Fill(values[varX],values[varY],values[varZ]); + if(((TH1*)h)->GetDimension()==3) ((TH3F*)h)->Fill(values[varX],values[varY],values[varZ]); + }; + } +} + + +//__________________________________________________________________ +void DstCommonMacros::AddHistClass(const Char_t* histClass) { + // + // Add a new histogram list + // + if(!gHistLists) + gHistLists = new TObjArray(); + gHistLists->SetOwner(); + gHistLists->SetName("histos"); + + if(gHistLists->FindObject(histClass)) { + cout << "Warning in AddHistClass: Cannot add histogram class " << histClass + << " because it already exists." << endl; + return; + } + THashList* hList=new THashList; + hList->SetOwner(kTRUE); + hList->SetName(histClass); + gHistLists->Add(hList); +} + + +//_________________________________________________________________ +Int_t DstCommonMacros::ValidateHistogramName(THashList* hList, const Char_t* name) { + // + // check whether a histogram with this name already exist, and how many of them + // + for(Int_t i=0; iGetEntries(); ++i) { + TString hname = hList->At(i)->GetName(); + TObjArray* arr=hname.Tokenize("#"); + TString nameRoot = arr->At(0)->GetName(); + if(nameRoot.CompareTo(name)==0) { + cout << "Warning in AddHistogram(): Histogram " << name << " already exists in this list (" << nameRoot.Data() << ")" << endl; + return -1; + } + } + Int_t nnames = 0; + for(Int_t il=0; ilGetEntries(); ++il) { + THashList* tlist = (THashList*)gHistLists->At(il); + for(Int_t ih=0; ihGetEntries(); ++ih) { + TString hname = tlist->At(ih)->GetName(); + TObjArray* arr=hname.Tokenize("#"); + TString nameRoot = arr->At(0)->GetName(); + if(nameRoot.CompareTo(name)==0) ++nnames; + } + } + return nnames; +} + + +//_________________________________________________________________ +void DstCommonMacros::AddHistogram(const Char_t* histClass, + const Char_t* name, const Char_t* title, Bool_t isProfile, + Int_t nXbins, Double_t xmin, Double_t xmax, Int_t varX, + Int_t nYbins, Double_t ymin, Double_t ymax, Int_t varY, + Int_t nZbins, Double_t zmin, Double_t zmax, Int_t varZ, + const Char_t* xLabels, const Char_t* yLabels, const Char_t* zLabels) { + // + // add a histogram + // + THashList* hList = (THashList*)gHistLists->FindObject(histClass); + if(hList->FindObject(name)) { + cout << "Warning in AddHistogram(): Histogram " << name << " already exists" << endl; + return; + } + /*Int_t nnames = ValidateHistogramName(hList, name); + if(nnames<0) return; + TString hname = Form("%s#v%d", name, nnames); + cout << "Assigned name " << hname.Data() << endl;*/ + TString hname = name; + + Int_t dimension = 1; + if(varY!=kNothing) dimension = 2; + if(varZ!=kNothing) dimension = 3; + + TString titleStr(title); + TObjArray* arr=titleStr.Tokenize(";"); + + TH1* h=0x0; + switch(dimension) { + case 1: + h=new TH1F(hname.Data(),arr->At(0)->GetName(),nXbins,xmin,xmax); + //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX) << endl; + h->SetUniqueID(UInt_t(varX)); + if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName()); + if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels); + gUsedVars[varX] = kTRUE; + hList->Add(h); + break; + case 2: + if(isProfile) + h=new TProfile(hname.Data(),arr->At(0)->GetName(),nXbins,xmin,xmax); + else + h=new TH2F(hname.Data(),arr->At(0)->GetName(),nXbins,xmin,xmax,nYbins,ymin,ymax); + //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX+kNVars*varY) << endl; + h->SetUniqueID(UInt_t(varX+kNVars*varY)); + if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName()); + if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels); + if(arr->At(2)) h->GetYaxis()->SetTitle(arr->At(2)->GetName()); + if(yLabels[0]!='\0') MakeAxisLabels(h->GetYaxis(), yLabels); + gUsedVars[varX] = kTRUE; + gUsedVars[varY] = kTRUE; + hList->Add(h); + break; + case 3: + if(isProfile) + h=new TProfile2D(hname.Data(),arr->At(0)->GetName(),nXbins,xmin,xmax,nYbins,ymin,ymax); + else + h=new TH3F(hname.Data(),arr->At(0)->GetName(),nXbins,xmin,xmax,nYbins,ymin,ymax,nZbins,zmin,zmax); + //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX+kNVars*varY+kNVars*kNVars*varZ) << endl; + h->SetUniqueID(UInt_t(varX+kNVars*varY+kNVars*kNVars*varZ)); + if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName()); + if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels); + if(arr->At(2)) h->GetYaxis()->SetTitle(arr->At(2)->GetName()); + if(yLabels[0]!='\0') MakeAxisLabels(h->GetYaxis(), yLabels); + if(arr->At(3)) h->GetZaxis()->SetTitle(arr->At(3)->GetName()); + if(zLabels[0]!='\0') MakeAxisLabels(h->GetZaxis(), zLabels); + gUsedVars[varX] = kTRUE; + gUsedVars[varY] = kTRUE; + gUsedVars[varZ] = kTRUE; + hList->Add(h); + break; + } +} + +//_________________________________________________________________ +void DstCommonMacros::AddHistogram(const Char_t* histClass, + const Char_t* name, const Char_t* title, Bool_t isProfile, + Int_t nXbins, Double_t* xbins, Int_t varX, + Int_t nYbins, Double_t* ybins, Int_t varY, + Int_t nZbins, Double_t* zbins, Int_t varZ, + const Char_t* xLabels, const Char_t* yLabels, const Char_t* zLabels) { + // + // add a histogram + // + THashList* hList = (THashList*)gHistLists->FindObject(histClass); + if(hList->FindObject(name)) { + cout << "Warning in AddHistogram(): Histogram " << name << " already exists" << endl; + return; + } + //Int_t nnames = ValidateHistogramName(hList, name); + //if(nnames<0) return; + //TString hname = Form("%s#v%d", name, nnames); + TString hname = name; + + Int_t dimension = 1; + if(varY!=kNothing) dimension = 2; + if(varZ!=kNothing) dimension = 3; + + TString titleStr(title); + TObjArray* arr=titleStr.Tokenize(";"); + + TH1* h=0x0; + switch(dimension) { + case 1: + h=new TH1F(hname.Data(),arr->At(0)->GetName(),nXbins,xbins); + //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX) << endl; + h->SetUniqueID(UInt_t(varX)); + if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName()); + if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels); + gUsedVars[varX] = kTRUE; + hList->Add(h); + break; + case 2: + if(isProfile) + h=new TProfile(hname.Data(),arr->At(0)->GetName(),nXbins,xbins); + else + h=new TH2F(hname.Data(),arr->At(0)->GetName(),nXbins,xbins,nYbins,ybins); + //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX+kNVars*varY) << endl; + h->SetUniqueID(UInt_t(varX+kNVars*varY)); + if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName()); + if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels); + if(arr->At(2)) h->GetYaxis()->SetTitle(arr->At(2)->GetName()); + if(yLabels[0]!='\0') MakeAxisLabels(h->GetYaxis(), yLabels); + gUsedVars[varX] = kTRUE; + gUsedVars[varY] = kTRUE; + hList->Add(h); + break; + case 3: + if(isProfile) + h=new TProfile2D(hname.Data(),arr->At(0)->GetName(),nXbins,xbins,nYbins,ybins); + else + h=new TH3F(hname.Data(),arr->At(0)->GetName(),nXbins,xbins,nYbins,ybins,nZbins,zbins); + //cout << "HISTOGRAM " << hname.Data() << ", ID: " << UInt_t(varX+kNVars*varY+kNVars*kNVars*varZ) << endl; + h->SetUniqueID(UInt_t(varX+kNVars*varY+kNVars*kNVars*varZ)); + if(arr->At(1)) h->GetXaxis()->SetTitle(arr->At(1)->GetName()); + if(xLabels[0]!='\0') MakeAxisLabels(h->GetXaxis(), xLabels); + if(arr->At(2)) h->GetYaxis()->SetTitle(arr->At(2)->GetName()); + if(yLabels[0]!='\0') MakeAxisLabels(h->GetYaxis(), yLabels); + if(arr->At(3)) h->GetZaxis()->SetTitle(arr->At(3)->GetName()); + if(zLabels[0]!='\0') MakeAxisLabels(h->GetZaxis(), zLabels); + gUsedVars[varX] = kTRUE; + gUsedVars[varY] = kTRUE; + gUsedVars[varZ] = kTRUE; + hList->Add(h); + break; + } +} + + +//____________________________________________________________________________________ +void DstCommonMacros::MakeAxisLabels(TAxis* ax, const Char_t* labels) { + // + // add bin labels to an axis + // + TString labelsStr(labels); + TObjArray* arr=labelsStr.Tokenize(";"); + for(Int_t ib=1; ib<=ax->GetNbins(); ++ib) { + if(ib>=arr->GetEntries()+1) break; + ax->SetBinLabel(ib, arr->At(ib-1)->GetName()); + } +} + + +//____________________________________________________________________________________ +void DstCommonMacros::InitFile(const Char_t* filename) { + // + // Open an existing ROOT file containing lists of histograms and initialize the global list pointer + // + //if(gHistFile) gHistFile->Close(); + TString histfilename=""; + if(gHistFile) histfilename = gHistFile->GetName(); + if(!histfilename.Contains(filename)) + gHistFile = new TFile(filename); // open file only if not already open + + if(!gHistFile) { + cout << "GlobalMacros.C::GetHistogram() : File " << filename << " not opened!!" << endl; + return; + } + if(gHistFile->IsZombie()) { + cout << "GlobalMacros.C::GetHistogram() : File " << filename << " not opened!!" << endl; + return; + } + TList* list1 = gHistFile->GetListOfKeys(); + TKey* key1 = (TKey*)list1->At(0); + gHistListsOld = (TDirectoryFile*)key1->ReadObj(); +} + + +//____________________________________________________________________________________ +void DstCommonMacros::CloseFile() { + // + // Close the opened file + // + gHistListsOld = 0x0; + if(gHistFile && gHistFile->IsOpen()) gHistFile->Close(); +} + + +//____________________________________________________________________________________ +TObject* DstCommonMacros::GetHistogram(const Char_t* listname, const Char_t* hname) { + // + // Retrieve a histogram from the list hlist + // + if(!gHistListsOld) { + cout << "GlobalMacros.C::GetHistogram() : The main list was not initialized." << endl; + cout << " A ROOT file must pe initialized first!!" << endl; + return 0x0; + } + TKey* listKey = gHistListsOld->FindKey(listname); + TDirectoryFile* hlist = (TDirectoryFile*)listKey->ReadObj(); + TKey* key = hlist->FindKey(hname); + return key->ReadObj(); +} + + +//_________________________________________________________________ +TChain* DstCommonMacros::GetChain(const Char_t* filename, Int_t howMany, Int_t offset, Long64_t& entries, + TChain* friendChain/*=0x0*/, const Char_t* friendChainFile/*=0x0*/) { + // + // read an ascii file containing a list of root files with reduced trees + // and build a TChain + // + cout << "Creating the data chain from " << filename << " ..." << flush; + TChain* chain = new TChain("DstTree"); + ifstream inBuf; + inBuf.open(filename); + Int_t index = 0; + while(inBuf.good()) { + Char_t str[512]; + inBuf.getline(str,512,'\n'); + + if(index=offset+howMany) break; + + TString strstr = str; + if(!strstr.Contains(".root")) continue; + cout << endl << "Adding file " << str << endl; + chain->Add(str); + if(friendChain) { + TObjArray* arr = strstr.Tokenize("/"); + strstr.ReplaceAll(arr->At(arr->GetEntries()-1)->GetName(),friendChainFile); + friendChain->Add(strstr.Data()); + } + ++index; + } + inBuf.close(); + entries = chain->GetEntries(); + Long64_t entriesFriend = (friendChain ? friendChain->GetEntries() : 0); + cout << "DstCommonMacros::GetChain() Chain entries = " << entries << endl; + if(friendChain) + cout << "DstCommonMacros::GetChain() Friend chain entries = " << entriesFriend << endl; + if(friendChain && (entries!=entriesFriend)) { + cout << "DstCommonMacros::GetChain() The friend chain does not have the same number of entries as the main chain!!!" << endl; + cout << " Check it out and retry !!" << endl; + return 0x0; + } + cout << " done" << endl; + return chain; +} + + +//__________________________________________________________________ +void DstCommonMacros::WriteOutput(TFile* save) { + // + // Write the histogram lists in the output file + // + cout << "Writing the output to " << save->GetName() << " ... " << flush; + //TFile* save=new TFile(filename,"RECREATE"); + TDirectory* mainDir = save->mkdir(gHistLists->GetName()); + mainDir->cd(); + for(Int_t i=0; iGetEntries(); ++i) { + THashList* list = (THashList*)gHistLists->At(i); + TDirectory* dir = mainDir->mkdir(list->GetName()); + dir->cd(); + list->Write(); + mainDir->cd(); + } + save->Close(); + cout << "done" << endl; +} + + +//__________________________________________________________________ +Bool_t DstCommonMacros::IsPairSelectedEM(Float_t* values) { + // + // pair selection on the pair resulting from the event mixing + // + if(values[kPairPt]gkJpsiPtCut[1]) return kFALSE; + return kTRUE; +} diff --git a/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/RunDstPbPbJpsiAnalysis.C b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/RunDstPbPbJpsiAnalysis.C new file mode 100644 index 00000000000..526356c88ef --- /dev/null +++ b/PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/RunDstPbPbJpsiAnalysis.C @@ -0,0 +1,928 @@ +// J/psi analysis in PbPb +// author: Ionut-Cristian Arsene, i.c.arsene@gsi.de +// 2012/Apr/11 + +#include +using namespace std; + +#include + +#ifndef ALICORRELATIONREDUCEDEVENT_H +#include "AliCorrelationReducedEvent.h" +#endif + +#include "DstCommonMacros.C" +using namespace DstCommonMacros; + +// function prototypes +void DefineHistograms(const Char_t* histClasses, const Char_t* output); +Bool_t IsEventSelected(AliCorrelationReducedEvent* event); +Bool_t IsLegSelected(AliCorrelationReducedTrack* track, Bool_t* legQualityMap); +Bool_t IsLegSelected(AliCorrelationReducedTrack* track); +Bool_t IsLegSelectedTRD(AliCorrelationReducedTrack* track, Int_t trdCut=-1, Int_t trdMinNtr=4, Float_t eleCut=0.7); +Bool_t IsLegSelectedTOF(AliCorrelationReducedTrack* track, Int_t tofCut=-1, Float_t lowExcl=0.3, Float_t highExcl=0.93); +Bool_t IsTrackUsedForMixing(UShort_t* idxArr, UShort_t size, UShort_t idx); +Bool_t IsPairSelected(AliCorrelationReducedPair* pair); +void RunDstPbPbJpsiAnalysis(const Char_t* output, const Char_t* inputfilename, + Int_t howMany/* = 5000*/, Int_t offset/* = 0*/); + +TFile* gSaveFile=0x0; +// event plane friend file +const Char_t* gkEPfile="dstTree_VZERO_TPC_recentered.root"; + +// mixing pool depth +const Int_t gkEventMixingPoolDepth = 100; + +// centrality binning for the event mixing +const Int_t gkNCentRanges = 3; +/*Double_t gkEMCentLims[gkNCentRanges+1] = { 0.0, 2.5, 5.0, 7.5, 10.0, + 12.5, 15.0, 17.5, 20.0, 22.5, + 25.0, 27.5, 30.0, 32.5, 35.0, + 37.5, 40.0, 50.0, 60.0, 70.0, 80.0};*/ +Double_t gkEMCentLims[gkNCentRanges+1] = { 0.0, 10.0, 40.0, 80.0}; +// vertex binning for the event mixing +const Int_t gkNVtxRanges = 1; +//Double_t gkEMVtxLims[gkNVtxRanges+1] = {-10.0, -6.0, -2.0, 2.0, 6.0, 10.0}; +Double_t gkEMVtxLims[gkNVtxRanges+1] = {-10.0, 10.0}; + +// event plane angle binning for the event mixing +const Int_t gkNPhiEPRanges = 1; +Double_t gkEMPhiEPLims[gkNPhiEPRanges+1] = {-1.5708, 1.5708}; +//Double_t gkEMPhiEPLims[gkNPhiEPRanges+1] = {-1.5708, -1.2566, -0.94248, -0.62832, -0.31415926, 0.0, 0.31415926, 0.62832, 0.94248, 1.2566, 1.5708}; +/*Double_t gkEMPhiEPLims[gkNPhiEPRanges+1] = {-1.5708, -1.4137, -1.2566, -1.0996, -0.94248, -0.78540, -0.62832, -0.47124, -0.31416, -0.15708, 0.0, + 0.15708, 0.31416, 0.47124, 0.62832, 0.7854, 0.94248, 1.09956, 1.2566, 1.4137, 1.5708};*/ +// which event plane to be used +const Int_t gkEPused = kVZERORP + 6*1 + 1; // VZEROC 2nd harmonic event plane +//const Int_t gkEPused = kVZERORP + 6*0 + 1; // VZEROA 2nd harmonic event plane +//const Int_t gkEPused = kTPCRP + 1; // TPC harmonic event plane + +// define detector cuts +const Int_t gkNDetCuts = 6; +const Char_t* gkLegDetCutNames[gkNDetCuts] = { + "TPC", + "TPC_TRD4_1_0.7", "TPC_TRD4_1_0.8", "TPC_TRD4_1_0.9", + "TPC_TOF", "TPC_TRD4_1_0.9_TOF" +}; + +//_________________________________________________________________________________________________ +void RunDstPbPbJpsiAnalysis(const Char_t* output, const Char_t* inputfilename, + Int_t howMany/* = 5000*/, Int_t offset/* = 0*/) { + // + // J/psi analysis in Pb-Pb + // + cout << "start ..." << endl; + TTimeStamp start; + cout << "creating chain ..." << endl; + // create the input chains ----------------------------------------- + Long64_t entries=0; + TChain* friendChain=0x0; + if(gkEPfile[0]!='\0') + friendChain = new TChain("DstFriendTree"); + TChain* chain = GetChain(inputfilename, howMany, offset, entries, friendChain, gkEPfile); + if(!chain) return; + if(gkEPfile[0]!='\0' && !friendChain) return; + AliCorrelationReducedEvent* event = new AliCorrelationReducedEvent(); + chain->SetBranchAddress("Event",&event); + AliCorrelationReducedEventFriend* eventF = 0x0; + if(gkEPfile[0]!='\0') { + eventF = new AliCorrelationReducedEventFriend(); + friendChain->SetBranchAddress("Event",&eventF); + } + + cout << "initialize event mixing lists ..." << endl; + + // define histograms ----------------------------------------------- + TString histClasses = ""; + histClasses += "Event_BeforeCuts;Event_AfterCuts;VZERORP;TPCRP;"; + histClasses += "OfflineTriggers;"; + histClasses += "TrackingFlags_Before;TrackQA_JpsiLegs_BeforeCuts_ITS_TPC_TRD_TOF;"; + for(Int_t iLegCut=0; iLegCutSetOwner(); + selectedNegLegs[iDetCut][iCent][iVtx][iPhi][i] = new TList(); + selectedNegLegs[iDetCut][iCent][iVtx][iPhi][i]->SetOwner(); + } + } // end loop over phi bins + } // end loop over vtx bins + } // end loop over centrality bins + } // end loop over leg cuts + DefineHistograms(histClasses.Data(), output); + + Float_t values[kNVars]; + Int_t trackIdMap[20000] = {-1}; + TClonesArray* trackList=0x0; + TClonesArray* pairList=0x0; + Bool_t legsQuality[2][gkNDetCuts] = {{kFALSE}}; + Bool_t allLegCutsOr[2] = {kFALSE}; + Bool_t pairQuality[gkNDetCuts] = {kFALSE}; + + Float_t oldVertex = -999.; Float_t oldBC=-999.; Float_t oldCentVZERO=-999.; Float_t oldCentSPD=-999.; Float_t oldCentTPC=-999.; + Double_t mixingTime = 0.0; + TTimeStamp startEventLoop; + + for(Int_t ie=0; ieGetEntry(ie); + friendChain->GetEntry(ie); + + gCurrentEvent = event; + gCurrentEventFriend = eventF; + if(ie%100==0) cout << "event " << ie << endl; + + FillEventInfo(event, values, eventF); + FillHistClass("Event_BeforeCuts", values); + // event cuts + if(!IsEventSelected(event)) continue; + // check wheter this event is the same as the previous + if(TMath::Abs(event->Vertex(2)-oldVertex)<0.0001 && TMath::Abs(event->CentralityVZERO()-oldCentVZERO)<0.001) { + cout << "This event is the a copy of the previous event" << endl; + cout << "OLD/NEW: vtxZ " << oldVertex << "/" << event->Vertex(2) + << " BC " << oldBC << "/" << event->BC() + << " CentVZERO " << oldCentVZERO << "/" << event->CentralityVZERO() + << " CentSPD " << oldCentSPD << "/" << event->CentralitySPD() + << " CentTPC " << oldCentTPC << "/" << event->CentralityTPC() << endl; + ++ie; + continue; + } + else { + oldVertex = event->Vertex(2); oldBC = event->BC(); + oldCentVZERO = event->CentralityVZERO(); oldCentSPD = event->CentralitySPD(); + oldCentTPC = event->CentralityTPC(); + } + + for(UShort_t ibit=0; ibit<64; ++ibit) { + FillEventOfflineTriggers(ibit, values); + FillHistClass("OfflineTriggers", values); + } + + // get the event vtx and centrality bins + Float_t centVZERO = values[kCentVZERO]; + Float_t vtxZ = values[kVtxZ]; + Int_t binCent = -1; + for(Int_t i=0; i=gkEMCentLims[i] && centVZERO=gkEMVtxLims[i] && vtxZGetTracks(); + TIter nextTrack(trackList); + for(Int_t i=0; i<20000; ++i) trackIdMap[i] = -1; + for(Int_t it=0; itNTracks(); ++it) { + track = (AliCorrelationReducedTrack*)nextTrack(); + if(track) + trackIdMap[track->TrackId()] = it; + } + + pairList = event->GetPairs(); + TIter nextPair(pairList); + AliCorrelationReducedPair* pair = 0x0; + AliCorrelationReducedTrack* leg1 = 0x0; + AliCorrelationReducedTrack* leg2 = 0x0; + values[kNpairsSelected] = 0.0; + + // reset arrays + for(Int_t iCut=0; iCutNDielectrons(); ++ip) { + pair = event->GetDielectronPair(ip); + // pair cuts + if(!IsPairSelected(pair)) continue; + + // reset flag arrays + for(Int_t iCut=0; iCutLegId(0)]!=-1) { + leg1 = event->GetTrack(trackIdMap[pair->LegId(0)]); + FillTrackInfo(leg1, values); + for(UShort_t iflag=0; iflagLegId(1)]!=-1) { + leg2 = event->GetTrack(trackIdMap[pair->LegId(1)]); + FillTrackInfo(leg2, values); + for(UShort_t iflag=0; iflagTrackId())) { + FillHistClass(Form("TrackQA_JpsiLegs_%s_AfterCuts_ITS_TPC_TRD_TOF", gkLegDetCutNames[iCut]), values); + for(UShort_t iflag=0; iflagTrackId())) { + FillHistClass(Form("TrackQA_JpsiLegs_%s_AfterCuts_ITS_TPC_TRD_TOF", gkLegDetCutNames[iCut]), values); + for(UShort_t iflag=0; iflagTrackId())) { + if(leg1->Charge()>0) selectedPosLegs[iCut][binCent][binVtxZ][binPhiEP][mixingPoolSize[iCut][binCent][binVtxZ][binPhiEP]]->Add(leg1->Clone()); + else selectedNegLegs[iCut][binCent][binVtxZ][binPhiEP][mixingPoolSize[iCut][binCent][binVtxZ][binPhiEP]]->Add(leg1->Clone()); + tracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP][nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP]] = leg1->TrackId(); + nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP] += 1; + } + if(!IsTrackUsedForMixing(tracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP], + nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP], + leg2->TrackId())) { + if(leg2->Charge()>0) selectedPosLegs[iCut][binCent][binVtxZ][binPhiEP][mixingPoolSize[iCut][binCent][binVtxZ][binPhiEP]]->Add(leg2->Clone()); + else selectedNegLegs[iCut][binCent][binVtxZ][binPhiEP][mixingPoolSize[iCut][binCent][binVtxZ][binPhiEP]]->Add(leg2->Clone()); + tracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP][nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP]] = leg2->TrackId(); + nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP] += 1; + } + eventHasCandidates[iCut][binCent][binVtxZ][binPhiEP] = kTRUE; + } // end if(pairQuality) + } + } // end loop over pairs + + // If this event has candidates, then store the event in the list for event mixing + for(Int_t iCut=0; iCutAdd(selectedPosLegs[iCut][iCent][iVtx][iPhi][mixingPoolSize[iCut][iCent][iVtx][iPhi]]); + list2[iCut][iCent][iVtx][iPhi]->Add(selectedNegLegs[iCut][iCent][iVtx][iPhi][mixingPoolSize[iCut][iCent][iVtx][iPhi]]); + /* + cout << "event mixing category (cut/cent/vtx) " << iCut << "/" << iCent << "/" << iVtx << ", pool size (pos/neg): " + << list1[iCut][iCent][iVtx]->GetEntries() + << "/" << list2[iCut][iCent][iVtx]->GetEntries() << "; ntracks this entry (pos/neg): " + << selectedPosLegs[iCut][iCent][iVtx][mixingPoolSize[iCut][iCent][iVtx]]->GetEntries() << "/" + << selectedNegLegs[iCut][iCent][iVtx][mixingPoolSize[iCut][iCent][iVtx]]->GetEntries() << endl; + */ + mixingPoolSize[iCut][iCent][iVtx][iPhi] += 1; + if(mixingPoolSize[iCut][iCent][iVtx][iPhi]==gkEventMixingPoolDepth) + fullListFound = kTRUE; + } // end loop over phi ranges + } // end loop over vertex ranges + } // end loop over centrality ranges + } // end loop over detector cuts + + FillHistClass("Event_AfterCuts", values); + FillHistClass("TPCRP", values); + FillHistClass("VZERORP", values); + } // end loop over events + + // Mixing the leftover events + cout << "Leftover mixing ..." << endl; + TTimeStamp startMix; + for(Int_t iCut=0; iCutIsPhysicsSelection()) return kFALSE; + if(event->VertexNContributors()<1) return kFALSE; + if(event->CentralityQuality()!=0) return kFALSE; + if(event->CentralityVZERO()<-0.0001) return kFALSE; + if(event->CentralityVZERO()>80.0) return kFALSE; + if(TMath::Abs(event->Vertex(2))>10.0) return kFALSE; + return kTRUE; +} + + +//__________________________________________________________________ +Bool_t IsLegSelected(AliCorrelationReducedTrack* track, Bool_t* legQuality) { + // + // track leg selection + // + legQuality[0] = IsLegSelected(track); + legQuality[1] = legQuality[0] && IsLegSelectedTRD(track, 1, 4, 0.7); + legQuality[2] = legQuality[0] && IsLegSelectedTRD(track, 1, 4, 0.80); + legQuality[3] = legQuality[0] && IsLegSelectedTRD(track, 1, 4, 0.90); + legQuality[4] = legQuality[0] && IsLegSelectedTOF(track, 2); + legQuality[5] = legQuality[3] && legQuality[4]; + Bool_t globalOr = kFALSE; + for(Int_t i=0; iCheckTrackStatus(kTPCrefit)) return kFALSE; + if(!track->CheckTrackStatus(kITSrefit)) return kFALSE; + if(!(track->ITSLayerHit(0) || track->ITSLayerHit(1))) return kFALSE; + + if(track->Pt()<0.7) return kFALSE; + if(TMath::Abs(track->Eta())>0.9) return kFALSE; + + if(track->TPCncls()<70) return kFALSE; + + if(track->TPCnSig(kElectron)>3.0) return kFALSE; + if(track->TPCnSig(kElectron)<-2.0) return kFALSE; + + if(TMath::Abs(track->DCAz())>3.0) return kFALSE; + if(TMath::Abs(track->DCAxy())>1.0) return kFALSE; + + //if(track->Pin()<2.5) { + if(TMath::Abs(track->TPCnSig(kProton))<3.5) return kFALSE; + if(TMath::Abs(track->TPCnSig(kPion))<3.5) return kFALSE; + if(track->TPCnSig(kProton)<-3.5) return kFALSE; + //if(track->TPCnSig(kProton)<-3.5) return kFALSE; + //} + + return kTRUE; +} + + +//__________________________________________________________________ +Bool_t IsLegSelectedTRD(AliCorrelationReducedTrack* track, + Int_t trdCut/*=-1*/, Int_t trdMinNtr/*=4*/, Float_t eleCut/*=0.7*/) { + // + // TRD leg selection + // + if(trdCut==1 && track->TRDntracklets(0)>=trdMinNtr) { + if(track->TRDpid(0)TRDntracklets(0)>=trdMinNtr && track->Pin()>1.0) { + if(track->TRDpid(0)CheckTrackStatus(kTRDpid) && track->TRDntracklets(0)>=trdMinNtr) { + if(track->TRDpid(0)CheckTrackStatus(kTRDpid) && track->TRDntracklets(1)>=trdMinNtr && track->Pin()>1.0) { + if(track->TRDpid(0)TOFbeta()>lowExcl && track->TOFbeta()CheckTrackStatus(kTOFpid) && TMath::Abs(track->TOFnSig(kElectron))>3.5) return kFALSE; + return kTRUE; +} + + +//__________________________________________________________________ +Bool_t IsPairSelected(AliCorrelationReducedPair* pair) { + // + // pair selection + // + if(!pair) return kFALSE; + if(pair->CandidateId()!=AliCorrelationReducedPair::kJpsiToEE) return kFALSE; + if(pair->PairType()!=1) return kFALSE; + if(TMath::Abs(pair->Rapidity())>0.9) return kFALSE; + return kTRUE; +} + + +//__________________________________________________________________ +void DefineHistograms(const Char_t* histClasses, const Char_t* output) { + // + // define the histograms + // + cout << "Defining histograms ..." << flush; + cout << "histogram classes: " << histClasses << endl; + + gSaveFile=new TFile(output,"RECREATE"); + + TString classesStr(histClasses); + TObjArray* arr=classesStr.Tokenize(";"); + + const Int_t kNRunBins = 3000; + Double_t runHistRange[2] = {137000.,140000.}; + + for(Int_t iclass=0; iclassGetEntries(); ++iclass) { + TString classStr = arr->At(iclass)->GetName(); + cout << "hist class: " << classStr.Data() << endl; + + // Event wise histograms + if(classStr.Contains("Event")) { + cout << "Event" << endl; + AddHistClass(classStr.Data()); + AddHistogram(classStr.Data(),"RunNo","Run numbers;Run",kFALSE, kNRunBins, runHistRange[0], runHistRange[1], kRunNo); + AddHistogram(classStr.Data(),"BC","Bunch crossing;BC",kFALSE,3000,0.,3000.,kBC); + AddHistogram(classStr.Data(),"IsPhysicsSelection","Physics selection flag;;",kFALSE, + 2,-0.5,1.5,kIsPhysicsSelection, 0,0.0,0.0,kNothing, 0,0.0,0.0,kNothing, "off;on"); + + AddHistogram(classStr.Data(),"VtxZ","Vtx Z;vtx Z (cm)",kFALSE,300,-15.,15.,kVtxZ); + AddHistogram(classStr.Data(),"VtxX","Vtx X;vtx X (cm)",kFALSE,300,-0.4,0.4,kVtxX); + AddHistogram(classStr.Data(),"VtxY","Vtx Y;vtx Y (cm)",kFALSE,300,-0.4,0.4,kVtxY); + + AddHistogram(classStr.Data(),"VtxZ_Run_prof", " vs run; Run; vtxZ (cm)", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, -20, 20., 1000., kVtxZ); + AddHistogram(classStr.Data(),"VtxX_Run_prof", " vs run; Run; vtxX (cm)", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, -20, 20., 1000., kVtxX); + AddHistogram(classStr.Data(),"VtxY_Run_prof", " vs run; Run; vtxY (cm)", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, -20, 20., 1000., kVtxY); + + AddHistogram(classStr.Data(),"CentVZERO","Centrality(VZERO);centrality VZERO (percents)",kFALSE, + 100, 0.0, 100.0, kCentVZERO); + AddHistogram(classStr.Data(),"CentQuality","Centrality quality;centrality quality",kFALSE, + 100, -50.5, 49.5, kCentQuality); + AddHistogram(classStr.Data(),"CentVZERO_Run_prof"," vs run;Run; centrality VZERO (%)",kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0.0, 100.0, kCentVZERO); + + AddHistogram(classStr.Data(),"NPairs","Number of candidates per event;# pairs",kFALSE, + 5000,0.,5000.,kNdielectrons); + AddHistogram(classStr.Data(),"NPairsSelected", "Number of selected pairs per event; #pairs", kFALSE, + 5000,0.,5000.,kNpairsSelected); + AddHistogram(classStr.Data(),"NTracksTotal","Number of total tracks per event;# tracks",kFALSE, + 1000,0.,30000.,kNtracksTotal); + AddHistogram(classStr.Data(),"NTracksSelected","Number of selected tracks per event;# tracks",kFALSE, + 1000,0.,30000.,kNtracksSelected); + AddHistogram(classStr.Data(),"SPDntracklets", "SPD #tracklets in |#eta|<1.0; tracklets", kFALSE, + 3000, -0.5, 2999.5, kSPDntracklets); + + AddHistogram(classStr.Data(),"Ndielectrons_Run_prof", " per run; Run; #tracks", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0., 10000., kNdielectrons); + AddHistogram(classStr.Data(),"NpairsSelected_Run_prof", " per run; Run; #tracks", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0., 10000., kNpairsSelected); + AddHistogram(classStr.Data(),"NTracksTotal_Run_prof", " per run; Run; #tracks", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0., 10000., kNtracksTotal); + AddHistogram(classStr.Data(),"NTracksSelected_Run_prof", " per run; Run; #tracks", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0., 10000., kNtracksSelected); + AddHistogram(classStr.Data(),"SPDntracklets_Run_prof", " per run; Run; #tracks", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0., 10000., kSPDntracklets); + + AddHistogram(classStr.Data(),"VtxZ_CentVZERO","Centrality(VZERO) vs vtx. Z;vtx Z (cm); centrality VZERO (%)",kFALSE, + 300,-15.,15.,kVtxZ, 100, 0.0, 100.0, kCentVZERO); + AddHistogram(classStr.Data(),"VtxZ_CentSPD","Centrality(SPD) vs vtx. Z;vtx Z (cm); centrality SPD (%)",kFALSE, + 300,-15.,15.,kVtxZ, 100, 0.0, 100.0, kCentSPD); + AddHistogram(classStr.Data(),"VtxZ_CentTPC","Centrality(TPC) vs vtx. Z;vtx Z (cm); centrality TPC (%)",kFALSE, + 300,-15.,15.,kVtxZ, 100, 0.0, 100.0, kCentTPC); + continue; + } // end if className contains "Event" + + // Offline trigger histograms + if(classStr.Contains("OfflineTriggers")) { + cout << "OfflineTriggers" << endl; + AddHistClass(classStr.Data()); + + TString triggerNames = ""; + for(Int_t i=0; i<64; ++i) {triggerNames += gkOfflineTriggerNames[i]; triggerNames+=";";} + + AddHistogram(classStr.Data(), "Triggers", "Offline triggers fired; ; ;", kFALSE, + 64, -0.5, 63.5, kOfflineTrigger, 2, -0.5, 1.5, kOfflineTriggerFired, 0, 0.0, 0.0, kNothing, triggerNames.Data(), "off;on"); + AddHistogram(classStr.Data(), "Triggers2", "Offline triggers fired; ; ;", kFALSE, + 64, -0.5, 63.5, kOfflineTriggerFired2, 0, 0.0, 0.0, kNothing, 0, 0.0, 0.0, kNothing, triggerNames.Data()); + AddHistogram(classStr.Data(), "CentVZERO_Triggers2", "Offline triggers fired vs centrality VZERO; ; centrality VZERO;", kFALSE, + 64, -0.5, 63.5, kOfflineTriggerFired2, 20, 0.0, 100.0, kCentVZERO, 0, 0.0, 0.0, kNothing, triggerNames.Data()); + AddHistogram(classStr.Data(), "VtxZ_Triggers2", "Offline triggers fired vs vtxZ; ; vtx Z (cm.);", kFALSE, + 64, -0.5, 63.5, kOfflineTriggerFired2, 200, -20.0, 20.0, kVtxZ, 0, 0.0, 0.0, kNothing, triggerNames.Data()); + continue; + } + + if(classStr.Contains("VZERORP")) { + cout << "VZERORP" << endl; + AddHistClass(classStr.Data()); + const Char_t* sname[3] = {"A","C","A&C"}; + for(Int_t ih=1; ih<2; ++ih) { + for(Int_t iS=0; iS<3; ++iS) { + AddHistogram(classStr.Data(), Form("QvecX_side%s_h%d_CentSPDVtxZ_prof",sname[iS],ih+1), + Form("Q_{x}, side %s, harmonic %d, vs centSPD and vtxZ; Centrality (percents); VtxZ (cm); ",sname[iS],ih+1), kTRUE, + 20, 0.0, 100.0, kCentSPD, 24, -12.0, 12.0, kVtxZ, 500, -10000.0, 10000.0, kVZEROQvecX+iS*6+ih); + AddHistogram(classStr.Data(), Form("QvecY_side%s_h%d_CentSPDVtxZ_prof",sname[iS],ih+1), + Form("Q_{y}, side %s, harmonic %d, vs centSPD and vtxZ; Centrality (percents); VtxZ (cm); ",sname[iS],ih+1), kTRUE, + 20, 0.0, 100.0, kCentSPD, 24, -12.0, 12.0, kVtxZ, 500, -10000.0, 10000.0, kVZEROQvecY+iS*6+ih); + AddHistogram(classStr.Data(), Form("QvecX_side%s_h%d_Run_prof", sname[iS], ih+1), + Form(", VZERO side %s, harmonic %d, vs run; Run; ", sname[iS], ih+1), kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, -100., 100., kVZEROQvecX+iS*6+ih); + AddHistogram(classStr.Data(), Form("QvecY_side%s_h%d_Run_prof", sname[iS], ih+1), + Form(", VZERO side %s, harmonic %d, vs run; Run; ", sname[iS], ih+1), kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, -100., 100., kVZEROQvecY+iS*6+ih); + AddHistogram(classStr.Data(), Form("RP_side%s_h%d_CentSPD",sname[iS],ih+1), + Form("VZERO reaction plane, side %s, harmonic %d, vs centrality SPD; #Psi (rad.); centSPD (percents)",sname[iS],ih+1), kFALSE, + 400, -4.0/Double_t(ih+1), 4.0/Double_t(ih+1), kVZERORP+iS*6+ih, 20, 0.0, 100.0, kCentSPD); + AddHistogram(classStr.Data(), Form("RP_side%s_h%d_VtxZ",sname[iS],ih+1), + Form("VZERO reaction plane, side %s, harmonic %d, vs vtxZ; #Psi (rad.); vtxZ (cm)",sname[iS],ih+1), kFALSE, + 400, -4.0/Double_t(ih+1), 4.0/Double_t(ih+1), kVZERORP+iS*6+ih, 24, -12.0, +12.0, kVtxZ); + } // end loop over VZERO sides + } // end loop over harmonics + continue; + } // end if for the VZERO reaction plane histograms + + if(classStr.Contains("TPCRP")) { + cout << "TPCRP" << endl; + AddHistClass(classStr.Data()); + for(Int_t ih=1; ih<2; ++ih) { + AddHistogram(classStr.Data(), Form("QvecX_TPC_h%d_Run_prof", ih+1), + Form(", TPC, harmonic %d, vs run; Run; ", ih+1), kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, -100., 100., kTPCQvecX+ih); + AddHistogram(classStr.Data(), Form("QvecY_TPC_h%d_Run_prof", ih+1), + Form(", TPC, harmonic %d, vs run; Run; ", ih+1), kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, -100., 100., kTPCQvecY+ih); + AddHistogram(classStr.Data(), Form("TPCRP_h%d", ih+1), + Form("TPC event plane, harmonic %d; #Psi (rad.)", ih+1), kFALSE, + 400, -4.0/Double_t(ih+1), 4.0/Double_t(ih+1), kTPCRP+ih); + AddHistogram(classStr.Data(), Form("TPCQvecX_h%d_CentSPDVtxZ_prof", ih+1), + Form("TPC Q_{x}, harmonic %d, vs centSPD and vtxZ; Centrality (percents); VtxZ (cm); ", ih+1), kTRUE, + 20, 0.0, 100.0, kCentSPD, 24, -12.0, 12.0, kVtxZ, 500, -600.0, 600.0, kTPCQvecX+ih); + AddHistogram(classStr.Data(), Form("TPCQvecY_h%d_CentSPDVtxZ_prof", ih+1), + Form("TPC Q_{y}, harmonic %d, vs centSPD and vtxZ; Centrality (percents); VtxZ (cm); ", ih+1), kTRUE, + 20, 0.0, 100.0, kCentSPD, 24, -12.0, 12.0, kVtxZ, 500, -600.0, 600.0, kTPCQvecY+ih); + } // end loop over harmonics + continue; + } // end if for the TPC reaction plane histograms + + TString trkFlagNames = ""; + for(Int_t iflag=0; iflag vs run; run;", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, 0.0, 50.0, kPt); + AddHistogram(classStr.Data(), "Eta_Run", "<#eta> vs run; run;", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, -1.5, 1.5, kEta); + AddHistogram(classStr.Data(), "Phi_Run", "<#varphi> vs run; run;", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, 0.0, 6.3, kPhi); + AddHistogram(classStr.Data(), "DCAxy_Run", " vs run; run;", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, -10.0, 10.0, kDcaXY); + AddHistogram(classStr.Data(), "DCAz_Run", " vs run; run;", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, -10.0, 10.0, kDcaZ); + + // correlations between parameters + AddHistogram(classStr.Data(), "Eta_Pt_prof", " vs #eta; #eta; p_{T} (GeV/c);", kTRUE, + 300, -1.5, +1.5, kEta, 100, 0.0, 10.0, kPt); + AddHistogram(classStr.Data(), "Phi_Pt_prof", " vs #varphi; #varphi (rad.); p_{T} (GeV/c)", kTRUE, + 300, 0.0, 6.3, kPhi, 100, 0.0, 10.0, kPt); + AddHistogram(classStr.Data(), "Eta_Phi", "#varphi vs #eta; #eta; #varphi (rad.);", kFALSE, + 200, -1.0, +1.0, kEta, 100, 0.0, 6.3, kPhi); + AddHistogram(classStr.Data(), "Pt_DCAxy", "DCAxy vs p_{T}; p_{T} (GeV/c); DCA_{xy} (cm)", kFALSE, + 100, 0.0, 10.0, kPt, 500, -2.0, 2.0, kDcaXY); + AddHistogram(classStr.Data(), "Pt_DCAz", "DCAz vs p_{T}; p_{T} (GeV/c); DCA_{z} (cm)", kFALSE, + 100, 0.0, 10.0, kPt, 500, -2.0, 2.0, kDcaZ); + AddHistogram(classStr.Data(), "Eta_DCAxy", "DCAxy vs #eta; #eta; DCA_{xy} (cm)", kFALSE, + 100, -1.0, 1.0, kEta, 500, -2.0, 2.0, kDcaXY); + AddHistogram(classStr.Data(), "Eta_DCAz", "DCAz vs #eta; #eta; DCA_{z} (cm)", kFALSE, + 100, -1.0, 1.0, kEta, 500, -2.0, 2.0, kDcaZ); + + if(classStr.Contains("ITS")) { + AddHistogram(classStr.Data(),"ITSncls", "ITS nclusters;# clusters ITS", kFALSE, + 7,-0.5,6.5,kITSncls); + AddHistogram(classStr.Data(),"ITSncls_Run", "ITS vs run;run;# clusters ITS", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 7,-0.5,6.5, kITSncls); + AddHistogram(classStr.Data(),"ITSncls_Cent_prof", "ITS vs centrality; centrality;# clusters ITS", kTRUE, + 20, 0.0, 100.0, kCentVZERO, 7,-0.5,6.5, kITSncls); + AddHistogram(classStr.Data(),"Pt_ITSncls","ITS nclusters vs p_{T};p_{T} (GeV/c);# clusters ITS",kFALSE, + 100, 0.0, 10.0, kPt, 7,-0.5,6.5, kITSncls); + AddHistogram(classStr.Data(),"Eta_Phi_ITSncls_prof","ITS vs (#eta,#phi);#eta;#phi (rad.);# clusters ITS",kTRUE, + 192, -1.2, 1.2, kEta, 126, 0.0, 6.3, kPhi, 7, -0.5, 6.5, kITSncls); + } // end if ITS histograms + + if(classStr.Contains("TPC")) { + AddHistogram(classStr.Data(),"TPCncls","TPC nclusters;# clusters TPC",kFALSE, + 160,-0.5,159.5,kTPCncls); + AddHistogram(classStr.Data(),"TPCncls_Run","TPC nclusters vs run;run;# clusters TPC",kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 160,-0.5,159.5, kTPCncls); + AddHistogram(classStr.Data(),"TPCncls_CentVZERO","TPC nclusters vs centrality;centrality;# clusters TPC",kFALSE, + 20, 0.0, 100.0, kCentVZERO, 160,-0.5,159.5, kTPCncls); + AddHistogram(classStr.Data(),"TPCncls_Pin","TPC nclusters vs inner param P; P (GeV/c); # clusters TPC",kFALSE, + 200,0.0,20,kPin,160,-0.5,159.5,kTPCncls); + AddHistogram(classStr.Data(),"Eta_Phi_TPCncls_prof","TPC vs (#eta,#phi);#eta;#phi (rad.);# clusters TPC",kTRUE, + 192, -1.2, 1.2, kEta, 126, 0.0, 6.3, kPhi, 160, -0.5, 159.5, kTPCncls); + AddHistogram(classStr.Data(),"TPCsignal_Pin","TPC dE/dx vs. inner param P;P (GeV/c);TPC dE/dx",kFALSE, + 1000,0.0,20.0,kPin,151,-0.5,150.5,kTPCsignal); + AddHistogram(classStr.Data(),"TPCsignal_TPCclusters","TPC dE/dx vs. TPC nclusters;TPC #clusters;TPC dE/dx",kFALSE, + 160,-0.5,159.5,kTPCncls,151,-0.5,150.5,kTPCsignal); + AddHistogram(classStr.Data(),"TPCnsigElectron_Pin","TPC N_{#sigma} electron vs. inner param P;P (GeV/c);N_{#sigma}",kFALSE, + 1000,0.0,20.0,kPin,100,-5.0,5.0,kTPCnSig+kElectron); + AddHistogram(classStr.Data(),"TPCnsigElectron_CentEta_prof","TPC N_{#sigma} electron vs. (#eta,centrality); #eta; centrality VZERO;N_{#sigma}",kTRUE, + 100,-1.0,1.0,kEta, 20, 0.0, 100.0, kCentVZERO, 100,-5.0,5.0,kTPCnSig+kElectron); + AddHistogram(classStr.Data(),"TPCnsigElectron_Run","TPC N_{#sigma} electron vs. run; run; N_{#sigma}",kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo,100,-5.0,5.0,kTPCnSig+kElectron); + } // end if TPC histograms + + if(classStr.Contains("TRD")) { + AddHistogram(classStr.Data(),"TRDntracklets","TRD ntracklets; #tracklets TRD",kFALSE, + 7,-0.5,6.5,kTRDntracklets); + AddHistogram(classStr.Data(),"TRDntrackletsPID","TRD ntracklets PID; #tracklets TRD",kFALSE, + 7,-0.5,6.5,kTRDntrackletsPID); + AddHistogram(classStr.Data(),"TRDprobabElectron","TRD electron probability; probability",kFALSE, + 500,0.0,1.0,kTRDpidProbabilities); + AddHistogram(classStr.Data(),"TRDprobabPion","TRD pion probability; probability",kFALSE, + 500,0.0,1.0,kTRDpidProbabilities+1); + AddHistogram(classStr.Data(),"TRDntracklets_Run","TRD vs run; run; #tracklets TRD",kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 7,-0.5,6.5,kTRDntracklets); + AddHistogram(classStr.Data(),"TRDntrackletsPID_Run","TRD vs run; run; #tracklets TRD",kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 7,-0.5,6.5,kTRDntrackletsPID); + AddHistogram(classStr.Data(),"TRDprobabEle_Run","TRD vs run; run; probability",kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 10,0.0,1.0,kTRDpidProbabilities); + AddHistogram(classStr.Data(),"TRDprobabPion_Run","TRD vs run; run; probability",kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 10,0.0,1.0,kTRDpidProbabilities+1); + AddHistogram(classStr.Data(),"Eta_Phi_TRDntracklets_prof","TRD vs (#eta,#phi);#eta;#phi (rad.);#tracklets TRD",kTRUE, + 192, -1.2, 1.2, kEta, 126, 0.0, 6.3, kPhi, 7, -0.5, 6.5, kTRDntracklets); + AddHistogram(classStr.Data(),"TRDntracklets_cent","TRD ntracklets vs centrality; #tracklets TRD; centrality (%)",kFALSE, + 7,-0.5,6.5,kTRDntracklets,20,0.0, 100.0, kCentSPD); + AddHistogram(classStr.Data(),"Eta_Phi_TRDntrackletsPID_prof","TRD vs (#eta,#phi);#eta;#phi (rad.);#tracklets TRD",kTRUE, + 192, -1.2, 1.2, kEta, 126, 0.0, 6.3, kPhi, 7, -0.5, 6.5, kTRDntrackletsPID); + AddHistogram(classStr.Data(),"TRDntrackletsPID_cent","TRD ntracklets PID vs centrality; #tracklets TRD; centrality (%)",kFALSE, + 7,-0.5,6.5,kTRDntrackletsPID,20,0.0, 100.0, kCentSPD); + AddHistogram(classStr.Data(),"TRDntracklets_TRDntrackletsPID","TRD ntracklets vs TRD ntracklets PID; #tracklets TRD; #tracklets TRD pid",kFALSE, + 7,-0.5,6.5,kTRDntracklets,7,-0.5, 6.5, kTRDntrackletsPID); + } // end if TRD histograms + + if(classStr.Contains("TOF")) { + AddHistogram(classStr.Data(),"TOFbeta_P","TOF #beta vs P;P (GeV/c);#beta",kFALSE, + 200,0.0,20.0,kP, 220,0.0,1.1,kTOFbeta); + AddHistogram(classStr.Data(),"Eta_Phi_TOFbeta_prof","TOF <#beta> vs (#eta,#phi);#eta;#phi (rad.);TOF #beta",kTRUE, + 192, -1.2, 1.2, kEta, 126, 0.0, 6.3, kPhi, 160, -0.5, 159.5, kTOFbeta); + AddHistogram(classStr.Data(),"TOFnsigElectron_P","TOF N_{#sigma} electron vs. P;P (GeV/c);N_{#sigma}",kFALSE, + 200,0.0,20.0,kP, 100,-5.0,5.0,kTOFnSig+kElectron); + AddHistogram(classStr.Data(),"TOFnSigElectron_Run","TOF vs run number;run;TOF n-#sigma_{e}",kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 160, -0.5, 159.5, kTOFnSig+kElectron); + } // end if TOF histograms + continue; + } // end if "TrackQA" + + Double_t massBinWidth = 0.04; // *GeV/c^2 + Double_t massRange[2] = {0.0,6.0}; + const Int_t nMassBins = TMath::Nint((massRange[1]-massRange[0])/massBinWidth); + Double_t massBinLims[nMassBins+1]; for(Int_t im=0; im<=nMassBins; ++im) massBinLims[im] = massBinWidth*im; + + // Histograms for pairs + if(classStr.Contains("PairQA")) { + AddHistClass(classStr.Data()); + + AddHistogram(classStr.Data(), "Mass_CentVZEROCep", + "Inv. mass vs (centVZERO, #Psi^{2}); centrality VZERO; #Psi^{2}; m (GeV/c^{2})", + kFALSE, gkNCentRanges, gkEMCentLims, kCentVZERO, gkNPhiEPRanges, gkEMPhiEPLims, gkEPused, nMassBins, massBinLims, kMass); + AddHistogram(classStr.Data(), "Pt", "Pt; p_{T} (GeV/c)", kFALSE, + 1000, 0.0, 10.0, kPairPt); + + if(classStr.Contains("SE")) { + AddHistogram(classStr.Data(), "CandidateId", "Candidate id; p_{T} (GeV/c)", kFALSE, + AliCorrelationReducedPair::kNMaxCandidateTypes+1, -0.5, Double_t(AliCorrelationReducedPair::kNMaxCandidateTypes)+0.5, kCandidateId); + AddHistogram(classStr.Data(), "PairType", "Pair type; pair type", kFALSE, + 4, -0.5, 3.5, kPairType); + AddHistogram(classStr.Data(), "Mass", "Invariant mass; m_{inv} (GeV/c^{2})", kFALSE, + nMassBins, massRange[0], massRange[1], kMass); + AddHistogram(classStr.Data(), "Rapidity", "Rapidity; y", kFALSE, + 240, -1.2, 1.2, kPairRap); + AddHistogram(classStr.Data(), "Eta", "Pseudo-rapidity #eta; #eta", kFALSE, + 240, -2.0, 2.0, kPairEta); + AddHistogram(classStr.Data(), "Phi", "Azimuthal distribution; #phi (rad.)", kFALSE, + 315, 0.0, 6.3, kPairPhi); + AddHistogram(classStr.Data(), "OpeningAngle", "Opening angle; op. angle (rad)", kFALSE, + 1000, 0.0, 3.2, kPairOpeningAngle); + AddHistogram(classStr.Data(), "Mass_Pt", "Pt vs invariant mass; m_{inv} (GeV/c^{2}); p_{T} (GeV/c)", kFALSE, + nMassBins, massRange[0], massRange[1], kMass, 100, 0.0, 10.0, kPairPt); + AddHistogram(classStr.Data(), "Mass_Run_prof", " vs run; run;m_{inv} (GeV/c^{2})", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, nMassBins, massRange[0], massRange[1], kMass); + AddHistogram(classStr.Data(), "Pt_Run_prof", " vs run; run;p_{T} (GeV/c)", kTRUE, + kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, 0.0, 10.0, kPairPt); + AddHistogram(classStr.Data(), "VZEROflow_sideA_v2_Mass", + "Pair v2 coefficient using VZERO-A reaction plane; inv.mass (GeV/c^{2});v_{2}(#Psi_{2}^{VZERO-A})", kTRUE, + nMassBins, massRange[0], massRange[1], kMass, 1000, -1., 1., kPairVZEROFlowVn+0*6+1); + AddHistogram(classStr.Data(), "VZEROflow_sideC_v2_Mass", + "Pair v2 coefficient using VZERO-C reaction plane; inv.mass (GeV/c^{2});v_{2}(#Psi_{2}^{VZERO-C})", kTRUE, + nMassBins, massRange[0], massRange[1], kMass, 1000, -1., 1., kPairVZEROFlowVn+1*6+1); + AddHistogram(classStr.Data(), "TPCflow_v2_Mass", + "Pair v2 coefficient using TPC reaction plane; inv.mass (GeV/c^{2});v_{2}(#Psi_{2}^{TPC})", kTRUE, + nMassBins, massRange[0], massRange[1], kMass, 1000, -1., 1., kPairTPCFlowVn+1); + } + continue; + } // end if for Pair classes of histograms + + } // end loop over histogram classes +} -- 2.43.0