From: miweber Date: Thu, 22 Mar 2012 19:55:06 +0000 (+0000) Subject: Added analysis task for * Event by event PID fluctuation analysis X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=a6fd3cfe1f879bab2cc24ae89e45a3ab12fc6479;ds=inline Added analysis task for * Event by event PID fluctuation analysis --- diff --git a/PWGCF/EBYE/PIDFluctuation/macros/AddAnalysisTaskPIDFluctuation.C b/PWGCF/EBYE/PIDFluctuation/macros/AddAnalysisTaskPIDFluctuation.C new file mode 100644 index 00000000000..9eb4b93d01b --- /dev/null +++ b/PWGCF/EBYE/PIDFluctuation/macros/AddAnalysisTaskPIDFluctuation.C @@ -0,0 +1,94 @@ +//__________________________________________________________________ + +AliAnalysisTaskPIDFluctuation * +AddAnalysisTaskPIDFluctuation(Int_t aodFilterBit, Float_t ptMin, Float_t ptMax, Float_t etaMin, Float_t etaMax) +{ + + /* init analysis name */ + TString analysisName = "PIDFluctuation"; + analysisName += "_"; + analysisName += GetAODFilterBitName(aodFilterBit); + analysisName += "_"; + analysisName += Form("pt_%.1f_%.1f", ptMin, ptMax); + analysisName += "_"; + analysisName += Form("eta_%.1f_%.1f", etaMin, etaMax); + + /* check analysis manager */ + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error("", "cannot get analysis manager"); + return NULL; + } + + /* check input event handler */ + if (!mgr->GetInputEventHandler()) { + Error("", "cannot get input event handler"); + return NULL; + } + + /* get common input data container */ + AliAnalysisDataContainer *inputc = mgr->GetCommonInputContainer(); + if (!inputc) { + Error("", "cannot get common input container"); + return NULL; + } + + /* create output data container */ + AliAnalysisDataContainer *outputc1 = mgr->CreateContainer(analysisName.Data(), TList::Class(), AliAnalysisManager::kOutputContainer, "PIDFluctuation.root"); + if (!outputc1) { + Error("", "cannot create output container \"Histos\""); + return NULL; + } + + /* create task and connect input/output */ + AliAnalysisTaskPIDFluctuation *task = new AliAnalysisTaskPIDFluctuation(analysisName.Data()); + mgr->AddTask(task); + mgr->ConnectInput(task, 0, inputc); + mgr->ConnectOutput(task, 1, outputc1); + + /* setup task */ + task->SetESDtrackCuts(GetESDtrackCuts(aodFilterBit)); + task->SetAODfilterBit(aodFilterBit); + task->SetEtaRange(etaMin, etaMax); + task->SetPtRange(ptMin, ptMax); + + task->Dump(); + return task; +} + +//__________________________________________________________________ + +AliESDtrackCuts * +GetESDtrackCuts(Int_t type) +{ + AliESDtrackCuts *trackCuts; + switch (type) { + case AliAODTrack::kTrkGlobal: + trackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(); + break; + case AliAODTrack::kTrkTPCOnlyConstrained: + case AliAODTrack::kTrkTPCOnly: + trackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts(); + trackCuts->SetMinNClustersTPC(70); + break; + } + return trackCuts; +} + +//__________________________________________________________________ + +Char_t * +GetAODFilterBitName(Int_t type) +{ + switch (type) { + case AliAODTrack::kTrkGlobal: + return "TrkGlobal"; + break; + case AliAODTrack::kTrkTPCOnlyConstrained: + return "TrkTPCOnlyConstrained"; + break; + case AliAODTrack::kTrkTPCOnly: + return "TrkTPCOnly"; + break; + } +} diff --git a/PWGCF/EBYE/PIDFluctuation/macros/SteerAODAnalysisTaskPIDFluctuation.C b/PWGCF/EBYE/PIDFluctuation/macros/SteerAODAnalysisTaskPIDFluctuation.C new file mode 100644 index 00000000000..05046e0099e --- /dev/null +++ b/PWGCF/EBYE/PIDFluctuation/macros/SteerAODAnalysisTaskPIDFluctuation.C @@ -0,0 +1,70 @@ +SteerAODAnalysisTaskPIDFluctuation(const Char_t *inputfilename, Int_t maxFiles = kMaxInt, Int_t maxEv = kMaxInt) +{ + + /* include path for ACLic */ + gSystem->AddIncludePath("-I$ALICE_ROOT/include"); + gSystem->AddIncludePath("-I$ALICE_ROOT/TOF"); + /* load libraries */ + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + /* build analysis task class */ + gROOT->LoadMacro("AliAnalysisTaskPIDFluctuation.cxx+g"); + + /* setup input chain */ + TString str = inputfilename; + const Char_t *filename; + TChain *chain = new TChain("aodTree"); + if (str.EndsWith(".xml")) { + TGrid::Connect("alien://"); + Info("", "reading data list from collection:"); + TAlienCollection coll(inputfilename, maxFiles); + coll.Reset(); + while (coll.Next()) { + filename = coll.GetTURL(); + Info("", Form("%s", filename)); + chain->Add(filename); + } + } + else if (str.EndsWith(".txt")) { + Info("", "reading data list from text file:"); + ifstream is(inputfilename); + Char_t buf[4096]; + while(!is.eof()) { + is.getline(buf, 4096); + if (is.eof()) break; + chain->Add(buf); + Info("", Form("%s", buf)); + } + is.close(); + } + else { + Info("", "single file:"); + filename = inputfilename; + Info("", Form("%s", filename)); + chain->Add(filename); + } + Info("", Form("chain is ready: %d events", chain->GetEntries())); + + /* create analysis manager */ + AliAnalysisManager *mgr = new AliAnalysisManager("PIDFluctuation"); + + /* define input event handler */ + AliAODInputHandler *aodh = new AliAODInputHandler(); + mgr->SetInputEventHandler(aodh); + + gROOT->LoadMacro("AddAnalysisTaskPIDFluctuation.C"); + AddAnalysisTaskPIDFluctuation(AliAODTrack::kTrkGlobal, 0., 10., -0.8, 0.8); + AddAnalysisTaskPIDFluctuation(AliAODTrack::kTrkTPCOnly, 0., 10., -0.8, 0.8); + AddAnalysisTaskPIDFluctuation(AliAODTrack::kTrkTPCOnly, 0.3, 1.5, -0.8, 0.8); + AddAnalysisTaskPIDFluctuation(AliAODTrack::kTrkTPCOnly, 0.3, 1.5, -0.4, 0.4); + + /* start analysis */ + mgr->SetDebugLevel(0); + if (!mgr->InitAnalysis()) return; + mgr->PrintStatus(); + mgr->StartAnalysis("local", chain, maxEv); + + /* create dummy file to tell we are done */ + gSystem->Exec("touch done"); + +} diff --git a/PWGCF/EBYE/PIDFluctuation/macros/SteerAnalysisTaskPIDFluctuation.C b/PWGCF/EBYE/PIDFluctuation/macros/SteerAnalysisTaskPIDFluctuation.C new file mode 100644 index 00000000000..b08b61b83cf --- /dev/null +++ b/PWGCF/EBYE/PIDFluctuation/macros/SteerAnalysisTaskPIDFluctuation.C @@ -0,0 +1,75 @@ +SteerAnalysisTaskPIDFluctuation(const Char_t *inputfilename, Int_t maxFiles = kMaxInt, Int_t maxEv = kMaxInt) +{ + + /* include path for ACLic */ + gSystem->AddIncludePath("-I$ALICE_ROOT/include"); + gSystem->AddIncludePath("-I$ALICE_ROOT/TOF"); + /* load libraries */ + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + /* build analysis task class */ + gROOT->LoadMacro("AliAnalysisTaskPIDFluctuation.cxx+g"); + + /* setup input chain */ + TString str = inputfilename; + const Char_t *filename; + TChain *chain = new TChain("esdTree"); + if (str.EndsWith(".xml")) { + TGrid::Connect("alien://"); + Info("", "reading data list from collection:"); + TAlienCollection coll(inputfilename, maxFiles); + coll.Reset(); + while (coll.Next()) { + filename = coll.GetTURL(); + Info("", Form("%s", filename)); + chain->Add(filename); + } + } + else if (str.EndsWith(".txt")) { + Info("", "reading data list from text file:"); + ifstream is(inputfilename); + Char_t buf[4096]; + while(!is.eof()) { + is.getline(buf, 4096); + if (is.eof()) break; + chain->Add(buf); + Info("", Form("%s", buf)); + } + is.close(); + } + else { + Info("", "single file:"); + filename = inputfilename; + Info("", Form("%s", filename)); + chain->Add(filename); + } + Info("", Form("chain is ready: %d events", chain->GetEntries())); + + /* create analysis manager */ + AliAnalysisManager *mgr = new AliAnalysisManager("EbyEFluctuationPID"); + + /* define input event handler */ + AliESDInputHandler *esdh = new AliESDInputHandler(); + esdh->SetReadFriends(kFALSE); + mgr->SetInputEventHandler(esdh); + + /* add tasks */ + gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C"); + AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(kFALSE); + gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C"); + AliCentralitySelectionTask *centralityTask = AddTaskCentrality(); + gROOT->LoadMacro("AddAnalysisTaskPIDFluctuation.C"); + AliAnalysisTaskPIDFluctuation *thisTask = AddAnalysisTaskPIDFluctuation(); + + printf("ready-steady-go\n"); + + /* start analysis */ + mgr->SetDebugLevel(0); + if (!mgr->InitAnalysis()) return; + mgr->PrintStatus(); + mgr->StartAnalysis("local", chain, maxEv); + + /* create dummy file to tell we are done */ + gSystem->Exec("touch done"); + +} diff --git a/PWGCF/EBYE/PIDFluctuation/macros/SteerESDAnalysisTaskPIDFluctuation.C b/PWGCF/EBYE/PIDFluctuation/macros/SteerESDAnalysisTaskPIDFluctuation.C new file mode 100644 index 00000000000..3097c04b473 --- /dev/null +++ b/PWGCF/EBYE/PIDFluctuation/macros/SteerESDAnalysisTaskPIDFluctuation.C @@ -0,0 +1,79 @@ +SteerESDAnalysisTaskPIDFluctuation(const Char_t *inputfilename, Int_t maxFiles = kMaxInt, Int_t maxEv = kMaxInt) +{ + + /* include path for ACLic */ + gSystem->AddIncludePath("-I$ALICE_ROOT/include"); + gSystem->AddIncludePath("-I$ALICE_ROOT/TOF"); + /* load libraries */ + gSystem->Load("libANALYSIS"); + gSystem->Load("libANALYSISalice"); + /* build analysis task class */ + gROOT->LoadMacro("AliAnalysisTaskPIDFluctuation.cxx+g"); + + /* setup input chain */ + TString str = inputfilename; + const Char_t *filename; + TChain *chain = new TChain("esdTree"); + if (str.EndsWith(".xml")) { + TGrid::Connect("alien://"); + Info("", "reading data list from collection:"); + TAlienCollection coll(inputfilename, maxFiles); + coll.Reset(); + while (coll.Next()) { + filename = coll.GetTURL(); + Info("", Form("%s", filename)); + chain->Add(filename); + } + } + else if (str.EndsWith(".txt")) { + Info("", "reading data list from text file:"); + ifstream is(inputfilename); + Char_t buf[4096]; + while(!is.eof()) { + is.getline(buf, 4096); + if (is.eof()) break; + chain->Add(buf); + Info("", Form("%s", buf)); + } + is.close(); + } + else { + Info("", "single file:"); + filename = inputfilename; + Info("", Form("%s", filename)); + chain->Add(filename); + } + Info("", Form("chain is ready: %d events", chain->GetEntries())); + + /* create analysis manager */ + AliAnalysisManager *mgr = new AliAnalysisManager("EbyEFluctuationPID"); + + /* define input event handler */ + AliESDInputHandler *esdh = new AliESDInputHandler(); + esdh->SetReadFriends(kFALSE); + mgr->SetInputEventHandler(esdh); + + /* add tasks */ + gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C"); + AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(kFALSE); + gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskCentrality.C"); + AliCentralitySelectionTask *centralityTask = AddTaskCentrality(); + gROOT->LoadMacro("AddAnalysisTaskPIDFluctuation.C"); + + AddAnalysisTaskPIDFluctuation(AliAODTrack::kTrkGlobal, 0., 10., -0.8, 0.8); + AddAnalysisTaskPIDFluctuation(AliAODTrack::kTrkTPCOnly, 0., 10., -0.8, 0.8); + AddAnalysisTaskPIDFluctuation(AliAODTrack::kTrkTPCOnly, 0.3, 1.5, -0.8, 0.8); + AddAnalysisTaskPIDFluctuation(AliAODTrack::kTrkTPCOnly, 0.3, 1.5, -0.4, 0.4); + + printf("ready-steady-go\n"); + + /* start analysis */ + mgr->SetDebugLevel(0); + if (!mgr->InitAnalysis()) return; + mgr->PrintStatus(); + mgr->StartAnalysis("local", chain, maxEv); + + /* create dummy file to tell we are done */ + gSystem->Exec("touch done"); + +} diff --git a/PWGCF/EBYE/PIDFluctuation/task/AliAnalysisTaskPIDFluctuation.cxx b/PWGCF/EBYE/PIDFluctuation/task/AliAnalysisTaskPIDFluctuation.cxx new file mode 100644 index 00000000000..30fa9b9b695 --- /dev/null +++ b/PWGCF/EBYE/PIDFluctuation/task/AliAnalysisTaskPIDFluctuation.cxx @@ -0,0 +1,747 @@ +#include "AliAnalysisTaskPIDFluctuation.h" +#include "AliInputEventHandler.h" +#include "AliAnalysisManager.h" +#include "AliESDEvent.h" +#include "AliESDtrack.h" +#include "AliAODEvent.h" +#include "AliAODTrack.h" +#include "AliCentrality.h" +#include "AliPIDResponse.h" +#include "AliESDpid.h" +#include "AliAODpidUtil.h" +#include "TFile.h" +#include "TKey.h" +#include "TList.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TH3F.h" +#include "TCanvas.h" +#include "AliESDtrackCuts.h" + +/* + * Event by event PID fluctuation analysis + * author: Roberto Preghenella (R+) + * email: preghenella@bo.infn.it + * + */ + +ClassImp(AliAnalysisTaskPIDFluctuation) + +//________________________________________________________________________ + +const Char_t *AliAnalysisTaskPIDFluctuation::fgkEventCounterName[AliAnalysisTaskPIDFluctuation::kNEventCounters] = { + "AllEvents", + "PhysicsSelection", + "PrimayVertex", + "PrimayVertexSPD", + "VertexAccepted", + "GoodCentrality", + "AcceptedEvents" +}; + +const Char_t *AliAnalysisTaskPIDFluctuation::fgkEventCounterTitle[AliAnalysisTaskPIDFluctuation::kNEventCounters] = { + "all events", + "physics selection", + "primary vertex", + "primary vertex SPD", + "vertex accepted", + "good centrality", + "accepted events" +}; + +const Char_t *AliAnalysisTaskPIDFluctuation::fgkSparseDataName[AliAnalysisTaskPIDFluctuation::kNSparseData] = { + "centV0M", + "centTRK", + "Nch", + "Nplus", + "Nminus", + "Npi", + "Npiplus", + "Npiminus", + "Nka", + "Nkaplus", + "Nkaminus", + "Npr", + "Nprplus", + "Nprminus" +}; + +const Char_t *AliAnalysisTaskPIDFluctuation::fgkSparseDataTitle[AliAnalysisTaskPIDFluctuation::kNSparseData] = { + "centrality percentile (V0M)", + "centrality percentile (TRK)", + "N_{charged}", + "N_{+}", + "N_{-}", + "N_{#pi}", + "N_{#pi^{+}}", + "N_{#pi^{-}}", + "N_{K}", + "N_{K^{+}}", + "N_{K^{-}}", + "N_{p+#bar{p}}", + "N_{p}", + "N_{#bar{p}}" +}; + +//________________________________________________________________________ + +AliAnalysisTaskPIDFluctuation::AliAnalysisTaskPIDFluctuation(const Char_t *name) : + AliAnalysisTaskSE(name), + fPIDMethod(kTPCTOF), + fESDtrackCuts(NULL), + fAODfilterBit(AliAODTrack::kTrkGlobal), + fEtaMin(-0.8), + fEtaMax(0.8), + fPtMin(0.3), + fPtMax(1.5), + fPID(NULL), + fHistoList(NULL) +{ + + /* + * default constructor + */ + + DefineOutput(1, TList::Class()); +} + +//________________________________________________________________________ + +AliAnalysisTaskPIDFluctuation::~AliAnalysisTaskPIDFluctuation() +{ + + /* + * default destructor + */ + + if (fPID) delete fPID; + if (fHistoList) delete fHistoList; + +} + +//________________________________________________________________________ + +void AliAnalysisTaskPIDFluctuation::UserCreateOutputObjects() +{ + + /* + * user create output objects + */ + + fHistoList = new TList(); + fHistoList->SetOwner(kTRUE); + + fHistoEventCounter = new TH1F("hHistoEventCounter", "", kNEventCounters, 0, kNEventCounters); + for (Int_t ievc = 0; ievc < kNEventCounters; ievc++) + fHistoEventCounter->GetXaxis()->SetBinLabel(ievc + 1, fgkEventCounterTitle[ievc]); + fHistoList->Add(fHistoEventCounter); + + fHistoAcceptedTracks = new TH2F("hHistoAcceptedTracks", ";p_{T} (GeV/c)", 10, 0., 100., 50, 0., 5.); + fHistoList->Add(fHistoAcceptedTracks); + + fHistoTOFMatchedTracks = new TH2F("hHistoTOFMatchedTracks", ";p_{T} (GeV/c)", 10, 0., 100., 50, 0., 5.); + fHistoList->Add(fHistoTOFMatchedTracks); + + fHistoTPCdEdx = new TH3F("hHistoTPCdEdx", ";centrality percentile;p_{TPCin} (GeV/c);dE/dx (au.)", 10, 0., 100., 50, 0., 5., 500, 0., 500.); + fHistoList->Add(fHistoTPCdEdx); + + fHistoTPCdEdx_inclusive = new TH3F("hHistoTPCdEdx_inclusive", ";centrality percentile;p_{TPCin} (GeV/c);dE/dx (au.)", 10, 0., 100., 50, 0., 5., 500, 0., 500.); + fHistoList->Add(fHistoTPCdEdx_inclusive); + + fHistoTOFbeta = new TH3F(Form("hHistoTOFbeta"), ";centrality percentile;p (GeV/c);v/c", 10, 0., 100., 50, 0., 5., 500, 0.1, 1.1); + fHistoList->Add(fHistoTOFbeta); + + /* loop over species */ + for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { + + fHistoTPCdEdx_selected[ipart] = new TH3F(Form("hHistoTPCdEdx_selected_%s", AliPID::ParticleName(ipart)), ";centrality percentile;p_{TPCin} (GeV/c);dE/dx (au.)", 10, 0., 100., 50, 0., 5., 500, 0., 500.); + fHistoList->Add(fHistoTPCdEdx_selected[ipart]); + + fHistoTOFbeta_selected[ipart] = new TH3F(Form("hHistoTOFbeta_selected_%s", AliPID::ParticleName(ipart)), ";centrality percentile;p (GeV/c);v/c", 10, 0., 100., 50, 0., 5., 500, 0.1, 1.1); + fHistoList->Add(fHistoTOFbeta_selected[ipart]); + + fHistoNSigmaTPC[ipart] = new TH3F(Form("hHistoNSigmaTPC_%s", AliPID::ParticleName(ipart)), Form(";centrality percentile;p_{T} (GeV/c); N_{#sigma-TPC}^{%s}", AliPID::ParticleLatexName(ipart)), 10, 0., 100., 50, 0., 5., 200, -10., 10.); + fHistoList->Add(fHistoNSigmaTPC[ipart]); + + fHistoNSigmaTOF[ipart] = new TH3F(Form("hHistoNSigmaTOF_%s", AliPID::ParticleName(ipart)), Form(";centrality percentile;p_{T} (GeV/c); N_{#sigma-TOF}^{%s}", AliPID::ParticleLatexName(ipart)), 10, 0., 100., 50, 0., 5., 200, -10., 10.); + fHistoList->Add(fHistoNSigmaTOF[ipart]); + + } /* end of loop over species */ + + Int_t fgSparseDataBins[kNSparseData] = { + 20, + 20, + 5000, + 2500, + 2500, + 3000, + 1500, + 1500, + 1000, + 500, + 500, + 500, + 250, + 250 + }; + Double_t fgSparseDataMin[kNSparseData] = { + 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0. + }; + Double_t fgSparseDataMax[kNSparseData] = { + 100., 100., 5000., 2500., 2500., 3000., 1500., 1500., 1000., 500., 500., 500., 250., 250. + }; + + + fHistoCorrelation = new THnSparseI("hHistoCorrelation", "", kNSparseData, fgSparseDataBins, fgSparseDataMin, fgSparseDataMax); + for (Int_t iaxis = 0; iaxis < kNSparseData; iaxis++) + fHistoCorrelation->GetAxis(iaxis)->SetTitle(fgkSparseDataTitle[iaxis]); + fHistoList->Add(fHistoCorrelation); + + PostData(1, fHistoList); +} + +//________________________________________________________________________ + +void AliAnalysisTaskPIDFluctuation::UserExec(Option_t *) +{ + + /* + * user exec + */ + + /* get ESD event */ + AliVEvent *event = InputEvent(); + if (!event) return; + + /* accept event */ + if (!AcceptEvent(event)) return; + + /* get centrality object and centrality */ + AliCentrality *centrality = event->GetCentrality(); + Float_t cent_v0m = centrality->GetCentralityPercentileUnchecked("V0M"); + Float_t cent_trk = centrality->GetCentralityPercentileUnchecked("TRK"); + + /* init PID */ + InitPID(event); + + Bool_t pidFlag[AliPID::kSPECIES]; + Int_t icharge; + Int_t chargedCounts[2]; + Int_t pidCounts[AliPID::kSPECIES][2]; + for (Int_t i = 0; i < 2; i++) { + chargedCounts[i] = 0; + for (Int_t ii = 0; ii < 5; ii++) { + pidCounts[ii][i] = 0; + } + } + + /* loop over tracks */ + for (Int_t itrk = 0; itrk < event->GetNumberOfTracks(); itrk++) { + + /* get track */ + AliVTrack *track = (AliVTrack *)event->GetTrack(itrk); + if (!track) continue; + + /* accept track */ + if (!AcceptTrack(track)) continue; + + /* get charge */ + icharge = track->Charge() > 0 ? 0 : 1; + chargedCounts[icharge]++; + + /* make PID */ + MakePID(track, pidFlag, cent_v0m); + for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) + if (pidFlag[ipart]) + pidCounts[ipart][icharge]++; + + } + + /* fill histogram */ + Double_t vsparse[kNSparseData]; + vsparse[kCent_V0M] = cent_v0m; + vsparse[kCent_TRK] = cent_trk; + vsparse[kNch] = chargedCounts[0] + chargedCounts[1]; + vsparse[kNch_plus] = chargedCounts[0]; + vsparse[kNch_minus] = chargedCounts[1]; + vsparse[kNpi] = pidCounts[AliPID::kPion][0] + pidCounts[AliPID::kPion][1]; + vsparse[kNpi_plus] = pidCounts[AliPID::kPion][0]; + vsparse[kNpi_minus] = pidCounts[AliPID::kPion][1]; + vsparse[kNka] = pidCounts[AliPID::kKaon][0] + pidCounts[AliPID::kKaon][1]; + vsparse[kNka_plus] = pidCounts[AliPID::kKaon][0]; + vsparse[kNka_minus] = pidCounts[AliPID::kKaon][1]; + vsparse[kNpr] = pidCounts[AliPID::kProton][0] + pidCounts[AliPID::kProton][1]; + vsparse[kNpr_plus] = pidCounts[AliPID::kProton][0]; + vsparse[kNpr_minus] = pidCounts[AliPID::kProton][1]; + fHistoCorrelation->Fill(vsparse); + + PostData(1, fHistoList); +} + +//___________________________________________________________ + +Bool_t +AliAnalysisTaskPIDFluctuation::AcceptEvent(AliVEvent *event) const +{ + /* + * accept event + */ + + /* fill histo event counter */ + fHistoEventCounter->Fill(kAllEvents); + + /* check physics selection */ + if (!(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB)) return kFALSE; + fHistoEventCounter->Fill(kPhysicsSelection); + + /* check primary vertex */ + const AliVVertex *vertex = event->GetPrimaryVertex(); + if (vertex->GetNContributors() < 1) { + /* get ESD vertex SPD */ + if (event->InheritsFrom("AliESDEvent")) { + AliESDEvent *esdevent = dynamic_cast(event); + vertex = esdevent->GetPrimaryVertexSPD(); + } + /* get AOD vertex SPD */ + else if (event->InheritsFrom("AliAODEvent")) { + AliAODEvent *aodevent = dynamic_cast(event); + vertex = aodevent->GetPrimaryVertexSPD(); + } + if (vertex->GetNContributors() < 1) return kFALSE; + fHistoEventCounter->Fill(kPrimaryVertexSPD); + } + fHistoEventCounter->Fill(kPrimaryVertex); + + /* check vertex position */ + if (TMath::Abs(vertex->GetZ()) > 10.) return kFALSE; + fHistoEventCounter->Fill(kVertexAccepted); + + /* get centrality object and check quality */ + AliCentrality *centrality = event->GetCentrality(); + if (centrality->GetQuality() != 0) return kFALSE; + fHistoEventCounter->Fill(kGoodCentrality); + + /* event accepted */ + fHistoEventCounter->Fill(kAcceptedEvents); + return kTRUE; +} + +//___________________________________________________________ + +Bool_t +AliAnalysisTaskPIDFluctuation::AcceptTrack(AliVTrack *track) const +{ + /* + * accept track + */ + + /* check ESD track cuts */ + if (track->InheritsFrom("AliESDtrack")) { + AliESDtrack *esdtrack = dynamic_cast(track); + if (!fESDtrackCuts->AcceptTrack(esdtrack)) return kFALSE; + } + /* check AOD filter bit */ + else if (track->InheritsFrom("AliAODTrack")) { + AliAODTrack *aodtrack = dynamic_cast(track); + if (!aodtrack->TestFilterBit(fAODfilterBit)) return kFALSE; + } + + /* check eta range */ + if (track->Eta() < fEtaMin || + track->Eta() > fEtaMax) return kFALSE; + /* check pt range */ + if (track->Pt() < fPtMin || + track->Pt() > fPtMax) return kFALSE; + + /* track accepted */ + return kTRUE; +} + +//___________________________________________________________ + +Bool_t +AliAnalysisTaskPIDFluctuation::HasTPCPID(AliVTrack *track) const +{ + /* + * has TPC PID + */ + + /* check PID signal */ + if (track->GetTPCsignal() <= 0. || + track->GetTPCsignalN() == 0) return kFALSE; + return kTRUE; + +} + +//___________________________________________________________ + +Bool_t +AliAnalysisTaskPIDFluctuation::HasTOFPID(AliVTrack *track) const +{ + /* + * has TOF PID + */ + + /* check TOF matched track */ + if (!(track->GetStatus() & AliESDtrack::kTOFout)|| + !(track->GetStatus() & AliESDtrack::kTIME)) return kFALSE; + return kTRUE; + +} + +//___________________________________________________________ + +Double_t +AliAnalysisTaskPIDFluctuation::MakeTPCPID(AliVTrack *track, Double_t *nSigma) const +{ + /* + * make TPC PID + * returns measured dEdx if PID available, otherwise -1. + * fills nSigma array with TPC nsigmas for e, mu, pi, K, p + */ + + /* check TPC PID */ + if (!HasTPCPID(track)) return -1.; + + /* get TPC info */ + Double_t ptpc = track->GetTPCmomentum(); + Double_t dEdx = track->GetTPCsignal(); + Int_t dEdxN = track->GetTPCsignalN(); + + /* loop over particles */ + for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { + Double_t bethe = fPID->GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)ipart); + Double_t diff = dEdx - bethe; + Double_t sigma = fPID->GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)ipart); + nSigma[ipart] = diff / sigma; + } + return dEdx; + +} + +//___________________________________________________________ + +Double_t +AliAnalysisTaskPIDFluctuation::MakeTOFPID(AliVTrack *track, Double_t *nSigma) const +{ + /* + * make TOF PID + * returns measured beta if PID available, otherwise -1. + * fills nSigma array with TOF nsigmas for e, mu, pi, K, p + */ + + /* check TOF PID */ + if (!HasTOFPID(track)) return -1.; + + /* get TOF info */ + Double_t p = track->P(); + Double_t time = track->GetTOFsignal() - fPID->GetTOFResponse().GetStartTime(p); + Double_t timei[5]; + track->GetIntegratedTimes(timei); + + /* loop over particles */ + for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { + Double_t timez = time - timei[ipart]; + Double_t sigma = fPID->GetTOFResponse().GetExpectedSigma(p, timei[ipart], AliPID::ParticleMass(ipart)); + nSigma[ipart] = timez / sigma; + } + + return timei[0] / time; + +} + +//___________________________________________________________ + +void +AliAnalysisTaskPIDFluctuation::MakePID(AliVTrack *track, Bool_t *pidFlag, Float_t centrality) const +{ + /* + * make PID + * fills PID QA plots + * fills pidFlag array with PID flags for e, mu, pi, K, p + */ + + /* cut definitions + (better put them as static variables so they can be changed from outside) */ + Double_t fgTPCPIDmomcut[AliPID::kSPECIES] = {0., 0., 0.5, 0.5, 0.7}; + Double_t fgTPCPIDsigmacut[AliPID::kSPECIES] = {0., 0., 2., 2., 2.}; + Double_t fgTPCPIDcompcut[AliPID::kSPECIES] = {0., 0., 3., 3., 3.}; + Double_t fgTOFPIDmomcut[AliPID::kSPECIES] = {0., 0., 1.5, 1.5, 2.0}; + Double_t fgTOFPIDsigmacut[AliPID::kSPECIES] = {0., 0., 2., 2., 2.}; + + /* get momentum information */ + Double_t p = track->P(); + Double_t pt = track->Pt(); + Double_t ptpc = track->GetTPCmomentum(); + + /* make pid and check if available */ + Double_t nsigmaTPC[AliPID::kSPECIES]; + Double_t nsigmaTOF[AliPID::kSPECIES]; + Double_t dEdx = MakeTPCPID(track, nsigmaTPC); + Double_t beta = MakeTOFPID(track, nsigmaTOF); + Bool_t hasTPCPID = dEdx > 0.; + Bool_t hasTOFPID = beta > 0.; + + /* check PID method */ + if (fPIDMethod == kTPConly) hasTOFPID = kFALSE; // inhibit TOF PID + if (fPIDMethod == kTOFonly) hasTPCPID = kFALSE; // inhibit TPC PID + + /* fill qa histos */ + fHistoAcceptedTracks->Fill(centrality, pt); + if (hasTPCPID) + fHistoTPCdEdx->Fill(centrality, ptpc, dEdx); + if (hasTOFPID) { + fHistoTOFMatchedTracks->Fill(centrality, pt); + fHistoTOFbeta->Fill(centrality, p, beta); + } + if (hasTPCPID && !hasTOFPID) + fHistoTPCdEdx_inclusive->Fill(centrality, ptpc, dEdx); + + /* loop over species */ + for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) { + + /* reset PID flag */ + pidFlag[ipart] = kFALSE; + + /* fill qa histos */ + if (hasTPCPID) + fHistoNSigmaTPC[ipart]->Fill(centrality, pt, nsigmaTPC[ipart]); + if (hasTOFPID) + fHistoNSigmaTOF[ipart]->Fill(centrality, pt, nsigmaTOF[ipart]); + + /* combined PID cuts */ + if (hasTPCPID && hasTOFPID) { + if (pt < fgTOFPIDmomcut[ipart] && + TMath::Abs(nsigmaTOF[ipart]) < fgTOFPIDsigmacut[ipart] && + TMath::Abs(nsigmaTPC[ipart]) < fgTPCPIDcompcut[ipart]) { + fHistoTOFbeta_selected[ipart]->Fill(centrality, p, beta); + pidFlag[ipart] = kTRUE; + } + } + /* TPC-only PID cuts */ + else if (hasTPCPID && !hasTOFPID) { + if (pt < fgTPCPIDmomcut[ipart] && + TMath::Abs(nsigmaTPC[ipart]) < fgTPCPIDsigmacut[ipart]) { + fHistoTPCdEdx_selected[ipart]->Fill(centrality, ptpc, dEdx); + pidFlag[ipart] = kTRUE; + } + } + /* TOF-only PID cuts */ + else if (!hasTPCPID && hasTOFPID) { + if (pt < fgTOFPIDmomcut[ipart] && + TMath::Abs(nsigmaTOF[ipart]) < fgTOFPIDsigmacut[ipart]) { + fHistoTOFbeta_selected[ipart]->Fill(centrality, p, beta); + pidFlag[ipart] = kTRUE; + } + } + + } /* end of loop over species */ + +} + +//___________________________________________________________ + +void +AliAnalysisTaskPIDFluctuation::InitPID(AliVEvent *event) +{ + /* + * init PID + */ + + /* create PID object if not there yet */ + if (!fPID) { + + /* instance object */ + Bool_t mcFlag = kFALSE; /*** WARNING: check whether is MC ***/ + if (event->InheritsFrom("AliESDEvent")) + fPID = new AliESDpid(mcFlag); + else if (event->InheritsFrom("AliAODEvent")) + fPID = new AliAODpidUtil(mcFlag); + + /* set OADB path */ + fPID->SetOADBPath("$ALICE_ROOT/OADB"); + } + + /* init ESD PID */ + Int_t recoPass = 2; /*** WARNING: need to set the recoPass somehow better ***/ + fPID->InitialiseEvent(event, recoPass); /* warning: this apparently sets TOF time + * resolution to some hardcoded value, + * therefore we have to set correct + * resolution value after this call */ + + /* set TOF response */ + Double_t tofReso = 85.; /* ps */ /*** WARNING: need to set tofReso somehow better ***/ + fPID->GetTOFResponse().SetTimeResolution(tofReso); + fPID->GetTOFResponse().SetTrackParameter(0, 0.007); + fPID->GetTOFResponse().SetTrackParameter(1, 0.007); + fPID->GetTOFResponse().SetTrackParameter(2, 0.0); + fPID->GetTOFResponse().SetTrackParameter(3, 30); + +} + +//___________________________________________________________ + +void +AliAnalysisTaskPIDFluctuation::MeasureNuDyn(const Char_t *filename, Int_t i1, Int_t i2, Int_t centralityEstimator) +{ + + /* + * measure nu + */ + + printf("MeasureNuDyn: running for %s vs. %s\n", fgkSparseDataName[i1], fgkSparseDataName[i2]); + + /* get data */ + TFile *filein = TFile::Open(filename); + /* output */ + TFile *fileout = TFile::Open(Form("MeasureNuDyn_%s_%s.%s", fgkSparseDataName[i1], fgkSparseDataName[i2], filename), "RECREATE"); + + /* loop over available containers */ + TList *keylist = filein->GetListOfKeys(); + for (Int_t ikey = 0; ikey < keylist->GetEntries(); ikey++) { + + /* get key and check */ + TKey *key = (TKey *)keylist->At(ikey); + TString contname = key->GetName(); + if (!contname.BeginsWith("PIDFluctuation")) continue; + + /* get data */ + TList *list = (TList *)filein->Get(contname.Data()); + THnSparse *hsparse = (THnSparse *)list->FindObject("hHistoCorrelation"); + + /* create output directory and cd there */ + fileout->mkdir(contname.Data()); + fileout->cd(contname.Data()); + + /* loop over centralities */ + Double_t centBin[kNCentralityBins + 1] = { + 0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90. + }; + TH1D *hNu = new TH1D("hNu", ";centrality percentile;#nu", kNCentralityBins, centBin); + TH1D *hNuStat = new TH1D("hNuStat", ";centrality percentile;#nu_{stat}", kNCentralityBins, centBin); + TH1D *hNuDyn = new TH1D("hNuDyn", ";centrality percentile;#nu_{dyn}", kNCentralityBins, centBin); + for (Int_t icent = 0; icent < kNCentralityBins; icent++) { + + /* select centrality range */ + hsparse->GetAxis(centralityEstimator)->SetRangeUser(centBin[icent] + 0.001, centBin[icent + 1] - 0.001); + /* projections */ + TH1D *hcent = hsparse->Projection(centralityEstimator); + TH1D *h1 = hsparse->Projection(i1); + TH1D *h2 = hsparse->Projection(i2); + TH2D *hcorr = hsparse->Projection(i2, i1); + TH1D *hnu = new TH1D("hnu", "", 2000, -1., 1.); + + Double_t n, a, b, nev, meanx, rmsx, meany, rmsy, meannu, rmsnu; + Double_t meanx_err, meany_err, meannu_err; + + /* compute mean values */ + nev = 0.; meanx = 0.; meany = 0.; + for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++) + for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) { + n = hcorr->GetBinContent(ibinx + 1, ibiny + 1); + if (n <= 0.) continue; + meanx += n * ibinx; + meany += n * ibiny; + nev += n; + } + meanx /= nev; + meany /= nev; + // printf("nev = %f, meanx = %f, meany = %f\n", nev, meanx, meany); + + /* compute RMS values */ + nev = 0.; rmsx = 0.; rmsy = 0.; + for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++) + for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) { + n = hcorr->GetBinContent(ibinx + 1, ibiny + 1); + if (n <= 0.) continue; + a = ibinx - meanx; + rmsx += n * a * a; + a = ibiny - meany; + rmsy += n * a * a; + nev += n; + } + rmsx /= nev; + rmsx = TMath::Sqrt(rmsx); + rmsy /= nev; + rmsy = TMath::Sqrt(rmsy); + // printf("nev = %f, rmsx = %f, rmsy = %f\n", nev, rmsx, rmsy); + meanx_err = rmsx / TMath::Sqrt(nev); + meany_err = rmsy / TMath::Sqrt(nev); + // printf("nev = %f, meanx_err = %f, meany_err = %f\n", nev, meanx_err, meany_err); + + /* compute mean nu */ + nev = 0.; meannu = 0.; + for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++) + for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) { + n = hcorr->GetBinContent(ibinx + 1, ibiny + 1); + if (n <= 0.) continue; + a = ibinx / meanx - ibiny / meany; + meannu += n * a * a; + hnu->Fill(a * a, n); + nev += n; + } + meannu /= nev; + // printf("nev = %f, meannu = %f\n", nev, meannu); + + /* compute RMS nu */ + nev = 0.; rmsnu = 0.; + for (Int_t ibinx = 0; ibinx < hcorr->GetNbinsX(); ibinx++) + for (Int_t ibiny = 0; ibiny < hcorr->GetNbinsY(); ibiny++) { + n = hcorr->GetBinContent(ibinx + 1, ibiny + 1); + if (n <= 0.) continue; + a = ibinx / meanx - ibiny / meany; + b = a * a - meannu; + rmsnu += n * b * b; + nev += n; + } + rmsnu /= nev; + rmsnu = TMath::Sqrt(rmsnu); + // printf("nev = %f, rmsnu = %f\n", nev, rmsnu); + meannu_err = rmsnu / TMath::Sqrt(nev); + // printf("nev = %f, meannu_err = %f\n", nev, meannu_err); + + /* final calculations */ + Double_t nu = meannu; + Double_t nu_err = meannu_err; + Double_t nu_stat = 1. / meanx + 1. / meany; + Double_t meanx4 = meanx * meanx * meanx * meanx; + Double_t meanx_err2 = meanx_err * meanx_err; + Double_t meany4 = meany * meany * meany * meany; + Double_t meany_err2 = meany_err * meany_err; + Double_t nu_stat_err = TMath::Sqrt(meanx_err2 / meanx4 + meany_err2 / meany4); + Double_t nu_dyn = nu - nu_stat; + Double_t nu_dyn_err = TMath::Sqrt(nu_err * nu_err + nu_stat_err * nu_stat_err); + + /* setup final plots */ + hNu->SetBinContent(icent + 1, nu); + hNu->SetBinError(icent + 1, nu_err); + hNuStat->SetBinContent(icent + 1, nu_stat); + hNuStat->SetBinError(icent + 1, nu_stat_err); + hNuDyn->SetBinContent(icent + 1, nu_dyn); + hNuDyn->SetBinError(icent + 1, nu_dyn_err); + + hcent->Write(Form("hcent_cent%d", icent)); + h1->Write(Form("h1_cent%d", icent)); + h2->Write(Form("h2_cent%d", icent)); + hcorr->Write(Form("hcorr_cent%d", icent)); + hnu->Write(Form("hnu_cent%d", icent)); + + /* clean-up */ + delete hcent; + delete h1; + delete h2; + delete hcorr; + delete hnu; + } + + hNu->Write(); + hNuStat->Write(); + hNuDyn->Write(); + } + + fileout->Close(); + filein->Close(); + +} + +//___________________________________________________________ diff --git a/PWGCF/EBYE/PIDFluctuation/task/AliAnalysisTaskPIDFluctuation.h b/PWGCF/EBYE/PIDFluctuation/task/AliAnalysisTaskPIDFluctuation.h new file mode 100644 index 00000000000..7713edb52aa --- /dev/null +++ b/PWGCF/EBYE/PIDFluctuation/task/AliAnalysisTaskPIDFluctuation.h @@ -0,0 +1,136 @@ +#ifndef ALIANALYSISTASKPIDFLUCTUATION_H +#define ALIANALYSISTASKPIDFLUCTUATION_H + +/* + * Event by event PID fluctuation analysis + * author: Roberto Preghenella (R+) + * email: preghenella@bo.infn.it + * + */ + +class TList; +class TH1F; +class TH2F; +class TH3F; +class AliVEvent; +class AliVTrack; +class AliESDtrackCuts; +class AliPIDResponse; + +#include "AliAnalysisTaskSE.h" +#include "AliPID.h" +#include "THnSparse.h" + +class AliAnalysisTaskPIDFluctuation : +public AliAnalysisTaskSE +{ + + public: + + AliAnalysisTaskPIDFluctuation(const Char_t *name = "PIDFluctuation"); // default constructor + virtual ~AliAnalysisTaskPIDFluctuation(); // default destructor + + void UserCreateOutputObjects(); + void UserExec(Option_t *option); + + /* setters */ + void SetPIDMethod(Int_t value) {fPIDMethod = value;}; // set PID method + void SetESDtrackCuts(AliESDtrackCuts *value) {fESDtrackCuts = value;}; // set ESD track cuts + void SetAODfilterBit(Int_t value) {fAODfilterBit = value;}; // set AOD filter bit + void SetEtaRange(Float_t etaMin, Float_t etaMax) {fEtaMin = etaMin; fEtaMax = etaMax;}; // set eta range + void SetPtRange(Float_t ptMin, Float_t ptMax) {fPtMin = ptMin; fPtMax = ptMax;}; // set pt range + + static const Int_t kNCentralityBins = 10; // N centrality bins + + enum EEventCounter_t { + kAllEvents, + kPhysicsSelection, + kPrimaryVertex, + kPrimaryVertexSPD, + kVertexAccepted, + kGoodCentrality, + kAcceptedEvents, + kNEventCounters + }; + + enum ESparseData_t { + kCent_V0M, + kCent_TRK, + kNch, + kNch_plus, + kNch_minus, + kNpi, + kNpi_plus, + kNpi_minus, + kNka, + kNka_plus, + kNka_minus, + kNpr, + kNpr_plus, + kNpr_minus, + kNSparseData + }; + + enum EPIDMethod_t { + kTPCTOF, + kTPConly, + kTOFonly, + kNPIDMethods + }; + + static void MeasureNuDyn(const Char_t *filename, Int_t i1, Int_t i2, Int_t centralityEstimator = kCent_V0M); + + private: + + AliAnalysisTaskPIDFluctuation(const AliAnalysisTaskPIDFluctuation &); // not implemented + AliAnalysisTaskPIDFluctuation &operator=(const AliAnalysisTaskPIDFluctuation &); // not implemented + + /*** event and track selection ***/ + Bool_t AcceptEvent(AliVEvent *event) const; // accept event + Bool_t AcceptTrack(AliVTrack *track) const; // accept track + + /*** PID functions ***/ + Bool_t HasTPCPID(AliVTrack *track) const; // has TPC PID + Bool_t HasTOFPID(AliVTrack *track) const; // has TOF PID + Double_t MakeTPCPID(AliVTrack *track, Double_t *nSigma) const; // make TPC PID + Double_t MakeTOFPID(AliVTrack *track, Double_t *nSigma) const; // make TOF PID + void MakePID(AliVTrack *track, Bool_t *pidFlag, Float_t centrality) const; // make PID + void InitPID(AliVEvent *event); // init PID + + /*** PID objects and flags ***/ + Int_t fPIDMethod; // PID method + AliESDtrackCuts *fESDtrackCuts; // ESD track cuts + Int_t fAODfilterBit; // AOD filter bit + Float_t fEtaMin; // eta min + Float_t fEtaMax; // eta max + Float_t fPtMin; // pt min + Float_t fPtMax; // pt max + AliPIDResponse *fPID; // PID + + /*** PID histos ***/ + TList *fHistoList; // histo list + TH1F *fHistoEventCounter; // event counter + TH2F *fHistoAcceptedTracks; // accepted tracks + TH2F *fHistoTOFMatchedTracks; // TOF-matched tracks + TH3F *fHistoTPCdEdx; // TPC dEdx + TH3F *fHistoTPCdEdx_inclusive; // TPC dEdx + TH3F *fHistoTOFbeta; // TOF beta + TH3F *fHistoTPCdEdx_selected[AliPID::kSPECIES]; // TPC dEdx + TH3F *fHistoTOFbeta_selected[AliPID::kSPECIES]; // TOF beta + TH3F *fHistoNSigmaTPC[AliPID::kSPECIES]; // nsigma TPC + TH3F *fHistoNSigmaTOF[AliPID::kSPECIES]; // nsigma TOF + + /*** correlation histos */ + THnSparseI *fHistoCorrelation; // correlation THnSparse + + /*** labels, names and titles ***/ + static const Char_t *fgkEventCounterName[kNEventCounters]; // event couter name + static const Char_t *fgkEventCounterTitle[kNEventCounters]; // event couter title + static const Char_t *fgkSparseDataName[kNSparseData]; // sparse data name + static const Char_t *fgkSparseDataTitle[kNSparseData]; // sparse data title + + + ClassDef(AliAnalysisTaskPIDFluctuation, 1); +}; + +#endif /* ALIANALYSISTASKPIDFLUCTUATION_H */