From 03d2384659259e05da22af91368f6ba257efbb1d Mon Sep 17 00:00:00 2001 From: pulvir Date: Tue, 14 Jun 2011 15:56:26 +0000 Subject: [PATCH] Classes for 'mini' subpackage for RSN analysis framework --- PWG2/RESONANCES/AliRsnMiniAnalysisTask.cxx | 663 ++++++++++++++++++ PWG2/RESONANCES/AliRsnMiniAnalysisTask.h | 155 ++++ PWG2/RESONANCES/AliRsnMiniAxis.cxx | 70 ++ PWG2/RESONANCES/AliRsnMiniAxis.h | 46 ++ PWG2/RESONANCES/AliRsnMiniEvent.cxx | 38 + PWG2/RESONANCES/AliRsnMiniEvent.h | 42 ++ PWG2/RESONANCES/AliRsnMiniOutput.cxx | 506 +++++++++++++ PWG2/RESONANCES/AliRsnMiniOutput.h | 115 +++ PWG2/RESONANCES/AliRsnMiniPair.cxx | 69 ++ PWG2/RESONANCES/AliRsnMiniPair.h | 105 +++ PWG2/RESONANCES/AliRsnMiniParticle.cxx | 52 ++ PWG2/RESONANCES/AliRsnMiniParticle.h | 56 ++ PWG2/RESONANCES/AliRsnMiniValue.cxx | 171 +++++ PWG2/RESONANCES/AliRsnMiniValue.h | 75 ++ .../macros/mini/AddAnalysisTaskRsnMini.C | 90 +++ PWG2/RESONANCES/macros/mini/AddTaskTender.C | 69 ++ .../macros/mini/AnalysisSetupRsnMini.C | 146 ++++ PWG2/RESONANCES/macros/mini/ConfigPhi.C | 85 +++ PWG2/RESONANCES/macros/mini/ConfigPhiMC.C | 58 ++ PWG2/RESONANCES/macros/mini/runLocal.C | 161 +++++ PWG2/RESONANCES/macros/mini/testOut.C | 30 + 21 files changed, 2802 insertions(+) create mode 100644 PWG2/RESONANCES/AliRsnMiniAnalysisTask.cxx create mode 100644 PWG2/RESONANCES/AliRsnMiniAnalysisTask.h create mode 100644 PWG2/RESONANCES/AliRsnMiniAxis.cxx create mode 100644 PWG2/RESONANCES/AliRsnMiniAxis.h create mode 100644 PWG2/RESONANCES/AliRsnMiniEvent.cxx create mode 100644 PWG2/RESONANCES/AliRsnMiniEvent.h create mode 100644 PWG2/RESONANCES/AliRsnMiniOutput.cxx create mode 100644 PWG2/RESONANCES/AliRsnMiniOutput.h create mode 100644 PWG2/RESONANCES/AliRsnMiniPair.cxx create mode 100644 PWG2/RESONANCES/AliRsnMiniPair.h create mode 100644 PWG2/RESONANCES/AliRsnMiniParticle.cxx create mode 100644 PWG2/RESONANCES/AliRsnMiniParticle.h create mode 100644 PWG2/RESONANCES/AliRsnMiniValue.cxx create mode 100644 PWG2/RESONANCES/AliRsnMiniValue.h create mode 100644 PWG2/RESONANCES/macros/mini/AddAnalysisTaskRsnMini.C create mode 100644 PWG2/RESONANCES/macros/mini/AddTaskTender.C create mode 100644 PWG2/RESONANCES/macros/mini/AnalysisSetupRsnMini.C create mode 100644 PWG2/RESONANCES/macros/mini/ConfigPhi.C create mode 100644 PWG2/RESONANCES/macros/mini/ConfigPhiMC.C create mode 100644 PWG2/RESONANCES/macros/mini/runLocal.C create mode 100644 PWG2/RESONANCES/macros/mini/testOut.C diff --git a/PWG2/RESONANCES/AliRsnMiniAnalysisTask.cxx b/PWG2/RESONANCES/AliRsnMiniAnalysisTask.cxx new file mode 100644 index 00000000000..688792fe592 --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniAnalysisTask.cxx @@ -0,0 +1,663 @@ +// +// Analysis task for 'mini' sub-package +// Contains all definitions needed for running an analysis: +// -- global event cut +// -- a list of track cuts (any number) +// -- definitions of output histograms +// -- values to be computed. +// Each one must be defined using the "CREATE" methods, which +// add directly a new element in the task collections, and don't +// need an external object to be passed to the task itself. +// + +#include + +#include +#include +#include +#include +#include + +#include "AliLog.h" +#include "AliAnalysisManager.h" +#include "AliESDtrackCuts.h" +#include "AliESDUtils.h" +#include "AliMultiplicity.h" +#include "AliInputEventHandler.h" +#include "AliTriggerAnalysis.h" +#include "AliMCEvent.h" +#include "AliAODEvent.h" +#include "AliMCParticle.h" +#include "AliAODMCParticle.h" + +#include "AliRsnCutSet.h" +#include "AliRsnMiniEvent.h" +#include "AliRsnMiniParticle.h" +#include "AliRsnMiniPair.h" + +#include "AliRsnMiniAnalysisTask.h" + +ClassImp(AliRsnMiniAnalysisTask) + +//__________________________________________________________________________________________________ +AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask() : + AliAnalysisTaskSE(), + fEvNum(0), + fUseCentrality(kFALSE), + fCentralityType("QUALITY"), + fNMix(0), + fMaxDiffMult(10), + fMaxDiffVz(1.0), + fMaxDiffAngle(1E20), + fOutput(0x0), + fHistograms("AliRsnMiniOutput", 0), + fValues("AliRsnMiniValue", 0), + fHEventStat(0x0), + fEventCuts(0x0), + fTrackCuts(0), + fRsnEvent(), + fTempTree(0x0), + fNMixed(0), + fTriggerAna(0x0) +{ +// +// Dummy constructor ALWAYS needed for I/O. +// +} + +//__________________________________________________________________________________________________ +AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const char *name) : + AliAnalysisTaskSE(name), + fEvNum(0), + fUseCentrality(kFALSE), + fCentralityType("QUALITY"), + fNMix(0), + fMaxDiffMult(10), + fMaxDiffVz(1.0), + fMaxDiffAngle(1E20), + fOutput(0x0), + fHistograms("AliRsnMiniOutput", 0), + fValues("AliRsnMiniValue", 0), + fHEventStat(0x0), + fEventCuts(0x0), + fTrackCuts(0), + fRsnEvent(), + fTempTree(0x0), + fNMixed(0), + fTriggerAna(0x0) +{ +// +// Default constructor. +// Define input and output slots here (never in the dummy constructor) +// Input slot #0 works with a TChain - it is connected to the default input container +// Output slot #1 writes into a TH1 container +// + + DefineOutput(1, TList::Class()); +} + +//__________________________________________________________________________________________________ +AliRsnMiniAnalysisTask::AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask& copy) : + AliAnalysisTaskSE(copy), + fEvNum(0), + fUseCentrality(copy.fUseCentrality), + fCentralityType(copy.fCentralityType), + fNMix(copy.fNMix), + fMaxDiffMult(copy.fMaxDiffMult), + fMaxDiffVz(copy.fMaxDiffVz), + fMaxDiffAngle(copy.fMaxDiffAngle), + fOutput(0x0), + fHistograms(copy.fHistograms), + fValues(copy.fValues), + fHEventStat(0x0), + fEventCuts(copy.fEventCuts), + fTrackCuts(copy.fTrackCuts), + fRsnEvent(), + fTempTree(0x0), + fNMixed(0), + fTriggerAna(copy.fTriggerAna) +{ +// +// Copy constructor. +// Implemented as requested by C++ standards. +// Can be used in PROOF and by plugins. +// +} + +//__________________________________________________________________________________________________ +AliRsnMiniAnalysisTask& AliRsnMiniAnalysisTask::operator=(const AliRsnMiniAnalysisTask& copy) +{ +// +// Assignment operator. +// Implemented as requested by C++ standards. +// Can be used in PROOF and by plugins. +// + + AliAnalysisTaskSE::operator=(copy); + + fUseCentrality = copy.fUseCentrality; + fCentralityType = copy.fCentralityType; + fNMix = copy.fNMix; + fMaxDiffMult = copy.fMaxDiffMult; + fMaxDiffVz = copy.fMaxDiffVz; + fMaxDiffAngle = copy.fMaxDiffAngle; + fHistograms = copy.fHistograms; + fValues = copy.fValues; + fEventCuts = copy.fEventCuts; + fTrackCuts = copy.fTrackCuts; + fTriggerAna = copy.fTriggerAna; + + return (*this); +} + +//__________________________________________________________________________________________________ +AliRsnMiniAnalysisTask::~AliRsnMiniAnalysisTask() +{ +// +// Destructor. +// Clean-up the output list, but not the histograms that are put inside +// (the list is owner and will clean-up these histograms). Protect in PROOF case. +// + + if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { + delete fOutput; + } +} + +//__________________________________________________________________________________________________ +Int_t AliRsnMiniAnalysisTask::AddTrackCuts(AliRsnCutSet *cuts) +{ +// +// Add a new cut set for a new criterion for track selection. +// A user can add as many as he wants, and each one corresponds +// to one of the available bits in the AliRsnMiniParticle mask. +// The only check is the following: if a cut set with the same name +// as the argument is there, this is not added. +// Return value is the array position of this set. +// + + TObject *obj = fTrackCuts.FindObject(cuts->GetName()); + + if (obj) { + AliInfo(Form("A cut set named '%s' already exists", cuts->GetName())); + return fTrackCuts.IndexOf(obj); + } else { + fTrackCuts.AddLast(cuts); + return fTrackCuts.IndexOf(cuts); + } +} + +//__________________________________________________________________________________________________ +void AliRsnMiniAnalysisTask::UserCreateOutputObjects() +{ +// +// Initialization of outputs. +// This is called once per worker node. +// + + // reset counter + fEvNum = 0; + + // initialize trigger analysis + fTriggerAna = new AliTriggerAnalysis; + + // create list and set it as owner of its content (MANDATORY) + fOutput = new TList(); + fOutput->SetOwner(); + + // initialize event statistics counter + fHEventStat = new TH1F("hEventStat", "Event statistics", 4, 0.0, 4.0); + fHEventStat->GetXaxis()->SetBinLabel(1, "CINT1B"); + fHEventStat->GetXaxis()->SetBinLabel(2, "V0AND"); + fHEventStat->GetXaxis()->SetBinLabel(3, "Candle"); + fHEventStat->GetXaxis()->SetBinLabel(4, "Accepted"); + fOutput->Add(fHEventStat); + + // create temporary tree for filtered events + AliRsnMiniEvent *mini = 0x0; + fTempTree = new TTree("TempTree", "Temporary tree for events"); + fTempTree->Branch("events", "AliRsnMiniEvent", &mini); + + // create one histogram per each stored definition (event histograms) + Int_t i, ndef = fHistograms.GetEntries(); + AliRsnMiniOutput *def = 0x0; + for (i = 0; i < ndef; i++) { + def = (AliRsnMiniOutput*)fHistograms[i]; + if (!def) continue; + if (!def->Init(GetName(), fOutput)) { + AliError(Form("Def '%s': failed initialization", def->GetName())); + continue; + } + } + + // setup PID response + AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager(); + AliInputEventHandler *inputHandler = (AliInputEventHandler*)man->GetInputEventHandler(); + fRsnEvent.SetPIDResponse(inputHandler->GetPIDResponse()); + + // post data for ALL output slots >0 here, to get at least an empty histogram + PostData(1, fOutput); +} + +//__________________________________________________________________________________________________ +void AliRsnMiniAnalysisTask::UserExec(Option_t *) +{ +// +// Main computation loop. +// To be precise, this loop only checks if events are OK and eventually +// stores them in the temporary tree used for the analysis. +// The computation of correlations is done only at the end, +// in the automatic function which is called every time +// + + // event counter + fEvNum++; + + // check current event + Char_t check = CheckCurrentEvent(); + if (!check) return; + + // if the check is successful, the mini-event is created and stored + AliRsnMiniEvent *miniEvent = 0x0; + fTempTree->SetBranchAddress("events", &miniEvent); + + // assign event-related values + miniEvent = new AliRsnMiniEvent; + miniEvent->Vz() = fInputEvent->GetPrimaryVertex()->GetZ(); + miniEvent->Angle() = 0.0; + miniEvent->Mult() = ComputeCentrality((check == 'E')); + AliDebugClass(1, Form("Event %d: vz = %f -- mult = %f -- angle = %f", fEvNum, miniEvent->Vz(), miniEvent->Mult(), miniEvent->Angle())); + + // fill all histograms for events only + Int_t id, ndef = fHistograms.GetEntries(); + AliRsnMiniOutput *def = 0x0; + for (id = 0; id < ndef; id++) { + def = (AliRsnMiniOutput*)fHistograms[id]; + if (!def) continue; + if (!def->IsEventOnly()) continue; + def->Fill(miniEvent, &fValues); + } + + // fill all histograms for mother only (if MC is present) + if (fRsnEvent.IsESD() && fMCEvent) + FillTrueMotherESD(miniEvent); + else if (fRsnEvent.IsAOD() && fRsnEvent.GetAODList()) + FillTrueMotherAOD(miniEvent); + + // loop on daughters + // and store only those that pass at least one cut + Int_t it, ic, ncuts = fTrackCuts.GetEntries(), ntot = fRsnEvent.GetAbsoluteSum(); + AliRsnDaughter cursor; + AliRsnMiniParticle miniParticle; + for (it = 0; it < ntot; it++) { + fRsnEvent.SetDaughter(cursor, it); + miniParticle.CopyDaughter(&cursor); + for (ic = 0; ic < ncuts; ic++) { + AliRsnCutSet *cuts = (AliRsnCutSet*)fTrackCuts[ic]; + if (cuts->IsSelected(&cursor)) miniParticle.SetCutBit(ic); + } + AliDebugClass(3, Form("-- Track %d: REC px = %f -- py = %f -- pz = %f", it, miniParticle.Px(0), miniParticle.Py(0), miniParticle.Pz(0))); + AliDebugClass(3, Form("-- Track %d: SIM px = %f -- py = %f -- pz = %f", it, miniParticle.Px(1), miniParticle.Py(1), miniParticle.Pz(1))); + AliDebugClass(3, Form("-- Track %d: Charge = %c -- PDG = %d -- mother = %d -- motherPDG = %d", it, miniParticle.Charge(), miniParticle.PDG(), miniParticle.Mother(), miniParticle.MotherPDG())); + AliDebugClass(2, Form("-- Track %d: Cutbit = %u", it, miniParticle.CutBits())); + if (miniParticle.CutBits()) { + miniEvent->AddParticle(miniParticle); + } + } + AliDebugClass(1, Form("Event %d: selected tracks = %d", fEvNum, miniEvent->Particles().GetEntriesFast())); + + // store event + fTempTree->Fill(); + + // process single event + ProcessEvents(miniEvent); + + // post outputs + PostData(1, fOutput); +} + +//__________________________________________________________________________________________________ +void AliRsnMiniAnalysisTask::FinishTaskOutput() +{ +// +// At the end of execution, when the temporary tree is filled, +// perform mixing with all found events +// + + Int_t i1, i2, imix, nEvents = fTempTree->GetEntries(); + AliRsnMiniEvent *evMix = 0x0, evMain; + fTempTree->SetBranchAddress("events", &evMix); + + // initialize mixing counter + fNMixed.Set(nEvents); + TString msg(); + + // loop on events + for (i1 = 0; i1 < nEvents; i1++) { + if (fNMixed[i1] >= fNMix) continue; + fTempTree->GetEntry(i1); + evMain = (*evMix); + for (i2 = 1; i2 < nEvents; i2++) { + imix = i1 + i2; + if (imix >= nEvents) imix -= nEvents; + if (imix == i1) continue; + if (fNMixed[i1] >= fNMix) break; + if (fNMixed[imix] >= fNMix) continue; + fTempTree->GetEntry(imix); + // exit if events are not matched + if (TMath::Abs(evMain.Vz() - evMix->Vz()) > fMaxDiffVz) continue; + if (TMath::Abs(evMain.Mult() - evMix->Mult()) > fMaxDiffMult) continue; + if (TMath::Abs(evMain.Angle() - evMix->Angle()) > fMaxDiffAngle) continue; + // found a match: increment counter for both events + fNMixed[i1]++; + fNMixed[imix]++; + ProcessEvents(&evMain, evMix); + // if mixed enough times, stop + } + } + + // message + for (Int_t ii = 0; ii < nEvents; ii++) cout << Form("Event #%6d mixed %3d times", ii, fNMixed[ii]) << endl; + + // post outputs + PostData(1, fOutput); +} + +//__________________________________________________________________________________________________ +void AliRsnMiniAnalysisTask::Terminate(Option_t *) +{ +// +// Draw result to screen, or perform fitting, normalizations +// Called once at the end of the query +// + + fOutput = dynamic_cast(GetOutputData(1)); + if (!fOutput) { + AliError("Could not retrieve TList fOutput"); + return; + } +} + +//__________________________________________________________________________________________________ +Char_t AliRsnMiniAnalysisTask::CheckCurrentEvent() +{ +// +// This method fills the statistic histogram which counts the CINT1B, V0AND and CANDLE +// and checks current event against eventually defined local cuts. +// Its return values can be: +// -- 'E' if the event is accepted and is ESD +// -- 'A' if the event is accepted and is AOD +// -- 0 if the event is not accepted +// + + // string to sum messages + TString msg(""); + + // check input type + Char_t output = 0; + if (fInputEvent->InheritsFrom(AliESDEvent::Class())) + output = 'E'; + else if (fInputEvent->InheritsFrom(AliAODEvent::Class())) + output = 'A'; + else { + AliError(Form("Bad input event class: %s", fInputEvent->ClassName())); + return 0; + } + + // set reference to input + fRsnEvent.SetRef(fInputEvent); + + // assign MC event, if present + // for ESD it is the 'fMCEvent' data member in mother class + // for AOD it is the same event, but a check is done to look for the list of MC particles + // since there is an exit point above, if an event is not ESD here, it is surely AOD + if ((output == 'E') && fMCEvent) + fRsnEvent.SetRefMC(fMCEvent); + else { + fRsnEvent.SetRefMC(fInputEvent); + } + + // check physics selection + Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB); + if (isSelected) { + msg += "PHSEL = YES"; + fHEventStat->Fill(0.1); + } else { + AliDebugClass(1, "Event does not pass physics selections"); + return 0; + } + + // check if it is V0AND + Bool_t v0A = fTriggerAna->IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0A); + Bool_t v0C = fTriggerAna->IsOfflineTriggerFired((AliESDEvent*)fInputEvent, AliTriggerAnalysis::kV0C); + if (v0A && v0C) { + msg += " -- VOAND = YES"; + fHEventStat->Fill(1.1); + } else { + msg += " -- VOAND = NO "; + } + + // check candle + static AliESDtrackCuts *cuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(); + Int_t ntracksLoop = fInputEvent->GetNumberOfTracks(); + Bool_t candle = kFALSE; + for (Int_t iTrack = 0; iTrackGetTrack(iTrack); + AliESDtrack *esdt = dynamic_cast(track); + AliAODTrack *aodt = dynamic_cast(track); + if (esdt && !cuts->AcceptTrack(esdt)) continue; + if (aodt && !aodt->TestFilterBit(5)) continue; + if (track->Pt() < 0.5) continue; + if(TMath::Abs(track->Eta()) > 0.8) continue; + candle = kTRUE; + break; + } + if (candle) { + msg += " -- CANDLE = YES"; + fHEventStat->Fill(2.1); + } else { + msg += " -- CANDLE = NO "; + } + + // if event cuts are defined, they are checked here: + // an exit point is provided in case they are not passed + if (fEventCuts) { + if (!fEventCuts->IsSelected(&fRsnEvent)) { + msg += " -- Local cuts = REJECTED"; + AliDebugClass(1, Form("Stats for event %d: %s", fEvNum, msg.Data())); + return 0; + } else { + msg += " -- Local cuts = ACCEPTED"; + } + } else { + msg += " -- Local cuts = NONE"; + } + + // if the above exit point is not taken, the event is accepted + AliDebugClass(1, Form("Stats for event %d: %s", fEvNum, msg.Data())); + fHEventStat->Fill(3.1); + return output; +} + +//__________________________________________________________________________________________________ +Double_t AliRsnMiniAnalysisTask::ComputeCentrality(Bool_t isESD) +{ +// +// Computes event centrality/multiplicity according to the criterion defined +// by two elements: (1) choice between multiplicity and centrality and +// (2) the string defining what criterion must be used for specific computation. +// + + if (fUseCentrality) { + AliCentrality *centrality = fInputEvent->GetCentrality(); + if (!centrality) { + AliError("Cannot compute centrality!"); + return -1.0; + } + return centrality->GetCentralityPercentile(fCentralityType.Data()); + } else { + if (!isESD) { + AliInfo("Can compute only number of tracks for multiplicity"); + fCentralityType = "TRACKS"; + } + if (!fCentralityType.CompareTo("TRACKS")) + return fInputEvent->GetNumberOfTracks(); + else if (!fCentralityType.CompareTo("QUALITY")) + return AliESDtrackCuts::GetReferenceMultiplicity((AliESDEvent*)fInputEvent, kTRUE); + else if (!fCentralityType.CompareTo("TRACKLETS")) { + const AliMultiplicity *mult = ((AliESDEvent*)fInputEvent)->GetMultiplicity(); + Float_t nClusters[6] = {0.0,0.0,0.0,0.0,0.0,0.0}; + for(Int_t ilay = 0; ilay < 6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay); + return AliESDUtils::GetCorrSPD2(nClusters[1], fInputEvent->GetPrimaryVertex()->GetZ()); + } else { + AliError(Form("String '%s' does not define a possible multiplicity/centrality computation", fCentralityType.Data())); + return -1.0; + } + } +} + +//__________________________________________________________________________________________________ +void AliRsnMiniAnalysisTask::FillTrueMotherESD(AliRsnMiniEvent *miniEvent) +{ +// +// Fills the histograms with true mother (ESD version) +// + + Bool_t okMatch; + Int_t id, ndef = fHistograms.GetEntries(); + Int_t ip, label1, label2, npart = fMCEvent->GetNumberOfTracks(); + static AliRsnMiniPair miniPair; + AliMCParticle *daughter1, *daughter2; + TLorentzVector p1, p2; + AliRsnMiniOutput *def = 0x0; + + for (id = 0; id < ndef; id++) { + def = (AliRsnMiniOutput*)fHistograms[id]; + if (!def) continue; + if (!def->IsMother()) continue; + for (ip = 0; ip < npart; ip++) { + AliMCParticle *part = (AliMCParticle*)fMCEvent->GetTrack(ip); + if (TMath::Abs(part->Particle()->GetPdgCode()) != def->GetMotherPDG()) continue; + // check that daughters match expected species + if (part->Particle()->GetNDaughters() < 2) continue; + label1 = part->Particle()->GetDaughter(0); + label2 = part->Particle()->GetDaughter(1); + daughter1 = (AliMCParticle*)fMCEvent->GetTrack(label1); + daughter2 = (AliMCParticle*)fMCEvent->GetTrack(label2); + okMatch = kFALSE; + if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(1)) { + okMatch = kTRUE; + p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0)); + p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1)); + } else if (TMath::Abs(daughter1->Particle()->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->Particle()->GetPdgCode()) == def->GetPDG(0)) { + okMatch = kTRUE; + p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1)); + p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0)); + } + if (!okMatch) continue; + // assign momenta to computation object + miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2); + miniPair.FillRef(def->GetMotherMass()); + // do computations + def->Fill(&miniPair, miniEvent, &fValues); + } + } +} + +//__________________________________________________________________________________________________ +void AliRsnMiniAnalysisTask::FillTrueMotherAOD(AliRsnMiniEvent *miniEvent) +{ +// +// Fills the histograms with true mother (AOD version) +// + + Bool_t okMatch; + TClonesArray *list = fRsnEvent.GetAODList(); + Int_t id, ndef = fHistograms.GetEntries(); + Int_t ip, label1, label2, npart = list->GetEntries(); + static AliRsnMiniPair miniPair; + AliAODMCParticle *daughter1, *daughter2; + TLorentzVector p1, p2; + AliRsnMiniOutput *def = 0x0; + + for (id = 0; id < ndef; id++) { + def = (AliRsnMiniOutput*)fHistograms[id]; + if (!def) continue; + if (!def->IsMother()) continue; + for (ip = 0; ip < npart; ip++) { + AliAODMCParticle *part = (AliAODMCParticle*)list->At(ip); + if (TMath::Abs(part->GetPdgCode()) != def->GetMotherPDG()) continue; + // check that daughters match expected species + if (part->GetNDaughters() < 2) continue; + label1 = part->GetDaughter(0); + label2 = part->GetDaughter(1); + daughter1 = (AliAODMCParticle*)list->At(label1); + daughter2 = (AliAODMCParticle*)list->At(label2); + okMatch = kFALSE; + if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(0) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(1)) { + okMatch = kTRUE; + p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(0)); + p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(1)); + } else if (TMath::Abs(daughter1->GetPdgCode()) == def->GetPDG(1) && TMath::Abs(daughter2->GetPdgCode()) == def->GetPDG(0)) { + okMatch = kTRUE; + p1.SetXYZM(daughter1->Px(), daughter1->Py(), daughter1->Pz(), def->GetMass(1)); + p2.SetXYZM(daughter2->Px(), daughter2->Py(), daughter2->Pz(), def->GetMass(0)); + } + if (!okMatch) continue; + // assign momenta to computation object + miniPair.Sum(0) = miniPair.Sum(1) = (p1 + p2); + miniPair.FillRef(def->GetMotherMass()); + // do computations + def->Fill(&miniPair, miniEvent, &fValues); + } + } +} + +//________________________________________________________________________ +void AliRsnMiniAnalysisTask::ProcessEvents(AliRsnMiniEvent *evMain, AliRsnMiniEvent *evMix) +{ +// +// This function processes all single tracks. +// Loops on the daughters defined in current event +// and on all of them which pass the single-track cuts defined, computes what needed +// + + // check if it is mixed + Bool_t isMix; + if (evMix) { + isMix = kTRUE; + } else { + evMix = evMain; + isMix = kFALSE; + } + AliDebugClass(1, Form("Event %d: processing %s", fEvNum, (isMix ? "mixing" : "single event"))); + + // 2 nested loops on tracks + // inner starts from next track or from first, depending if is mixing or not + Int_t i1, i2, idef, start; + Int_t n1 = evMain->Particles().GetEntries(); + Int_t n2 = evMix ->Particles().GetEntries(); + Int_t ndef = fHistograms.GetEntries(); + AliRsnMiniParticle *daughter1 = 0x0, *daughter2 = 0x0; + AliRsnMiniOutput *def = 0x0; + for (i1 = 0; i1 < n1; i1++) { + daughter1 = (AliRsnMiniParticle*)evMain->Particles()[i1]; + if (!daughter1) continue; + if (isMix) start = 0; else start = i1 + 1; + for (i2 = start; i2 < n2; i2++) { + daughter2 = (AliRsnMiniParticle*)evMix->Particles()[i2]; + if (!daughter2) continue; + // loop on definitions + for (idef = 0; idef < ndef; idef++) { + def = (AliRsnMiniOutput*)fHistograms[idef]; + if (!def) continue; + if (def->IsEventOnly()) continue; + if (def->IsMother()) continue; + if (isMix && !def->IsTrackPairMix()) continue; + if (!isMix && def->IsTrackPairMix()) continue; + // check if definitions here match the current pair + def->Fill(daughter1, daughter2, evMain, &fValues); + } // end loop on definitions + } // end loop on daughter #2 + } // end loop on daughter #1 +} diff --git a/PWG2/RESONANCES/AliRsnMiniAnalysisTask.h b/PWG2/RESONANCES/AliRsnMiniAnalysisTask.h new file mode 100644 index 00000000000..5af424a47d2 --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniAnalysisTask.h @@ -0,0 +1,155 @@ +#ifndef ALIRSNMINIANALYSISTASK_H +#define ALIRSNMINIANALYSISTASK_H + +// +// Analysis task for 'mini' sub-package +// Contains all definitions needed for running an analysis: +// -- global event cut +// -- a list of track cuts (any number) +// -- definitions of output histograms +// -- values to be computed. +// + +#include "AliAnalysisTaskSE.h" + +#include +#include + +#include "AliRsnEvent.h" +#include "AliRsnMiniValue.h" +#include "AliRsnMiniOutput.h" + +class TList; +class AliTriggerAnalysis; +class AliRsnMiniEvent; +class AliRsnCutSet; + +class AliRsnMiniAnalysisTask : public AliAnalysisTaskSE { + +public: + + AliRsnMiniAnalysisTask(); + AliRsnMiniAnalysisTask(const char *name); + AliRsnMiniAnalysisTask(const AliRsnMiniAnalysisTask ©); + AliRsnMiniAnalysisTask& operator=(const AliRsnMiniAnalysisTask ©); + virtual ~AliRsnMiniAnalysisTask(); + + virtual void UserCreateOutputObjects(); + virtual void UserExec(Option_t *option); + virtual void Terminate(Option_t *); + virtual void FinishTaskOutput(); + + void SetNMix(Int_t nmix) {fNMix = nmix;} + void SetMaxDiffMult (Double_t val) {fMaxDiffMult = val;} + void SetMaxDiffVz (Double_t val) {fMaxDiffVz = val;} + void SetMaxDiffAngle(Double_t val) {fMaxDiffAngle = val;} + void SetEventCuts(AliRsnCutSet *cuts) {fEventCuts = cuts;} + void UseCentrality(const char *type) {fUseCentrality = kTRUE; fCentralityType = type; fCentralityType.ToUpper();} + void UseMultiplicity(const char *type) {fUseCentrality = kFALSE; fCentralityType = type; fCentralityType.ToUpper();} + Int_t AddTrackCuts(AliRsnCutSet *cuts); + + TClonesArray *Outputs() {return &fHistograms;} + TClonesArray *Values() {return &fValues;} + + Int_t ValueID(AliRsnMiniValue::EType type, Bool_t useMC = kFALSE); + Int_t CreateValue(AliRsnMiniValue::EType type, Bool_t useMC = kFALSE); + AliRsnMiniOutput *CreateOutput(const char *name, AliRsnMiniOutput::EOutputType type, AliRsnMiniOutput::EComputation src); + AliRsnMiniOutput *CreateOutput(const char *name, const char *outType, const char *compType); + + void ProcessEvents(AliRsnMiniEvent *evMain, AliRsnMiniEvent *evMix = 0x0); + +private: + + void InitHistograms(TClonesArray *array); + Char_t CheckCurrentEvent(); + Double_t ComputeCentrality(Bool_t isESD); + void FillTrueMotherESD(AliRsnMiniEvent *event); + void FillTrueMotherAOD(AliRsnMiniEvent *event); + void StoreTrueMother(AliRsnMiniPair *pair, AliRsnMiniEvent *event); + + Int_t fEvNum; //! absolute event counter + Bool_t fUseCentrality; // if true, use centrality for event, otherwise use multiplicity + TString fCentralityType; // definition used to choose what centrality or multiplicity to use + + Int_t fNMix; // mixing --> required number of mixes + Double_t fMaxDiffMult; // mixing --> max difference in multiplicity + Double_t fMaxDiffVz; // mixing --> max difference in Vz of prim vert + Double_t fMaxDiffAngle; // mixing --> max difference in reaction plane angle + + TList *fOutput; // output list + TClonesArray fHistograms; // list of histogram definitions + TClonesArray fValues; // list of values to be computed + TH1F *fHEventStat; // histogram of event statistics + + AliRsnCutSet *fEventCuts; // cuts on events + TObjArray fTrackCuts; // list of single track cuts + AliRsnEvent fRsnEvent; //! interface object to the event + TTree *fTempTree; //! tree to contain converted events (temporary) + TArrayI fNMixed; //! array to keep trace of how many times an event was mixed + AliTriggerAnalysis *fTriggerAna; //! trigger analysis + + ClassDef(AliRsnMiniAnalysisTask, 1); // AliRsnMiniAnalysisTask +}; + +inline Int_t AliRsnMiniAnalysisTask::CreateValue(AliRsnMiniValue::EType type, Bool_t useMC) +{ +// +// Create a new value in the task, +// and returns its ID, which is needed for setting up histograms. +// If that value was already initialized, returns its ID and does not recreate it. +// + + Int_t valID = ValueID(type, useMC); + if (valID >= 0 && valID < fValues.GetEntries()) { + AliInfo(Form("Value '%s' is already created in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID)); + } else { + valID = fValues.GetEntries(); + AliInfo(Form("Creating value '%s' in slot #%d", AliRsnMiniValue::ValueName(type, useMC), valID)); + new (fValues[valID]) AliRsnMiniValue(type, useMC); + } + + return valID; +} + +inline Int_t AliRsnMiniAnalysisTask::ValueID(AliRsnMiniValue::EType type, Bool_t useMC) +{ +// +// Searches if a value computation is initialized +// + + const char *name = AliRsnMiniValue::ValueName(type, useMC); + TObject *obj = fValues.FindObject(name); + if (obj) + return fValues.IndexOf(obj); + else + return -1; +} + +inline AliRsnMiniOutput* AliRsnMiniAnalysisTask::CreateOutput +(const char *name, AliRsnMiniOutput::EOutputType type, AliRsnMiniOutput::EComputation src) +{ +// +// Create a new histogram definition in the task, +// which is then returned to the user for its configuration +// + + Int_t n = fHistograms.GetEntries(); + AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, type, src); + + return newDef; +} + +inline AliRsnMiniOutput* AliRsnMiniAnalysisTask::CreateOutput +(const char *name, const char *outType, const char *compType) +{ +// +// Create a new histogram definition in the task, +// which is then returned to the user for its configuration +// + + Int_t n = fHistograms.GetEntries(); + AliRsnMiniOutput *newDef = new (fHistograms[n]) AliRsnMiniOutput(name, outType, compType); + + return newDef; +} +#endif diff --git a/PWG2/RESONANCES/AliRsnMiniAxis.cxx b/PWG2/RESONANCES/AliRsnMiniAxis.cxx new file mode 100644 index 00000000000..c3b692aec05 --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniAxis.cxx @@ -0,0 +1,70 @@ +// +// All implementations related to definition of an axis +// which is used in the output histogams. +// Simpler than TAxis, it defines an array of edges +// which is then ported to the output histogram definition. +// currently ported only in mini-package, but it could +// become a default also for general package. +// + +#include "TMath.h" +#include "AliRsnMiniAxis.h" + +ClassImp(AliRsnMiniAxis) + +//_____________________________________________________________________________ +void AliRsnMiniAxis::Set(Int_t nbins, Double_t min, Double_t max) +{ +// +// Set binning for the axis in equally spaced bins +// where the number of bins, minimum and maximum are given. +// + + if (!nbins) { + fBins.Set(0); + return; + } + + fBins.Set(nbins + 1); + + Double_t mymax = TMath::Max(min, max); + Double_t mymin = TMath::Min(min, max); + + Int_t k = 0; + Double_t binSize = (mymax - mymin) / ((Double_t)nbins); + + fBins[0] = mymin; + for (k = 1; k <= nbins; k++) fBins[k] = fBins[k - 1] + binSize; +} + +//_____________________________________________________________________________ +void AliRsnMiniAxis::Set(Double_t min, Double_t max, Double_t step) +{ +// +// Set binning for the axis in equally spaced bins +// where the bin size, minimum and maximum are given. +// + + Double_t dblNbins = TMath::Abs(max - min) / step; + Int_t intNbins = ((Int_t)dblNbins) + 1; + + Set(intNbins, min, max); +} + +//_____________________________________________________________________________ +void AliRsnMiniAxis::Set(Int_t nbins, Double_t *array) +{ +// +// Set binning for the axis in unequally spaced bins +// using the same way it is done in TAxis +// + + if (!nbins) { + fBins.Set(0); + return; + } + + Int_t i; + fBins.Set(nbins); + for (i = 0; i < nbins; i++) fBins[i] = array[i]; +} diff --git a/PWG2/RESONANCES/AliRsnMiniAxis.h b/PWG2/RESONANCES/AliRsnMiniAxis.h new file mode 100644 index 00000000000..40c739876fc --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniAxis.h @@ -0,0 +1,46 @@ +#ifndef ALIRSNMINIAXIS_H +#define ALIRSNMINIAXIS_H + +// +// All implementations related to definition of an axis +// which is used in the output histogams. +// Simpler than TAxis, it defines an array of edges +// which is then ported to the output histogram definition. +// currently ported only in mini-package, but it could +// become a default also for general package. +// + +#include "TObject.h" +#include "TArrayD.h" + +class AliRsnMiniAxis : public TObject { + +public: + + AliRsnMiniAxis(Int_t valID = -1) : fValueID(valID), fBins(0) { } + AliRsnMiniAxis(Int_t valID, Int_t nbins, Double_t min, Double_t max) : fValueID(valID), fBins(0) {Set(nbins, min, max);} + AliRsnMiniAxis(Int_t valID, Double_t min, Double_t max, Double_t step) : fValueID(valID), fBins(0) {Set(min, max, step);} + AliRsnMiniAxis(Int_t valID, Int_t nbins, Double_t *bins) : fValueID(valID), fBins(0) {Set(nbins, bins);} + AliRsnMiniAxis(const AliRsnMiniAxis& copy) : TObject(copy), fValueID(copy.fValueID), fBins(copy.fBins) { } + AliRsnMiniAxis& operator=(const AliRsnMiniAxis& copy) {fValueID = copy.fValueID; fBins = copy.fBins; return (*this);} + + void SetValueID(Int_t id) {fValueID = id;} + Int_t GetValueID() const {return fValueID;} + + Int_t NBins() {return fBins.GetSize() - 1;} + TArrayD *Bins() {return &fBins;} + Double_t *BinArray() {return fBins.GetArray();} + + void Set(Int_t nbins, Double_t min, Double_t max); + void Set(Int_t nbins, Double_t *bins); + void Set(Double_t min, Double_t max, Double_t step); + +private: + + Int_t fValueID; // index of used value in task collection + TArrayD fBins; // bins + + ClassDef(AliRsnMiniAxis,1) +}; + +#endif diff --git a/PWG2/RESONANCES/AliRsnMiniEvent.cxx b/PWG2/RESONANCES/AliRsnMiniEvent.cxx new file mode 100644 index 00000000000..3de40a132d5 --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniEvent.cxx @@ -0,0 +1,38 @@ +// +// Mini-Event +// Contains only the useful quantities computed on the event +// which can be used for event mixing, or for direct output +// when doing analysis w.r. to multiplicity or event plane, for example. +// +// Author: A. Pulvirenti +// + +#include "AliRsnMiniParticle.h" +#include "AliRsnMiniEvent.h" + +ClassImp(AliRsnMiniEvent) + +//__________________________________________________________________________________________________ +void AliRsnMiniEvent::AddParticle(AliRsnMiniParticle copy) +{ +// +// Add a new particle to the list and returns a pointer to it, +// in order to allow to se its parameters. +// + + Int_t n = fParticles.GetEntries(); + new (fParticles[n]) AliRsnMiniParticle(copy); +} + +//__________________________________________________________________________________________________ +AliRsnMiniParticle* AliRsnMiniEvent::LeadingParticle() +{ +// +// Return the leading particle +// + + if (fLeading < 0) return 0x0; + if (fLeading >= fParticles.GetEntriesFast()) return 0x0; + + return (AliRsnMiniParticle*)fParticles[fLeading]; +} diff --git a/PWG2/RESONANCES/AliRsnMiniEvent.h b/PWG2/RESONANCES/AliRsnMiniEvent.h new file mode 100644 index 00000000000..611fd75a4da --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniEvent.h @@ -0,0 +1,42 @@ +#ifndef ALIRSNMINIEVENT_H +#define ALIRSNMINIEVENT_H + +// +// Mini-Event +// Contains only the useful quantities computed on the event +// which can be used for event mixing, or for direct output +// when doing analysis w.r. to multiplicity or event plane, for example. +// + +#include +#include + +class AliRsnMiniParticle; + +class AliRsnMiniEvent : public TObject { +public: + + AliRsnMiniEvent() : fVz(0.0), fMult(0.0), fAngle(0.0), fLeading(-1), fParticles("AliRsnMiniParticle", 0) {} + ~AliRsnMiniEvent() {fParticles.Delete();} + + Float_t& Vz() {return fVz;} + Float_t& Mult() {return fMult;} + Float_t& Angle() {return fAngle;} + TClonesArray& Particles() {return fParticles;} + + AliRsnMiniParticle* LeadingParticle(); + void AddParticle(AliRsnMiniParticle copy); + +private: + + Float_t fVz; // z-position of vertex + Float_t fMult; // multiplicity or centrality + Float_t fAngle; // angle of reaction plane to main reference frame + + Int_t fLeading; // index of leading particle + TClonesArray fParticles; // list of selected particles + + ClassDef(AliRsnMiniEvent,1) +}; + +#endif diff --git a/PWG2/RESONANCES/AliRsnMiniOutput.cxx b/PWG2/RESONANCES/AliRsnMiniOutput.cxx new file mode 100644 index 00000000000..12cb7e302a8 --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniOutput.cxx @@ -0,0 +1,506 @@ +// +// Mini-Output +// All the definitions needed for building a RSN histogram +// including: +// -- properties of resonance (mass, PDG code if needed) +// -- properties of daughters (assigned mass, charges) +// -- definition of output histogram +// + +#include "Riostream.h" + +#include "TH1.h" +#include "TH2.h" +#include "TH3.h" +#include "TList.h" +#include "THnSparse.h" +#include "TString.h" +#include "TClonesArray.h" + +#include "AliRsnMiniParticle.h" +#include "AliRsnMiniPair.h" +#include "AliRsnMiniEvent.h" + +#include "AliLog.h" +#include "AliRsnCutSet.h" +#include "AliRsnMiniAxis.h" +#include "AliRsnMiniOutput.h" +#include "AliRsnMiniValue.h" + +ClassImp(AliRsnMiniOutput) + +//__________________________________________________________________________________________________ +AliRsnMiniOutput::AliRsnMiniOutput() : + TNamed(), + fOutputType(kTypes), + fComputation(kComputations), + fMotherPDG(0), + fMotherMass(0.0), + fPairCuts(0x0), + fOutputID(-1), + fAxes("AliRsnMiniAxis", 0), + fComputed(0), + fPair(), + fList(0x0) +{ +// +// Constructor +// + + fCutID[0] = fCutID[1] = -1; + fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown; + fCharge[0] = fCharge[1] = 0; +} + +//__________________________________________________________________________________________________ +AliRsnMiniOutput::AliRsnMiniOutput(const char *name, EOutputType type, EComputation src) : + TNamed(name, ""), + fOutputType(type), + fComputation(src), + fMotherPDG(0), + fMotherMass(0.0), + fPairCuts(0x0), + fOutputID(-1), + fAxes("AliRsnMiniAxis", 0), + fComputed(0), + fPair(), + fList(0x0) +{ +// +// Constructor +// + + fCutID[0] = fCutID[1] = -1; + fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown; + fCharge[0] = fCharge[1] = 0; +} + +//__________________________________________________________________________________________________ +AliRsnMiniOutput::AliRsnMiniOutput(const char *name, const char *outType, const char *compType) : + TNamed(name, ""), + fOutputType(kTypes), + fComputation(kComputations), + fMotherPDG(0), + fMotherMass(0.0), + fPairCuts(0x0), + fOutputID(-1), + fAxes("AliRsnMiniAxis", 0), + fComputed(0), + fPair(), + fList(0x0) +{ +// +// Constructor, with a more user friendly implementation, where +// the user sets the type of output and computations through conventional strings: +// +// Output: +// -- "HIST" --> common histogram (up to 3 dimensions) +// -- "SPARSE" --> sparse histogram +// +// Computation: +// -- "EVENT" --> event-only computations +// -- "PAIR" --> track pair computations (default) +// -- "MIX" --> event mixing (like track pair, but different events) +// -- "TRUE" --> true pairs (like track pair, but checking that come from same mother) +// -- "MOTHER" --> mother (loop on MC directly for mothers --> denominator of efficiency) +// + + TString input; + + // understand output type + input = outType; + input.ToUpper(); + if (!input.CompareTo("HIST")) + fOutputType = kHistogram; + else if (!input.CompareTo("SPARSE")) + fOutputType = kHistogramSparse; + else + AliWarning(Form("String '%s' does not define a meaningful output type", outType)); + + // understand computation type + input = compType; + input.ToUpper(); + if (!input.CompareTo("EVENT")) + fComputation = kEventOnly; + else if (!input.CompareTo("PAIR")) + fComputation = kTrackPair; + else if (!input.CompareTo("MIX")) + fComputation = kTrackPairMix; + else if (!input.CompareTo("TRUE")) + fComputation = kTruePair; + else if (!input.CompareTo("MOTHER")) + fComputation = kMother; + else + AliWarning(Form("String '%s' does not define a meaningful computation type", compType)); + + fCutID[0] = fCutID[1] = -1; + fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown; + fCharge[0] = fCharge[1] = 0; +} + +//__________________________________________________________________________________________________ +AliRsnMiniOutput::AliRsnMiniOutput(const AliRsnMiniOutput ©) : + TNamed(copy), + fOutputType(copy.fOutputType), + fComputation(copy.fComputation), + fMotherPDG(copy.fMotherPDG), + fMotherMass(copy.fMotherMass), + fPairCuts(copy.fPairCuts), + fOutputID(copy.fOutputID), + fAxes(copy.fAxes), + fComputed(copy.fComputed), + fPair(), + fList(copy.fList) +{ +// +// Copy constructor +// + + Int_t i; + for (i = 0; i < 2; i++) { + fCutID[i] = copy.fCutID[i]; + fDaughter[i] = copy.fDaughter[i]; + fCharge[i] = copy.fCharge[i]; + } +} + +//__________________________________________________________________________________________________ +AliRsnMiniOutput& AliRsnMiniOutput::operator=(const AliRsnMiniOutput ©) +{ +// +// Assignment operator +// + + fOutputType = copy.fOutputType; + fComputation = copy.fComputation; + fMotherPDG = copy.fMotherPDG; + fMotherMass = copy.fMotherMass; + fPairCuts = copy.fPairCuts; + fOutputID = copy.fOutputID; + fAxes = copy.fAxes; + fComputed = copy.fComputed; + fList = copy.fList; + + Int_t i; + for (i = 0; i < 2; i++) { + fCutID[i] = copy.fCutID[i]; + fDaughter[i] = copy.fDaughter[i]; + fCharge[i] = copy.fCharge[i]; + } + + return (*this); +} + + +//__________________________________________________________________________________________________ +void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t min, Double_t max) +{ +// +// Create a new axis reference +// + + Int_t size = fAxes.GetEntries(); + new (fAxes[size]) AliRsnMiniAxis(i, nbins, min, max); +} + +//__________________________________________________________________________________________________ +void AliRsnMiniOutput::AddAxis(Int_t i, Double_t min, Double_t max, Double_t step) +{ +// +// Create a new axis reference +// + + Int_t size = fAxes.GetEntries(); + new (fAxes[size]) AliRsnMiniAxis(i, min, max, step); +} + +//__________________________________________________________________________________________________ +void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t *values) +{ +// +// Create a new axis reference +// + + Int_t size = fAxes.GetEntries(); + new (fAxes[size]) AliRsnMiniAxis(i, nbins, values); +} + +//__________________________________________________________________________________________________ +Bool_t AliRsnMiniOutput::Init(const char *prefix, TList *list) +{ +// +// Initialize properly the histogram and add it to the argument list +// + + if (!list) { + AliError("Required an output list"); + return kFALSE; + } + + fList = list; + + switch (fOutputType) { + case kHistogram: + CreateHistogram(Form("hist_%s_%s", prefix, GetName())); + return kTRUE; + case kHistogramSparse: + CreateHistogramSparse(Form("hsparse_%s_%s", prefix, GetName())); + return kTRUE; + default: + AliError("Wrong output histogram definition"); + return kFALSE; + } +} + +//__________________________________________________________________________________________________ +void AliRsnMiniOutput::CreateHistogram(const char *name) +{ +// +// Initialize the 'default' TH1 output object. +// In case one of the expected axes is NULL, the initialization fails. +// + + Int_t size = fAxes.GetEntries(); + AliInfo(Form("Histogram name = '%s', with %d axes", name, size)); + + // we expect to have maximum 3 axes in this case + AliRsnMiniAxis *xAxis = 0x0, *yAxis = 0x0, *zAxis = 0x0; + if (size >= 1) xAxis = (AliRsnMiniAxis*)fAxes[0]; + if (size >= 2) yAxis = (AliRsnMiniAxis*)fAxes[1]; + if (size >= 3) zAxis = (AliRsnMiniAxis*)fAxes[2]; + + // create histogram depending on the number of axes + TH1 *h1 = 0x0; + if (xAxis && yAxis && zAxis) { + h1 = new TH3F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray(), zAxis->NBins(), zAxis->BinArray()); + } else if (xAxis && yAxis) { + h1 = new TH2F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray()); + } else if (xAxis) { + h1 = new TH1F(name, "", xAxis->NBins(), xAxis->BinArray()); + } else { + AliError("No axis was initialized"); + return; + } + + // switch the correct computation of errors + if (h1 && fList) { + h1->Sumw2(); + fList->Add(h1); + fOutputID = fList->IndexOf(h1); + } +} + +//________________________________________________________________________________________ +void AliRsnMiniOutput::CreateHistogramSparse(const char *name) +{ +// +// Initialize the THnSparse output object. +// In case one of the expected axes is NULL, the initialization fails. +// + + // retrieve binnings and sizes of all axes + // since the check for null values is done in Init(), + // we assume that here they must all be well defined + Int_t size = fAxes.GetEntries(); + Int_t i, *nbins = new Int_t[size]; + for (i = 0; i < size; i++) { + AliRsnMiniAxis *axis = (AliRsnMiniAxis*)fAxes[i]; + nbins[i] = axis->NBins(); + } + + // create fHSparseogram + THnSparseF *h1 = new THnSparseF(name, "", size, nbins); + + // update the various axes using the definitions given in the array of axes here + for (i = 0; i < size; i++) { + AliRsnMiniAxis *axis = (AliRsnMiniAxis*)fAxes[i]; + h1->GetAxis(i)->Set(nbins[i], axis->BinArray()); + } + + // clear heap + delete [] nbins; + + // add to list + if (h1 && fList) { + h1->Sumw2(); + fList->Add(h1); + fOutputID = fList->IndexOf(h1); + } +} + +//________________________________________________________________________________________ +Bool_t AliRsnMiniOutput::Fill(AliRsnMiniEvent *event, TClonesArray *valueList) +{ +// +// Compute values for event-based computations (does not use the pair) +// + + // check computation type + if (fComputation != kEventOnly) { + AliError("This method can be called only for event-based computations"); + return kFALSE; + } + + // get computed values + if (!ComputeValues(event, valueList)) return kFALSE; + + // call filling method + return FillHistogram(); +} + +//________________________________________________________________________________________ +Bool_t AliRsnMiniOutput::Fill(const AliRsnMiniPair *pair, AliRsnMiniEvent *event, TClonesArray *valueList) +{ +// +// Compute values for mother-based computations +// + + // check computation type + if (fComputation != kMother) { + AliError("This method can be called only for mother-based computations"); + return kFALSE; + } + + // copy passed pair info + fPair = (*pair); + + // check pair against cuts + if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE; + + // get computed values + if (!ComputeValues(event, valueList)) return kFALSE; + + // call filling method + return FillHistogram(); +} + +//________________________________________________________________________________________ +Bool_t AliRsnMiniOutput::Fill +(AliRsnMiniParticle *p1, AliRsnMiniParticle *p2, AliRsnMiniEvent *event, TClonesArray *valueList) +{ +// +// Compute values for pair-based comptuations +// + + // check computation type + if (fComputation != kTrackPair && fComputation != kTrackPairMix && fComputation != kTruePair) { + AliError("This method can be called only for pair-based computations"); + return kFALSE; + } + + // require that the two particles are well defined + if (!p1 || !p2) { + AliError("Required two particles"); + return kFALSE; + } + + // match check #1: first particle with first def, second particle with second def + // match check #2: first particle with second def, second particle with first def + if (p1->Charge() == fCharge[0] && p1->HasCutBit(fCutID[0]) && p2->Charge() == fCharge[1] && p2->HasCutBit(fCutID[1])) { + fPair.Fill(p1, p2, GetMass(0), GetMass(1), fMotherMass); + } else if (p1->Charge() == fCharge[1] && p1->HasCutBit(fCutID[1]) && p2->Charge() == fCharge[0] && p2->HasCutBit(fCutID[0])) { + fPair.Fill(p2, p1, GetMass(0), GetMass(1), fMotherMass); + } else { + AliDebugClass(1, "Pair does not match definition in any order"); + return kFALSE; + } + + // if required, check that this is a true pair + if (fComputation == kTruePair) { + //cout << fPair.Mother() << ' ' << fPair.MotherPDG() << ' '; + //cout << p1->PDG() << " (" << p1->Mother() << ") "; + //cout << p2->PDG() << " (" << p2->Mother() << ")" << endl; + if (fPair.Mother() < 0) { + return kFALSE; + } else if (TMath::Abs(fPair.MotherPDG()) != fMotherPDG) { + return kFALSE; + } + //cout << fPair.Mother() << ' ' << fPair.MotherPDG() << ' ' << p1->PDG() << ' ' << p2->PDG() << endl; + Bool_t decayMatch = kFALSE; + if (p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1])) + decayMatch = kTRUE; + if (p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1])) + decayMatch = kTRUE; + //cout << " DECAY MATCH = " << (decayMatch ? " OK " : " NO ") << endl; + if (!decayMatch) return kFALSE; + } + + // check pair against cuts + if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE; + + // get computed values + if (!ComputeValues(event, valueList)) return kFALSE; + + // call filling method + return FillHistogram(); +} + +//________________________________________________________________________________________ +Bool_t AliRsnMiniOutput::ComputeValues(AliRsnMiniEvent *event, TClonesArray *valueList) +{ +// +// Using the arguments and the internal 'fPair' data member, +// compute all values to be stored in the histogram +// + + // check size of computed array + Int_t size = fAxes.GetEntries(); + if (fComputed.GetSize() != size) fComputed.Set(size); + + Int_t i, ival, nval = valueList->GetEntries(); + + for (i = 0; i < size; i++) { + fComputed[i] = 1E20; + AliRsnMiniAxis *axis = (AliRsnMiniAxis*)fAxes[i]; + if (!axis) { + AliError("Null axis"); + continue; + } + ival = axis->GetValueID(); + if (ival < 0 || ival >= nval) { + AliError(Form("Required value #%d, while maximum is %d", ival, nval)); + continue; + } + AliRsnMiniValue *val = (AliRsnMiniValue*)valueList->At(ival); + if (!val) { + AliError(Form("Value in position #%d is NULL", ival)); + continue; + } + // if none of the above exit points is taken, compute value + fComputed[i] = val->Eval(&fPair, event); + } + + return kTRUE; +} + +//________________________________________________________________________________________ +Bool_t AliRsnMiniOutput::FillHistogram() +{ +// +// Fills the internal histogram using the current values stored in the +// 'fComputed' array, in the order as they are stored, up to the max +// dimension of the initialized histogram itself. +// + + // retrieve object from list + if (!fList) { + AliError("List pointer is NULL"); + return kFALSE; + } + TObject *obj = fList->At(fOutputID); + + if (obj->InheritsFrom(TH1F::Class())) { + ((TH1F*)obj)->Fill(fComputed[0]); + } else if (obj->InheritsFrom(TH2F::Class())) { + ((TH2F*)obj)->Fill(fComputed[0], fComputed[1]); + } else if (obj->InheritsFrom(TH3F::Class())) { + ((TH3F*)obj)->Fill(fComputed[0], fComputed[1], fComputed[2]); + } else if (obj->InheritsFrom(THnSparseF::Class())) { + ((THnSparseF*)obj)->Fill(fComputed.GetArray()); + } else { + AliError("No output initialized"); + return kFALSE; + } + + return kTRUE; +} diff --git a/PWG2/RESONANCES/AliRsnMiniOutput.h b/PWG2/RESONANCES/AliRsnMiniOutput.h new file mode 100644 index 00000000000..c3f94fa4039 --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniOutput.h @@ -0,0 +1,115 @@ +#ifndef ALIRSNMINIOUTPUT_H +#define ALIRSNMINIOUTPUT_H + +// +// Mini-Output +// All the definitions needed for building a RSN histogram +// including: +// -- properties of resonance (mass, PDG code if needed) +// -- properties of daughters (assigned mass, charges) +// -- definition of output histogram +// + +#include "AliRsnDaughter.h" +#include "AliRsnMiniParticle.h" + +class THnSparse; +class TList; +class TH1; + +class TList; +class TClonesArray; +class AliRsnMiniAxis; +class AliRsnMiniPair; +class AliRsnMiniEvent; + +typedef AliRsnDaughter::ESpecies RSNPID; + +class AliRsnMiniOutput : public TNamed { +public: + + enum EOutputType { + kHistogram, + kHistogramSparse, + kTypes + }; + + enum EComputation { + kEventOnly, + kTrackPair, + kTrackPairMix, + kTruePair, + kMother, + kComputations + }; + + AliRsnMiniOutput(); + AliRsnMiniOutput(const char *name, EOutputType type, EComputation src = kTrackPair); + AliRsnMiniOutput(const char *name, const char *outType, const char *compType); + AliRsnMiniOutput(const AliRsnMiniOutput ©); + AliRsnMiniOutput& operator=(const AliRsnMiniOutput ©); + + Bool_t IsEventOnly() const {return (fComputation == kEventOnly);} + Bool_t IsTrackPair() const {return (fComputation == kTrackPair);} + Bool_t IsTrackPairMix() const {return (fComputation == kTrackPairMix);} + Bool_t IsTruePair() const {return (fComputation == kTruePair);} + Bool_t IsMother() const {return (fComputation == kMother);} + Bool_t IsDefined() const {return (IsEventOnly() || IsTrackPair() || IsTrackPairMix() || IsTruePair() || IsMother());} + + EOutputType GetOutputType() const {return fOutputType;} + EComputation GetComputation() const {return fComputation;} + Int_t GetCutID(Int_t i) const {if (i <= 0) return fCutID [0]; else return fCutID [1];} + RSNPID GetDaughter(Int_t i) const {if (i <= 0) return fDaughter[0]; else return fDaughter[1];} + Double_t GetMass(Int_t i) const {return AliRsnDaughter::SpeciesMass(GetDaughter(i));} + Int_t GetPDG(Int_t i) const {return AliRsnDaughter::SpeciesPDG(GetDaughter(i));} + Int_t GetCharge(Int_t i) const {if (i <= 0) return fCharge[0]; else return fCharge[1];} + Int_t GetMotherPDG() const {return fMotherPDG;} + Double_t GetMotherMass() const {return fMotherMass;} + + void SetOutputType(EOutputType type) {fOutputType = type;} + void SetComputation(EComputation src) {fComputation = src;} + void SetCutID(Int_t i, Int_t value) {if (i <= 0) fCutID [0] = value; else fCutID [1] = value;} + void SetDaughter(Int_t i, RSNPID value) {if (i <= 0) fDaughter[0] = value; else fDaughter[1] = value;} + void SetCharge(Int_t i, Char_t value) {if (i <= 0) fCharge[0] = value; else fCharge[1] = value;} + void SetMotherPDG(Int_t pdg) {fMotherPDG = pdg;} + void SetMotherMass(Double_t mass) {fMotherMass = mass;} + void SetPairCuts(AliRsnCutSet *set) {fPairCuts = set;} + + void AddAxis(Int_t id, Int_t nbins, Double_t min, Double_t max); + void AddAxis(Int_t id, Double_t min, Double_t max, Double_t step); + void AddAxis(Int_t id, Int_t nbins, Double_t *values); + AliRsnMiniAxis* GetAxis(Int_t i) {if (i >= 0 && i < fAxes.GetEntries()) return (AliRsnMiniAxis*)fAxes[i]; return 0x0;} + Double_t* GetAllComputed() {return fComputed.GetArray();} + + AliRsnMiniPair& Pair() {return fPair;} + Bool_t Init(const char *prefix, TList *list); + Bool_t Fill(AliRsnMiniParticle *p1, AliRsnMiniParticle *p2, AliRsnMiniEvent *event, TClonesArray *valueList); + Bool_t Fill(const AliRsnMiniPair *pair, AliRsnMiniEvent *event, TClonesArray *valueList); + Bool_t Fill(AliRsnMiniEvent *event, TClonesArray *valueList); + +private: + + void CreateHistogram(const char *name); + void CreateHistogramSparse(const char *name); + Bool_t ComputeValues(AliRsnMiniEvent *event, TClonesArray *valueList); + Bool_t FillHistogram(); + + EOutputType fOutputType; // type of output + EComputation fComputation; // type of computation + Int_t fCutID[2]; // ID of cut set used to select tracks + RSNPID fDaughter[2]; // species of daughters + Char_t fCharge[2]; // required track charge + Int_t fMotherPDG; // PDG code of resonance + Double_t fMotherMass; // nominal resonance mass + AliRsnCutSet *fPairCuts; // cuts on the pair + + Int_t fOutputID; // index of output object in container list + TClonesArray fAxes; // definitions for the axes of each value + TArrayD fComputed; //! temporary container for all computed values + AliRsnMiniPair fPair; //! minipair for computations + TList *fList; //! pointer to the TList containing the output + + ClassDef(AliRsnMiniOutput,1) // AliRsnMiniOutput class +}; + +#endif diff --git a/PWG2/RESONANCES/AliRsnMiniPair.cxx b/PWG2/RESONANCES/AliRsnMiniPair.cxx new file mode 100644 index 00000000000..15025cee3f4 --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniPair.cxx @@ -0,0 +1,69 @@ +#include "AliRsnMiniParticle.h" +#include "AliRsnMiniPair.h" + +ClassImp(AliRsnMiniPair) + +//__________________________________________________________________________________________________ +void AliRsnMiniPair::Fill +(AliRsnMiniParticle *p1, AliRsnMiniParticle *p2, Double_t m1, Double_t m2, Double_t refMass) +{ +// +// Fill this object with data coming +// from arguments +// + + p1->Set4Vector(fP1[0], m1, kFALSE); + p2->Set4Vector(fP2[0], m2, kFALSE); + p1->Set4Vector(fP1[1], m1, kTRUE ); + p2->Set4Vector(fP2[1], m2, kTRUE ); + + fMother = -1; + if (p1->Mother() == p2->Mother()) { + fMother = p1->Mother(); + fMotherPDG = p1->MotherPDG(); + } + + Int_t i; + for (i = 0; i < 2; i++) { + fSum[i] = fP1[i] + fP2[i]; + fRef[i].SetXYZM(fSum[i].X(), fSum[i].Y(), fSum[i].Z(), refMass); + } +} + +//__________________________________________________________________________________________________ +Double_t AliRsnMiniPair::CosThetaStar(Bool_t useMC) +{ +// +// Return cosine of angle of one daughter to the resonance momentum in its rest frame +// + + TLorentzVector &mother = fSum[ID(useMC)]; + TLorentzVector &daughter0 = fP1[ID(useMC)]; + TLorentzVector &daughter1 = fP2[ID(useMC)]; + TVector3 momentumM(mother.Vect()); + TVector3 normal(mother.Y() / momentumM.Mag(), -mother.X() / momentumM.Mag(), 0.0); + + // Computes first the invariant mass of the mother + Double_t mass0 = daughter0.M(); + Double_t mass1 = daughter1.M(); + Double_t p0 = daughter0.Vect().Mag(); + Double_t p1 = daughter1.Vect().Mag(); + Double_t E0 = TMath::Sqrt(mass0 * mass0 + p0 * p0); + Double_t E1 = TMath::Sqrt(mass1 * mass1 + p1 * p1); + Double_t MotherMass = TMath::Sqrt((E0 + E1) * (E0 + E1) - (p0 * p0 + 2.0 * daughter0.Vect().Dot(daughter1.Vect()) + p1 * p1)); + MotherMass = mother.M(); + + // Computes components of beta + Double_t betaX = -mother.X() / mother.E(); + Double_t betaY = -mother.Y() / mother.E(); + Double_t betaZ = -mother.Z() / mother.E(); + + // Computes Lorentz transformation of the momentum of the first daughter + // into the rest frame of the mother and theta* + daughter0.Boost(betaX, betaY, betaZ); + TVector3 momentumD = daughter0.Vect(); + + Double_t cosThetaStar = normal.Dot(momentumD) / momentumD.Mag(); + + return cosThetaStar; +} diff --git a/PWG2/RESONANCES/AliRsnMiniPair.h b/PWG2/RESONANCES/AliRsnMiniPair.h new file mode 100644 index 00000000000..50cc35aeaa2 --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniPair.h @@ -0,0 +1,105 @@ +#ifndef ALIRSNMINIPAIR_H +#define ALIRSNMINIPAIR_H + +// +// This object is used as lightweight temporary container +// of all information needed from any input object and +// useful for resonance analysis. +// + +#include +#include + +class AliRsnMiniParticle; + +class AliRsnMiniPair : public TObject { +public: + + AliRsnMiniPair() : fMother(-1), fMotherPDG(0) { } + + Int_t& Mother() {return fMother;} + Int_t& MotherPDG() {return fMotherPDG;} + void Fill(AliRsnMiniParticle *p1, AliRsnMiniParticle *p2, Double_t m1, Double_t m2, Double_t refMass); + void FillRef(Double_t mass); + + Int_t ID(Bool_t mc) const {if (mc) return 1; else return 0;} + + TLorentzVector& P1 (Bool_t mc) {return fP1 [ID(mc)];} + TLorentzVector& P2 (Bool_t mc) {return fP2 [ID(mc)];} + TLorentzVector& Sum(Bool_t mc) {return fSum[ID(mc)];} + TLorentzVector& Ref(Bool_t mc) {return fRef[ID(mc)];} + + Double_t Pt(Bool_t mc) const {return fSum[ID(mc)].Pt();} + Double_t Pz(Bool_t mc) const {return fSum[ID(mc)].Pz();} + Double_t Eta(Bool_t mc) const {return fSum[ID(mc)].Eta();} + Double_t InvMass(Bool_t mc) const {return fSum[ID(mc)].M();} + Double_t InvMassRes() const; + Double_t Mt(Bool_t mc) const {return fRef[ID(mc)].Mt();} + Double_t Y(Bool_t mc) const {return fRef[ID(mc)].Rapidity();} + Double_t PtRatio(Bool_t mc) const; + Double_t DipAngle(Bool_t mc) const; + Double_t CosThetaStar(Bool_t mc); + +private: + + TLorentzVector fP1 [2]; // 1st daughter momentum + TLorentzVector fP2 [2]; // 2nd daughter momentum + TLorentzVector fSum[2]; // sum of momenta + TLorentzVector fRef[2]; // same as 'fSum' but with nominal resonance mass + + Int_t fMother; // label of mothers (when common) + Int_t fMotherPDG; // PDG code of mother (when common) + + ClassDef(AliRsnMiniPair,1) +}; + +inline void AliRsnMiniPair::FillRef(Double_t mass) +{ +// +// Fill ref 4-vectors using the passed mass and the values in 'sum' +// + + Int_t i; + for (i = 0; i < 2; i++) { + fRef[i].SetXYZM(fSum[i].X(), fSum[i].Y(), fSum[i].Z(), mass); + } +} + +inline Double_t AliRsnMiniPair::InvMassRes() const +{ +// +// Return invariant mass resolution +// + + if (fSum[1].M() <= 0.0) return 1E20; + + return (fSum[0].M() - fSum[1].M()) / fSum[1].M(); +} + +inline Double_t AliRsnMiniPair::PtRatio(Bool_t mc) const +{ +// +// Return ratio of transverse momenta of daughters +// + + Double_t num = TMath::Abs(fP1[ID(mc)].Perp() - fP2[ID(mc)].Perp()); + Double_t den = TMath::Abs(fP1[ID(mc)].Perp() + fP2[ID(mc)].Perp()); + + if (den <= 0.0) return 1E20; + + return num / den; +} + +inline Double_t AliRsnMiniPair::DipAngle(Bool_t mc) const +{ +// +// Opening angle in a Z-T space +// + + const TLorentzVector &p1 = fP1[ID(mc)]; + const TLorentzVector &p2 = fP2[ID(mc)]; + + return ((p1.Perp() * p2.Perp() + p1.Z() * p2.Z()) / p1.Mag() / p2.Mag()); +} + +#endif diff --git a/PWG2/RESONANCES/AliRsnMiniParticle.cxx b/PWG2/RESONANCES/AliRsnMiniParticle.cxx new file mode 100644 index 00000000000..5bdabb68400 --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniParticle.cxx @@ -0,0 +1,52 @@ +// +// This object is used as lightweight temporary container +// of all information needed from any input object and +// useful for resonance analysis. +// Lists of such objects are stored in a buffer, in order +// to allow an event mixing. +// + +#include "AliRsnDaughter.h" +#include "AliRsnMiniParticle.h" + +ClassImp(AliRsnMiniParticle) + +//__________________________________________________________________________________________________ +void AliRsnMiniParticle::CopyDaughter(AliRsnDaughter *daughter) +{ +// +// Sets data members from the passed object +// + + // reset what could not be initialized + fPDG = 0; + fMother = -1; + fMotherPDG = 0; + fCutBits = 0x0; + fPsim[0] = fPrec[0] = fPsim[1] = fPrec[1] = fPsim[2] = fPrec[2] = 0.0; + + // charge + if (daughter->IsPos()) + fCharge = '+'; + else if (daughter->IsNeg()) + fCharge = '-'; + else + fCharge = '0'; + + // rec info + if (daughter->GetRef()) { + fPrec[0] = daughter->GetRef()->Px(); + fPrec[1] = daughter->GetRef()->Py(); + fPrec[2] = daughter->GetRef()->Pz(); + } + + // MC info + if (daughter->GetRefMC()) { + fPsim[0] = daughter->GetRefMC()->Px(); + fPsim[1] = daughter->GetRefMC()->Py(); + fPsim[2] = daughter->GetRefMC()->Pz(); + fPDG = daughter->GetPDG(); + fMother = daughter->GetMother(); + fMotherPDG = daughter->GetMotherPDG(); + } +} diff --git a/PWG2/RESONANCES/AliRsnMiniParticle.h b/PWG2/RESONANCES/AliRsnMiniParticle.h new file mode 100644 index 00000000000..4c869fd65ac --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniParticle.h @@ -0,0 +1,56 @@ +#ifndef ALIRSNMINIPARTICLE_H +#define ALIRSNMINIPARTICLE_H + +// +// This object is used as lightweight temporary container +// of all information needed from any input object and +// useful for resonance analysis. +// + +#include +#include +#include + +class AliRsnDaughter; + +class AliRsnMiniParticle : public TObject { +public: + + AliRsnMiniParticle() : fCharge(0), fPDG(0), fMother(0), fMotherPDG(0), fCutBits(0x0) {Int_t i = 3; while (i--) fPsim[i] = fPrec[i] = 0.0;} + + Char_t& Charge() {return fCharge;} + Float_t& PsimX() {return fPsim[0];} + Float_t& PsimY() {return fPsim[1];} + Float_t& PsimZ() {return fPsim[2];} + Float_t& PrecX() {return fPrec[0];} + Float_t& PrecY() {return fPrec[1];} + Float_t& PrecZ() {return fPrec[2];} + Float_t& Px(Bool_t mc) {return (mc ? fPsim[0] : fPrec[0]);} + Float_t& Py(Bool_t mc) {return (mc ? fPsim[1] : fPrec[1]);} + Float_t& Pz(Bool_t mc) {return (mc ? fPsim[2] : fPrec[2]);} + Short_t& PDG() {return fPDG;} + Short_t PDGAbs() {return TMath::Abs(fPDG);} + Short_t& Mother() {return fMother;} + Short_t& MotherPDG() {return fMotherPDG;} + UShort_t& CutBits() {return fCutBits;} + Bool_t HasCutBit(Int_t i) {UShort_t bit = 1 << i; return ((fCutBits & bit) != 0);} + void SetCutBit(Int_t i) {UShort_t bit = 1 << i; fCutBits |= bit;} + void ClearCutBit(Int_t i) {UShort_t bit = 1 << i; fCutBits &= (~bit);} + + void Set4Vector(TLorentzVector &v, Float_t mass, Bool_t mc) {v.SetXYZM(Px(mc), Py(mc), Pz(mc), mass);} + void CopyDaughter(AliRsnDaughter *daughter); + +private: + + Char_t fCharge; // track charge *character*: '+', '-', '0' (whatever else = undefined) + Float_t fPsim[3]; // MC momentum of the track + Float_t fPrec[3]; // reconstructed momentum of the track + Short_t fPDG; // particle PDG code + Short_t fMother; // index of mother in its container + Short_t fMotherPDG; // PDG code of mother + UShort_t fCutBits; // list of bits used to know what cuts were passed by this track + + ClassDef(AliRsnMiniParticle,1) +}; + +#endif diff --git a/PWG2/RESONANCES/AliRsnMiniValue.cxx b/PWG2/RESONANCES/AliRsnMiniValue.cxx new file mode 100644 index 00000000000..894bf56f40e --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniValue.cxx @@ -0,0 +1,171 @@ +/************************************************************************** + * Copyright(c) 1998-1999, 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. * + **************************************************************************/ + +//////////////////////////////////////////////////////////////////////////////// +// +// This class contains all code which is used to compute any of the values +// which can be of interest within a resonance analysis. Besides the obvious +// invariant mass, it allows to compute other utility values on all possible +// targets, in order to allow a wide spectrum of binning and checks. +// When needed, this object can also define a binning in the variable which +// it is required to compute, which is used for initializing axes of output +// histograms (see AliRsnFunction). +// The value computation requires this object to be passed the object whose +// informations will be used. This object can be of any allowed input type +// (track, pair, event), then this class must inherit from AliRsnTarget. +// Then, when value computation is attempted, a check on target type is done +// and computation is successful only if expected target matches that of the +// passed object. +// In some cases, the value computation can require a support external object, +// which must then be passed to this class. It can be of any type inheriting +// from TObject. +// +// authors: A. Pulvirenti (alberto.pulvirenti@ct.infn.it) +// M. Vala (martin.vala@cern.ch) +// +//////////////////////////////////////////////////////////////////////////////// + +#include "Riostream.h" + +#include "AliLog.h" + +#include "AliRsnMiniPair.h" +#include "AliRsnMiniEvent.h" + +#include "AliRsnMiniValue.h" + +ClassImp(AliRsnMiniValue) + +//_____________________________________________________________________________ +AliRsnMiniValue::AliRsnMiniValue(EType type, Bool_t useMC) : + TNamed(ValueName(type, useMC), ""), + fType(type), + fUseMCInfo(useMC) +{ +// +// Constructor +// +} + +//_____________________________________________________________________________ +AliRsnMiniValue::AliRsnMiniValue(const AliRsnMiniValue& copy) : + TNamed(copy), + fType(copy.fType), + fUseMCInfo(copy.fUseMCInfo) +{ +// +// Copy constructor +// +} + +//_____________________________________________________________________________ +AliRsnMiniValue& AliRsnMiniValue::operator=(const AliRsnMiniValue& copy) +{ +// +// Assignment operator. +// Works like copy constructor. +// + + TNamed::operator=(copy); + fType = copy.fType; + fUseMCInfo = copy.fUseMCInfo; + + return (*this); +} + +//_____________________________________________________________________________ +const char* AliRsnMiniValue::TypeName(EType type) +{ +// +// This method returns a string to give a name to each possible +// computation value. +// + + switch (type) { + case kVz: return "EventVz"; + case kMult: return "EventMult"; + case kPlaneAngle: return "EventPlane"; + case kLeadingPt: return "EventLeadingPt"; + case kPt: return "Pt"; + case kPz: return "Pz"; + case kInvMass: return "InvMass"; + case kInvMassRes: return "InvMassResolution"; + case kEta: return "Eta"; + case kMt: return "Mt"; + case kY: return "Y"; + case kPtRatio: return "PtRatio"; + case kDipAngle: return "DipAngle"; + case kCosThetaStar: return "CosThetaStar"; + case kAngleLeading: return "AngleToLeading"; + default: return "Undefined"; + } +} + +//_____________________________________________________________________________ +Float_t AliRsnMiniValue::Eval(AliRsnMiniPair *pair, AliRsnMiniEvent *event) +{ +// +// Evaluation of the required value. +// In this implementation, fills the member 4-vectors with data +// coming from the object passed as argument, and then returns the value +// + + if (!pair && fType > kEventCuts) { + AliError("Null pair passed!"); + return 1E20; + } + + // compute value depending on types in the enumeration + // if the type does not match any available choice, or if + // the computation is not doable due to any problem + // (not initialized support object, wrong values, risk of floating point errors) + // the method returng kFALSE and sets the computed value to a meaningless number + switch (fType) { + // ---- event values ------------------------------------------------------------------------- + case kVz: + return event->Vz(); + case kMult: + return event->Mult(); + case kPlaneAngle: + return event->Angle(); + case kLeadingPt: + return 0.0; + // ---- pair values -------------------------------------------------------------------------- + case kPt: + return pair->Pt(fUseMCInfo); + case kInvMass: + return pair->InvMass(fUseMCInfo); + case kEta: + return pair->Eta(fUseMCInfo); + case kInvMassRes: + return pair->InvMassRes(); + case kMt: + return pair->Mt(fUseMCInfo); + case kY: + return pair->Y(fUseMCInfo); + case kPtRatio: + return pair->PtRatio(fUseMCInfo); + case kDipAngle: + return pair->DipAngle(fUseMCInfo); + case kCosThetaStar: + return pair->CosThetaStar(fUseMCInfo); + case kAngleLeading: + AliWarning("This method is not yet implemented"); + return 1E20; + default: + AliError("Invalid value type"); + return 1E20; + } +} diff --git a/PWG2/RESONANCES/AliRsnMiniValue.h b/PWG2/RESONANCES/AliRsnMiniValue.h new file mode 100644 index 00000000000..01849384191 --- /dev/null +++ b/PWG2/RESONANCES/AliRsnMiniValue.h @@ -0,0 +1,75 @@ +#ifndef ALIRSNMINIVALUE_H +#define ALIRSNMINIVALUE_H + +/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * See cxx source for full Copyright notice */ + +//////////////////////////////////////////////////////////////////////////////// +// +// Values which depend on 4-momentum of the pair. +// +//////////////////////////////////////////////////////////////////////////////// + +class AliRsnMiniPair; +class AliRsnMiniEvent; + +class AliRsnMiniValue : public TNamed { +public: + + enum EType { + kVz, // event Z position of primary vertex + kMult, // event multiplicity or centrality (depends on task settings) + kPlaneAngle, // event reaction plane angle + kLeadingPt, // event leading particle momentum + kEventCuts, // -- limit of event cuts ---------------------------------------------------- + kPt, // pair transverse momentum + kPz, // pair longitudinal momentum + kInvMass, // pair invariant mass (with reconstructed momenta) + kInvMassRes, // pair invariant mass resolution + kEta, // pair pseudo-rapidity + kMt, // pair transverse mass (need a reference mass) + kY, // pair rapidity (need a reference mass) + kPtRatio, // ratio |pt1 - pt2|/(pt1 + pt2) of daughter transverse momenta + kDipAngle, // inverse cosine of the angle between daughter vector momenta + kCosThetaStar, // polarization angle + kAngleLeading, // angle to leading particle + kTypes // -- general limit ---------------------------------------------------------- + }; + + AliRsnMiniValue(EType type = kTypes, Bool_t useMC = kFALSE); + AliRsnMiniValue(const AliRsnMiniValue& copy); + AliRsnMiniValue& operator=(const AliRsnMiniValue& copy); + virtual ~AliRsnMiniValue() { } + + void SetType(EType type) {fType = type;} + EType GetType() const {return fType;} + const char* GetTypeName() const {return TypeName(fType);} + Bool_t IsEventValue() const {return (fType < kEventCuts);} + + Float_t Eval(AliRsnMiniPair *pair, AliRsnMiniEvent *event = 0x0); + + static const char* TypeName(EType type); + static const char* ValueName(EType type, Bool_t useMC); + +protected: + + EType fType; // type from enumeration + Bool_t fUseMCInfo; // switch to use rec/sim momentum + + ClassDef(AliRsnMiniValue, 1) // AliRsnMiniValue class +}; + +inline const char* AliRsnMiniValue::ValueName(EType type, Bool_t useMC) +{ +// +// Define a criterion to name these object. +// They are not managed by the user, since each object is a singleton +// + + if (useMC) + return Form("MC_%s", TypeName(type)); + else + return TypeName(type); +} + +#endif diff --git a/PWG2/RESONANCES/macros/mini/AddAnalysisTaskRsnMini.C b/PWG2/RESONANCES/macros/mini/AddAnalysisTaskRsnMini.C new file mode 100644 index 00000000000..4ca8ec1a588 --- /dev/null +++ b/PWG2/RESONANCES/macros/mini/AddAnalysisTaskRsnMini.C @@ -0,0 +1,90 @@ +// +// General macro to configure the RSN analysis task. +// It calls all configs desired by the user, by means +// of the boolean switches defined in the first lines. +// --- +// Inputs: +// 1) flag to know if running on MC or data +// 2) path where all configs are stored +// --- +// Returns: +// kTRUE --> initialization successful +// kFALSE --> initialization failed (some config gave errors) +// + +Bool_t usePhi = 1; +Bool_t usePhiMC = 1; + +AliRsnMiniAnalysisTask * AddAnalysisTaskRsnMini +( + Bool_t isMC, + const char *path, + Int_t nmix = 0 +) +{ + // + // -- INITIALIZATION ---------------------------------------------------------------------------- + // + + // retrieve analysis manager + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + + // create the task and connect with physics selection + AliRsnMiniAnalysisTask *task = new AliRsnMiniAnalysisTask("RSN"); + mgr->AddTask(task); + + // settings + task->UseMultiplicity("TRACKS"); + + // set mixing + task->SetNMix(nmix); + task->SetMaxDiffVz(2.0); + task->SetMaxDiffMult(10.0); + task->SetMaxDiffAngle(1E20); + + // + // -- EVENT CUTS (same for all configs) --------------------------------------------------------- + // + + // cut on primary vertex: + // - 2nd argument --> |Vz| range + // - 3rd argument --> minimum required number of contributors + // - 4th argument --> tells if TPC stand-alone vertexes must be accepted + AliRsnCutPrimaryVertex *cutVertex = new AliRsnCutPrimaryVertex("cutVertex", 10.0, 0, kFALSE); + + // set the check for pileup + cutVertex->SetCheckPileUp(kTRUE); + + // define and fill cut set + AliRsnCutSet *eventCuts = new AliRsnCutSet("eventCuts", AliRsnTarget::kEvent); + eventCuts->AddCut(cutVertex); + eventCuts->SetCutScheme(cutVertex->GetName()); + + // set cuts in task + task->SetEventCuts(eventCuts); + + // + // -- CONFIGS ----------------------------------------------------------------------------------- + // + + if (usePhi) { + gROOT->LoadMacro(Form("%s/ConfigPhi.C", path)); + if (!ConfigPhi(task, isMC, "default")) return 0x0; + } + + if (isMC && usePhiMC) { + gROOT->LoadMacro(Form("%s/ConfigPhiMC.C", path)); + if (!ConfigPhiMC(task, "default")) return 0x0; + } + + // + // -- CONTAINERS -------------------------------------------------------------------------------- + // + + const char *file = AliAnalysisManager::GetCommonFileName(); + AliAnalysisDataContainer *output = mgr->CreateContainer("RsnOut", TList::Class(), AliAnalysisManager::kOutputContainer, file); + mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer()); + mgr->ConnectOutput(task, 1, output); + + return task; +} diff --git a/PWG2/RESONANCES/macros/mini/AddTaskTender.C b/PWG2/RESONANCES/macros/mini/AddTaskTender.C new file mode 100644 index 00000000000..df4f8a99eb3 --- /dev/null +++ b/PWG2/RESONANCES/macros/mini/AddTaskTender.C @@ -0,0 +1,69 @@ +AliAnalysisTask *AddTaskTender(Bool_t useV0 = kTRUE) +{ + // get the current analysis manager + Bool_t checkEvtSelection = useV0; + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) { + Error("AddTask_tender_Tender", "No analysis manager found."); + return 0; + } + // currently don't accept AOD input + //if (!mgr->GetInputEventHandler()->InheritsFrom(AliESDInputHandler::Class())) { + // Error("AddTask_tender_Tender","The analysis tender only works with ESD input!"); + // return 0; + // } + + + //========= Add tender to the ANALYSIS manager and set default storage ===== + AliTender *tender=new AliTender("AnalysisTender"); + tender->SetCheckEventSelection(checkEvtSelection); + tender->SetDefaultCDBStorage("raw://"); + mgr->AddTask(tender); + if (checkEvtSelection) { + if (mgr->GetTasks()->First() != (TObject*)tender) { + ::Error("When setting the tender to check the event selection, it has to be the first wagon ! Aborting."); + return NULL; + } + } + + //========= Attach VZERO supply ====== + if (useV0) { + AliVZEROTenderSupply *vzeroSupply=new AliVZEROTenderSupply("VZEROtender"); + vzeroSupply->SetDebug(kFALSE); + tender->AddSupply(vzeroSupply); + } + //========= Attach TPC supply ====== + AliTPCTenderSupply *tpcSupply=new AliTPCTenderSupply("TPCtender"); + tpcSupply->SetDebugLevel(2); + //tpcSupply->SetMip(50.); + tender->AddSupply(tpcSupply); + + //========= Attach TOF supply ====== + AliTOFTenderSupply *tofTender = new AliTOFTenderSupply("TOFtender"); + tender->AddSupply(tofTender); + + //========= Attach TRD supply ====== + AliTRDTenderSupply *trdSupply=new AliTRDTenderSupply("TRDtender"); + tender->AddSupply(trdSupply); + + //========= Attach PID supply ====== + tender->AddSupply(new AliPIDTenderSupply("PIDtender")); + + //========= Attach Primary Vertex supply ====== + tender->AddSupply(new AliVtxTenderSupply("PriVtxtender")); + + //================================================ + // data containers + //================================================ + + // define output containers, please use 'username'_'somename' + AliAnalysisDataContainer *coutput1 = + mgr->CreateContainer("tender_event", AliESDEvent::Class(), + AliAnalysisManager::kExchangeContainer,"default_tender"); + + // connect containers + mgr->ConnectInput (tender, 0, mgr->GetCommonInputContainer() ); + mgr->ConnectOutput (tender, 1, coutput1); + + return tender; +} diff --git a/PWG2/RESONANCES/macros/mini/AnalysisSetupRsnMini.C b/PWG2/RESONANCES/macros/mini/AnalysisSetupRsnMini.C new file mode 100644 index 00000000000..c58fde57631 --- /dev/null +++ b/PWG2/RESONANCES/macros/mini/AnalysisSetupRsnMini.C @@ -0,0 +1,146 @@ +// +// This macro sets all the aspects of configuration of an Analysis Train run +// which are always the same for all kinds of analysis (local, PROOF, AliEn) +// +// Inputs: +// +// - nmix = number of mixings to do (if > 0, initialize mixing stuff) +// - options = a set of keywords which drive some configurations +// - outputFileName = name of file produced by train +// - configPath = a path where all required config macros are stored +// +// Notes: +// +// - in case the source is an ESD, and if inputs are a MC production +// the MC input handler is created by default +// +// Returns: +// +// - if successful: the name of the expected input TTre (esdTree or aodTree) +// - if failed : NULL +// +TString Setup +( + Int_t nmix, + const char *options, + const char *outputFileName, + const char *macroPath = "." +) +{ + // prepare output + TString out(""); + + // + // === EXAMINE OPTIONS ========================================================================== + // + + // this is done using the utility 'RsnOptions.C' + // which provides a unique way to interpret them + + TString opt(options); + opt.ToUpper(); + + Bool_t isMC = opt.Contains("MC") || (!opt.Contains("DATA")); + Bool_t isESD = opt.Contains("ESD"); + Bool_t useTender = opt.Contains("TENDER"); + Bool_t noV0 = opt.Contains("NOV0"); + + // + // === LOAD LIBRARIES =========================================================================== + // + + // load analysis libraries + gSystem->Load("libSTEERBase.so"); + gSystem->Load("libESD.so"); + gSystem->Load("libAOD.so"); + gSystem->Load("libANALYSIS.so"); + gSystem->Load("libANALYSISalice.so"); + gSystem->Load("libEventMixing.so"); + gSystem->Load("libCORRFW.so"); + + // tender-related libraries + if (useTender) { + ::Info("AnalysisSetup", "Loading tender libraries"); + gSystem->Load("libTENDER.so"); + gSystem->Load("libTENDERSupplies.so"); + } + + // load development RSN library + if (!AliAnalysisAlien::SetupPar("PWG2resonances.par")) return ""; + + // + // === CREATE ANALYSIS MANAGER ================================================================== + // + + AliAnalysisManager *mgr = new AliAnalysisManager("RsnAnalysisManager"); + mgr->SetCommonFileName(outputFileName); + ::Info("AnalysisSetup", "Common file name: %s", outputFileName); + + // + // === INPUT / OUTPUT HANDLER CONFIGURATION ===================================================== + // + + // create multi input event handler + AliMultiInputEventHandler *multiHandler = new AliMultiInputEventHandler(); + + if (isESD) { + out = "esdTree"; + ::Info("AnalysisSetup", "Creating ESD handler"); + AliESDInputHandler *esdHandler = new AliESDInputHandler(); + mgr->SetInputEventHandler(esdHandler); + if (isMC) { + ::Info("AnalysisSetup", "Creating MC handler"); + AliMCEventHandler *mcHandler = new AliMCEventHandler(); + mgr->SetMCtruthEventHandler(mcHandler); + } + } else { + out = "aodTree"; + ::Info("AnalysisSetup", "Creating AOD handler"); + AliAODInputHandler *aodHandler = new AliAODInputHandler(); + mgr->SetInputEventHandler(aodHandler); + } + + // + // === TENDER TASK (ESD only -- optional) ======================================================= + // + + if (isESD && useTender) { + ::Info("AnalysisSetup", "Adding tender (and then accepting V0 info)", options); + gROOT->LoadMacro(Form("%s/AddTaskTender.C", macroPath)); + AddTaskTender(); + noV0 = kFALSE; + } + + // + // === PHYSICS SELECTION (ESD only) ============================================================= + // + + if (isESD) { + // setup physics selection + ::Info("AnalysisSetup", "Add physics selection by default on ESD analysis"); + gROOT->LoadMacro("$(ALICE_ROOT)/ANALYSIS/macros/AddTaskPhysicsSelection.C"); + AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(isMC); + if (noV0) { + ::Info("AnalysisSetup", "Skip of V0 info is required"); + physSelTask->GetPhysicsSelection()->SetSkipV0(kTRUE); + } + } + + // + // === PID RESPONSE ============================================================================= + // + + gROOT->LoadMacro("$(ALICE_ROOT)/ANALYSIS/macros/AddTaskPIDResponse.C"); + AddTaskPIDResponse(isMC); + + // + // === OTHER TASKS ============================================================================== + // + + // add RSN task + gROOT->LoadMacro(Form("%s/AddAnalysisTaskRsnMini.C", macroPath)); + if (!AddAnalysisTaskRsnMini(isMC, macroPath, nmix)) return kFALSE; + + ::Info("AnalysisSetup", "Setup successful"); + return out; +} diff --git a/PWG2/RESONANCES/macros/mini/ConfigPhi.C b/PWG2/RESONANCES/macros/mini/ConfigPhi.C new file mode 100644 index 00000000000..7918e423f07 --- /dev/null +++ b/PWG2/RESONANCES/macros/mini/ConfigPhi.C @@ -0,0 +1,85 @@ +// +// *** Configuration script for phi->KK analysis with 2010 runs *** +// +// A configuration script for RSN package needs to define the followings: +// +// (1) decay tree of each resonance to be studied, which is needed to select +// true pairs and to assign the right mass to all candidate daughters +// (2) cuts at all levels: single daughters, tracks, events +// (3) output objects: histograms or trees +// +Bool_t ConfigPhi(AliRsnMiniAnalysisTask *task, Bool_t isMC, const char *suffix) +{ + // + // -- Define track cuts ------------------------------------------------------------------------- + // + + // integrated kaon cut + AliRsnCutKaonForPhi2010PP *cutStd = new AliRsnCutKaonForPhi2010PP("cutStd"); + // cut set + AliRsnCutSet *cutSetStd = new AliRsnCutSet("kaonForPhi", AliRsnTarget::kDaughter); + cutSetStd->AddCut(cutStd); + cutSetStd->SetCutScheme(cutStd->GetName()); + // add to task + Int_t iCutStd = task->AddTrackCuts(cutSetStd); + + // + // -- Values ------------------------------------------------------------------------------------ + // + + // invariant mass + Int_t imID = task->CreateValue(AliRsnMiniValue::kInvMass, kFALSE); + + // transverse momentum + Int_t ptID = task->CreateValue(AliRsnMiniValue::kPt, kFALSE); + + // + // -- Pair cuts --------------------------------------------------------------------------------- + // + + AliRsnCutMiniPair *cutY = new AliRsnCutMiniPair("cutRapidity", AliRsnCutMiniPair::kRapidityRange); + cutY->SetRangeD(-0.5, 0.5); + + AliRsnCutSet *cutsPair = new AliRsnCutSet("pairCuts", AliRsnTarget::kMother); + cutsPair->AddCut(cutY); + cutsPair->SetCutScheme(cutY->GetName()); + + // + // -- Create all needed outputs ----------------------------------------------------------------- + // + + // use an array for more compact writing, which are different on mixing and charges + // [0] = unlike + // [1] = mixing + // [2] = like ++ + // [3] = like -- + Bool_t use [5] = { 1 , 1 , 1 , 1 , 1 }; + TString name [5] = {"Unlike", "Mixing", "LikePP", "LikeMM", "Trues" }; + TString comp [5] = {"PAIR" , "MIX" , "PAIR" , "PAIR" , "TRUE" }; + TString output [5] = {"HIST" , "HIST" , "HIST" , "HIST" , "HIST" }; + Char_t charge1 [5] = {'+' , '+' , '+' , '-' , '+' }; + Char_t charge2 [5] = {'-' , '-' , '+' , '-' , '-' }; + Int_t cutID [5] = { iCutStd, iCutStd, iCutStd, iCutStd, iCutStd }; + + for (Int_t i = 0; i < 5; i++) { + if (!use[i]) continue; + // create output + AliRsnMiniOutput *out = task->CreateOutput(Form("phi_%s_%s", name[i].Data(), suffix), output[i].Data(), comp[i].Data()); + // selection settings + out->SetCutID(0, cutID[i]); + out->SetCutID(1, cutID[i]); + out->SetDaughter(0, AliRsnDaughter::kKaon); + out->SetDaughter(1, AliRsnDaughter::kKaon); + out->SetCharge(0, charge1[i]); + out->SetCharge(1, charge2[i]); + out->SetMotherPDG(333); + out->SetMotherMass(1.019455); + // pair cuts + out->SetPairCuts(cutsPair); + // binnings + out->AddAxis(imID, 500, 0.9, 1.4); + out->AddAxis(ptID, 200, 0.0, 10.0); + } + + return kTRUE; +} diff --git a/PWG2/RESONANCES/macros/mini/ConfigPhiMC.C b/PWG2/RESONANCES/macros/mini/ConfigPhiMC.C new file mode 100644 index 00000000000..6ac78573c2c --- /dev/null +++ b/PWG2/RESONANCES/macros/mini/ConfigPhiMC.C @@ -0,0 +1,58 @@ +// +// *** Configuration script for phi->KK analysis with 2010 runs *** +// +// A configuration script for RSN package needs to define the followings: +// +// (1) decay tree of each resonance to be studied, which is needed to select +// true pairs and to assign the right mass to all candidate daughters +// (2) cuts at all levels: single daughters, tracks, events +// (3) output objects: histograms or trees +// +Bool_t ConfigPhiMC(AliRsnMiniAnalysisTask *task, const char *suffix) +{ + // + // -- Define track cuts ------------------------------------------------------------------------- + // + + /*** EMPTY FOR TRUE PAIRS COMPUTATION ***/ + + // + // -- Values ------------------------------------------------------------------------------------ + // + + // invariant mass + Int_t imID = task->CreateValue(AliRsnMiniValue::kInvMass, kTRUE); + + // transverse momentum + Int_t ptID = task->CreateValue(AliRsnMiniValue::kPt, kTRUE); + + // + // -- Pair cuts --------------------------------------------------------------------------------- + // + + AliRsnCutMiniPair *cutY = new AliRsnCutMiniPair("cutRapidityMC", AliRsnCutMiniPair::kRapidityRangeMC); + cutY->SetRangeD(-0.5, 0.5); + + AliRsnCutSet *cutsPair = new AliRsnCutSet("pairCutsMC", AliRsnTarget::kMother); + cutsPair->AddCut(cutY); + cutsPair->SetCutScheme(cutY->GetName()); + + // + // -- Create all needed outputs ----------------------------------------------------------------- + // + + // create output + AliRsnMiniOutput *out = task->CreateOutput(Form("phi_TrueMC_%s", suffix), "HIST", "MOTHER"); + // selection settings + out->SetDaughter(0, AliRsnDaughter::kKaon); + out->SetDaughter(1, AliRsnDaughter::kKaon); + out->SetMotherPDG(333); + out->SetMotherMass(1.019455); + // pair cuts + out->SetPairCuts(cutsPair); + // binnings + out->AddAxis(imID, 500, 0.9, 1.4); + out->AddAxis(ptID, 200, 0.0, 10.0); + + return kTRUE; +} diff --git a/PWG2/RESONANCES/macros/mini/runLocal.C b/PWG2/RESONANCES/macros/mini/runLocal.C new file mode 100644 index 00000000000..9579c46bc64 --- /dev/null +++ b/PWG2/RESONANCES/macros/mini/runLocal.C @@ -0,0 +1,161 @@ +void runLocal +( + Int_t nReadFiles = 0, + Int_t nSkipFiles = 0, + Int_t nmix = 10, + //const char *inputSource = "file-collections/AOD048_LHC11a10b.txt", + //const char *options = "aod", + //const char *inputSource = "file-collections/AOD049_LHC10h_pass2.txt", + //const char *options = "aod", + const char *inputSource = "file-collections/ESD_LHC10d1.txt", + const char *options = "esd_mc", + const char *outName = "test.root", + const char *macroPath = ".", + const char *setupName = "AnalysisSetupRsnMini.C" +) +{ + // + // === PREPARATION ============================================================================== + // + +// AliLog::SetClassDebugLevel("AliRsnMiniOutput" , 2); +// AliLog::SetClassDebugLevel("AliRsnMiniAnalysisTask", 1); + + // execute the general setup from the apposite macro + // it returns also a TString value with the input tree name + gROOT->LoadMacro(setupName); + TString out = Setup(nmix, options, outName, macroPath); + cout << "OUT = " << out.Data() << endl; + if (out.Length() < 1) return; + + // + // === BUILD INPUT LIST ========================================================================= + // + + // check extension of input to distinguish between XML and TXT + TString sInput(inputSource); + sInput.ToLower(); + Bool_t isTXT = (!strcmp(sInput(sInput.Length() - 3, 3).Data(), "txt")); + cout << "Input = " << (isTXT ? "TXT" : "XML") << endl; + + // if input is XML, connect to AliEn + if (!isTXT) TGrid::Connect("alien://"); + + // create TChain of input events + TChain *analysisChain = 0x0; + if (isTXT) analysisChain = CreateChainFromText(inputSource, out.Data(), nReadFiles, nSkipFiles); + else analysisChain = CreateChainFromXML (inputSource, out.Data(), nReadFiles, nSkipFiles); + if (!analysisChain) { + Error("runLocal", "Analysis chain not properly initialized"); + return; + } + + // + // === ANALYSIS TASK CREATION AND INCLUSION ===================================================== + // + + AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager(); + if (!mgr) return; + + // initialize and start analysis + if (!mgr->InitAnalysis()) { + ::Error("runLocal.C", "Failed to init analysis"); + return; + } + mgr->PrintStatus(); + if (isTXT) mgr->StartAnalysis("local", analysisChain); + else mgr->StartAnalysis("alien", analysisChain); +} + +//_________________________________________________________________________________________________ +TChain* CreateChainFromXML +(const char *xmlFileName, const char *treeName, Int_t nread, Int_t nskip) +{ +// +// Create a TChain with all required files listed into an XML collection. +// Necessary to run analysis in AliEn jobs. +// --- +// Arguments: +// - xmlFileName = input list +// - treeName = "esdTree" or "aodTree" +// - nread = how many files to read (0 = all) +// - nskip = how many files to skip from beginning +// + + // if nread argument is 0, it is disabled + if (nread == 0) nread = 1000000000; + + // initialize output object + TChain *chain = new TChain(treeName); + + // initialize the AliEn collection + TAlienCollection *myCollection = TAlienCollection::Open(xmlFileName); + if (!myCollection) { + Error("CreateChainFromXML", "Cannot create an AliEn collection from %s", xmlFileName); + return 0x0; + } + + // loop on collection + myCollection->Reset(); + while (myCollection->Next()) { + // skip until reached required number of offset + if (nskip > 0) {--nskip; continue;} + + // stop if required number of read files is reached + // otherwise update the counter + if (nread <= 0) break; + nread--; + + // recovery file and add it + Info("CreateChainFromXML", Form("Adding: %s", myCollection->GetTURL(""))); + chain->Add(myCollection->GetTURL("")); + } + + return chain; +} + +//__________________________________________________________________________________________________ +TChain* CreateChainFromText(const char *fileName, const char *treeName, Int_t nread, Int_t nskip) +{ +// +// Create a TChain with all required files listed into a text file. +// Necessary to run analysis in local jobs. +// --- +// Arguments: +// - xmlFileName = input file list +// - treeName = "esdTree" or "aodTree" +// - nread = how many files to read (0 = all) +// - nskip = how many files to skip from beginning +// + + // if third argument is 0, it is interpreted + // as "read all lines" + Bool_t readAll = (nread <= 0); + + // initialize output object + TChain* target = new TChain(treeName); + + // open text file + ifstream fileIn(fileName); + + // loop on collection + TString line; + while (fileIn.good()) { + fileIn >> line; + if (line.IsNull()) continue; + + // skip until reached required number of offset + if (nskip > 0) {--nskip; continue;} + + // stop if required number of read files is reached + // otherwise update the counter + if (!readAll && nread <= 0) break; + nread--; + + // add file + Info("CreateChainFromText", "Adding '%s'", line.Data()); + target->Add(line.Data()); + } + + return target; +} diff --git a/PWG2/RESONANCES/macros/mini/testOut.C b/PWG2/RESONANCES/macros/mini/testOut.C new file mode 100644 index 00000000000..d1cbf45acd1 --- /dev/null +++ b/PWG2/RESONANCES/macros/mini/testOut.C @@ -0,0 +1,30 @@ +void testOut(const char *file = "test.root", Double_t scale = 0.1, Int_t p1 = 0, Int_t p2 = -1) +{ + TFile *f = TFile::Open(file); + TList *l = (TList*)f->Get("RsnOut"); + + l->Print(); + + TH2F *hhPM = (TH2F*)RsnOut->FindObject("hist_RSN_phi_Unlike_default"); + TH2F *hhPP = (TH2F*)RsnOut->FindObject("hist_RSN_phi_LikePP_default"); + TH2F *hhMM = (TH2F*)RsnOut->FindObject("hist_RSN_phi_LikeMM_default"); + TH2F *hhMX = (TH2F*)RsnOut->FindObject("hist_RSN_phi_Mixing_default"); + + TH1D *hPM = 0x0; if (hhPM) hPM = hhPM->ProjectionX(Form("px1_%d_%d", p1, p2), p1, p2); + TH1D *hPP = 0x0; if (hhPP) hPP = hhPP->ProjectionX(Form("px2_%d_%d", p1, p2), p1, p2); + TH1D *hMM = 0x0; if (hhMM) hMM = hhMM->ProjectionX(Form("px3_%d_%d", p1, p2), p1, p2); + TH1D *hMX = 0x0; if (hhMX) hMX = hhMX->ProjectionX(Form("px4_%d_%d", p1, p2), p1, p2); + + hPM->SetLineColor(kBlack); + + hPP->SetLineColor(kGreen); + hPP->Add(hMM); + + hMX->SetLineColor(kRed); + hMX->Scale(scale); + + TCanvas *c1 = new TCanvas("c1"); + hPM->Draw(); + hPP->Draw("same"); + hMX->Draw("same"); +} -- 2.39.3