--- /dev/null
+//
+// 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 <Riostream.h>
+
+#include <TList.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TH3.h>
+#include <THnSparse.h>
+
+#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<TList*>(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; iTrack<ntracksLoop; iTrack++) {
+ AliVTrack *track = (AliVTrack*)fInputEvent->GetTrack(iTrack);
+ AliESDtrack *esdt = dynamic_cast<AliESDtrack*>(track);
+ AliAODTrack *aodt = dynamic_cast<AliAODTrack*>(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
+}
--- /dev/null
+#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 <TString.h>
+#include <TClonesArray.h>
+
+#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
--- /dev/null
+//
+// 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];
+}
--- /dev/null
+#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
--- /dev/null
+//
+// 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];
+}
--- /dev/null
+#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 <TArrayI.h>
+#include <TClonesArray.h>
+
+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
--- /dev/null
+//
+// 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;
+}
--- /dev/null
+#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
--- /dev/null
+#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;
+}
--- /dev/null
+#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 <TObject.h>
+#include <TLorentzVector.h>
+
+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
--- /dev/null
+//
+// 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();
+ }
+}
--- /dev/null
+#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 <TMath.h>
+#include <TObject.h>
+#include <TLorentzVector.h>
+
+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
--- /dev/null
+/**************************************************************************
+ * 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;
+ }
+}
--- /dev/null
+#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
--- /dev/null
+//
+// 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;
+}
--- /dev/null
+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;
+}
--- /dev/null
+//
+// 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;
+}
--- /dev/null
+//
+// *** 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;
+}
--- /dev/null
+//
+// *** 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;
+}
--- /dev/null
+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;
+}
--- /dev/null
+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");
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff