[u/mrichter/AliRoot.git] / PWGLF / RESONANCES / AliRsnMiniOutput.cxx

//
// 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 "TMath.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),
   fSel1(0),
   fSel2(0)
{
//
// 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),
   fSel1(0),
   fSel2(0)
{
//
// 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),
   fSel1(0),
   fSel2(0)
{
//
// 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)
//    -- "ROTATE1" --> rotated background (rotate first track)
//    -- "ROTATE2" --> rotated background (rotate second track)
//    -- "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("ROTATE1"))
      fComputation = kTrackPairRotated1;
   else if (!input.CompareTo("ROTATE2"))
      fComputation = kTrackPairRotated2;
   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 &copy) :
   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),
   fSel1(0),
   fSel2(0)
{
//
// 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 &copy)
{
//
// Assignment operator
//
   if (this == &copy)
      return *this;
   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];
   }

   fSel1.Set(0);
   fSel2.Set(0);

   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;
   Int_t size = fAxes.GetEntries();
   if (size < 1) {
      AliWarning(Form("[%s] Cannot initialize histogram with less than 1 axis", GetName()));
      return kFALSE;
   }

   switch (fOutputType) {
      case kHistogram:
         if (size <= 3) {
            CreateHistogram(Form("%s_%s", prefix, GetName()));
         } else {
            AliInfo(Form("[%s] Added %d > 3 axes. Creating a sparse histogram", GetName(), size));
            fOutputType = kHistogramSparse;
            CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
         }
         return kTRUE;
      case kHistogramSparse:
         CreateHistogramSparse(Form("%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.
//

   Int_t size = fAxes.GetEntries();
   AliInfo(Form("Sparse histogram name = '%s', with %d axes", name, size));

   // 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 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::FillEvent(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;
   }

   // compute & fill
   ComputeValues(event, valueList);
   FillHistogram();
   return kTRUE;
}

//________________________________________________________________________________________
Bool_t AliRsnMiniOutput::FillMother(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;

   // compute & fill
   ComputeValues(event, valueList);
   FillHistogram();
   return kTRUE;
}

//________________________________________________________________________________________
Int_t AliRsnMiniOutput::FillPair(AliRsnMiniEvent *event1, AliRsnMiniEvent *event2, TClonesArray *valueList, Bool_t refFirst)
{
//
// Loops on the passed mini-event, and for each pair of particles
// which satisfy the charge and cut requirements defined here, add an entry.
// Returns the number of successful fillings.
// Last argument tells if the reference event for event-based values is the first or the second.
//

   // check computation type
   Bool_t okComp = kFALSE;
   if (fComputation == kTrackPair)         okComp = kTRUE;
   if (fComputation == kTrackPairMix)      okComp = kTRUE;
   if (fComputation == kTrackPairRotated1) okComp = kTRUE;
   if (fComputation == kTrackPairRotated2) okComp = kTRUE;
   if (fComputation == kTruePair)          okComp = kTRUE;
   if (!okComp) {
      AliError(Form("[%s] This method can be called only for pair-based computations", GetName()));
      return kFALSE;
   }

   // loop variables
   Int_t i1, i2, start, nadded = 0;
   AliRsnMiniParticle *p1, *p2;

   // it is necessary to know if criteria for the two daughters are the same
   // and if the two events are the same or not (mixing)
   //Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fCutID[0] == fCutID[1]));
   Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fDaughter[0] == fDaughter[1]));
   Bool_t sameEvent = (event1->ID() == event2->ID());

   TString selList1  = "";
   TString selList2  = "";
   Int_t   n1 = event1->CountParticles(fSel1, fCharge[0], fCutID[0]);
   Int_t   n2 = event2->CountParticles(fSel2, fCharge[1], fCutID[1]);
   for (i1 = 0; i1 < n1; i1++) selList1.Append(Form("%d ", fSel1[i1]));
   for (i2 = 0; i2 < n2; i2++) selList2.Append(Form("%d ", fSel2[i2]));
   AliDebugClass(1, Form("[%10s] Part #1: [%s] -- evID %6d -- charge = %c -- cut ID = %d --> %4d tracks (%s)", GetName(), (event1 == event2 ? "def" : "mix"), event1->ID(), fCharge[0], fCutID[0], n1, selList1.Data()));
   AliDebugClass(1, Form("[%10s] Part #2: [%s] -- evID %6d -- charge = %c -- cut ID = %d --> %4d tracks (%s)", GetName(), (event1 == event2 ? "def" : "mix"), event2->ID(), fCharge[1], fCutID[1], n2, selList2.Data()));
   if (!n1 || !n2) {
      AliDebugClass(1, "No pairs to mix");
      return 0;
   }

   // external loop
   for (i1 = 0; i1 < n1; i1++) {
      p1 = event1->GetParticle(fSel1[i1]);
      //p1 = event1->GetParticle(i1);
      //if (p1->Charge() != fCharge[0]) continue;
      //if (!p1->HasCutBit(fCutID[0])) continue;
      // define starting point for inner loop
      // if daughter selection criteria (charge, cuts) are the same
      // and the two events coincide, internal loop must start from
      // the first track *after* current one;
      // otherwise it starts from the beginning
      start = ((sameEvent && sameCriteria) ? i1 + 1 : 0);
      AliDebugClass(2, Form("Start point = %d", start));
      // internal loop
      for (i2 = start; i2 < n2; i2++) {
         p2 = event2->GetParticle(fSel2[i2]);
         //p2 = event2->GetParticle(i2);
         //if (p2->Charge() != fCharge[1]) continue;
         //if (!p2->HasCutBit(fCutID[1])) continue;
         // avoid to mix a particle with itself
         if (sameEvent && (p1->Index() == p2->Index())) {
            AliDebugClass(2, "Skipping same index");
            continue;
         }
         // sum momenta
         fPair.Fill(p1, p2, GetMass(0), GetMass(1), fMotherMass);
         // do rotation if needed
         if (fComputation == kTrackPairRotated1) fPair.InvertP(kTRUE);
         if (fComputation == kTrackPairRotated2) fPair.InvertP(kFALSE);
         // if required, check that this is a true pair
         if (fComputation == kTruePair) {
            if (fPair.Mother() < 0)  {
               continue;
            } else if (TMath::Abs(fPair.MotherPDG()) != fMotherPDG) {
               continue;
            }
            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;
            if (!decayMatch) continue;
         }
         // check pair against cuts
         if (fPairCuts) {
            if (!fPairCuts->IsSelected(&fPair)) continue;
         }
         // get computed values & fill histogram
         nadded++;
         if (refFirst) ComputeValues(event1, valueList); else ComputeValues(event2, valueList);
         FillHistogram();
      } // end internal loop
   } // end external loop

   AliDebugClass(1, Form("Pairs added in total = %4d", nadded));
   return nadded;
}

//________________________________________________________________________________________
void 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);
   }
}

//________________________________________________________________________________________
void 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;
   }
   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");
   }
}
Commit	Line	Data
03d23846	1	//
	2	// Mini-Output
	3	// All the definitions needed for building a RSN histogram
	4	// including:
	5	// -- properties of resonance (mass, PDG code if needed)
	6	// -- properties of daughters (assigned mass, charges)
	7	// -- definition of output histogram
61f275d1	8	//
03d23846	9
	10	#include "Riostream.h"
	11
	12	#include "TH1.h"
	13	#include "TH2.h"
	14	#include "TH3.h"
	15	#include "TList.h"
45aa62b9	16	#include "TMath.h"
03d23846	17	#include "THnSparse.h"
	18	#include "TString.h"
	19	#include "TClonesArray.h"
	20
	21	#include "AliRsnMiniParticle.h"
	22	#include "AliRsnMiniPair.h"
	23	#include "AliRsnMiniEvent.h"
	24
	25	#include "AliLog.h"
	26	#include "AliRsnCutSet.h"
	27	#include "AliRsnMiniAxis.h"
	28	#include "AliRsnMiniOutput.h"
	29	#include "AliRsnMiniValue.h"
	30
	31	ClassImp(AliRsnMiniOutput)
	32
	33	//__________________________________________________________________________________________________
	34	AliRsnMiniOutput::AliRsnMiniOutput() :
	35	TNamed(),
	36	fOutputType(kTypes),
	37	fComputation(kComputations),
	38	fMotherPDG(0),
	39	fMotherMass(0.0),
	40	fPairCuts(0x0),
	41	fOutputID(-1),
	42	fAxes("AliRsnMiniAxis", 0),
	43	fComputed(0),
	44	fPair(),
a2455d2a	45	fList(0x0),
	46	fSel1(0),
	47	fSel2(0)
03d23846	48	{
	49	//
	50	// Constructor
	51	//
	52
	53	fCutID[0] = fCutID[1] = -1;
	54	fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
	55	fCharge[0] = fCharge[1] = 0;
	56	}
	57
	58	//__________________________________________________________________________________________________
	59	AliRsnMiniOutput::AliRsnMiniOutput(const char *name, EOutputType type, EComputation src) :
	60	TNamed(name, ""),
	61	fOutputType(type),
	62	fComputation(src),
	63	fMotherPDG(0),
	64	fMotherMass(0.0),
	65	fPairCuts(0x0),
	66	fOutputID(-1),
	67	fAxes("AliRsnMiniAxis", 0),
	68	fComputed(0),
	69	fPair(),
a2455d2a	70	fList(0x0),
	71	fSel1(0),
	72	fSel2(0)
03d23846	73	{
	74	//
	75	// Constructor
	76	//
	77
	78	fCutID[0] = fCutID[1] = -1;
	79	fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
	80	fCharge[0] = fCharge[1] = 0;
	81	}
	82
	83	//__________________________________________________________________________________________________
	84	AliRsnMiniOutput::AliRsnMiniOutput(const char name, const char outType, const char *compType) :
	85	TNamed(name, ""),
	86	fOutputType(kTypes),
	87	fComputation(kComputations),
	88	fMotherPDG(0),
	89	fMotherMass(0.0),
	90	fPairCuts(0x0),
	91	fOutputID(-1),
	92	fAxes("AliRsnMiniAxis", 0),
	93	fComputed(0),
	94	fPair(),
a2455d2a	95	fList(0x0),
	96	fSel1(0),
	97	fSel2(0)
03d23846	98	{
	99	//
	100	// Constructor, with a more user friendly implementation, where
	101	// the user sets the type of output and computations through conventional strings:
	102	//
	103	// Output:
6aaeb33c	104	// -- "HIST" --> common histogram (up to 3 dimensions)
6aaeb33c	105	// -- "SPARSE" --> sparse histogram
61f275d1	106	//
61f275d1	107	// Computation:
6aaeb33c	108	// -- "EVENT" --> event-only computations
	109	// -- "PAIR" --> track pair computations (default)
	110	// -- "MIX" --> event mixing (like track pair, but different events)
	111	// -- "ROTATE1" --> rotated background (rotate first track)
	112	// -- "ROTATE2" --> rotated background (rotate second track)
	113	// -- "TRUE" --> true pairs (like track pair, but checking that come from same mother)
	114	// -- "MOTHER" --> mother (loop on MC directly for mothers --> denominator of efficiency)
03d23846	115	//
	116
	117	TString input;
61f275d1	118
03d23846	119	// understand output type
	120	input = outType;
	121	input.ToUpper();
	122	if (!input.CompareTo("HIST"))
	123	fOutputType = kHistogram;
	124	else if (!input.CompareTo("SPARSE"))
	125	fOutputType = kHistogramSparse;
	126	else
	127	AliWarning(Form("String '%s' does not define a meaningful output type", outType));
61f275d1	128
03d23846	129	// understand computation type
	130	input = compType;
	131	input.ToUpper();
	132	if (!input.CompareTo("EVENT"))
	133	fComputation = kEventOnly;
	134	else if (!input.CompareTo("PAIR"))
	135	fComputation = kTrackPair;
	136	else if (!input.CompareTo("MIX"))
	137	fComputation = kTrackPairMix;
6aaeb33c	138	else if (!input.CompareTo("ROTATE1"))
	139	fComputation = kTrackPairRotated1;
	140	else if (!input.CompareTo("ROTATE2"))
	141	fComputation = kTrackPairRotated2;
03d23846	142	else if (!input.CompareTo("TRUE"))
	143	fComputation = kTruePair;
	144	else if (!input.CompareTo("MOTHER"))
	145	fComputation = kMother;
	146	else
	147	AliWarning(Form("String '%s' does not define a meaningful computation type", compType));
61f275d1	148
03d23846	149	fCutID[0] = fCutID[1] = -1;
	150	fDaughter[0] = fDaughter[1] = AliRsnDaughter::kUnknown;
	151	fCharge[0] = fCharge[1] = 0;
	152	}
	153
	154	//__________________________________________________________________________________________________
	155	AliRsnMiniOutput::AliRsnMiniOutput(const AliRsnMiniOutput &copy) :
	156	TNamed(copy),
	157	fOutputType(copy.fOutputType),
	158	fComputation(copy.fComputation),
	159	fMotherPDG(copy.fMotherPDG),
	160	fMotherMass(copy.fMotherMass),
	161	fPairCuts(copy.fPairCuts),
	162	fOutputID(copy.fOutputID),
	163	fAxes(copy.fAxes),
	164	fComputed(copy.fComputed),
	165	fPair(),
a2455d2a	166	fList(copy.fList),
	167	fSel1(0),
	168	fSel2(0)
03d23846	169	{
	170	//
	171	// Copy constructor
	172	//
	173
	174	Int_t i;
	175	for (i = 0; i < 2; i++) {
	176	fCutID[i] = copy.fCutID[i];
	177	fDaughter[i] = copy.fDaughter[i];
	178	fCharge[i] = copy.fCharge[i];
	179	}
	180	}
	181
	182	//__________________________________________________________________________________________________
61f275d1	183	AliRsnMiniOutput &AliRsnMiniOutput::operator=(const AliRsnMiniOutput &copy)
03d23846	184	{
	185	//
	186	// Assignment operator
	187	//
61f275d1	188	if (this == &copy)
	189	return *this;
	190	fOutputType = copy.fOutputType;
	191	fComputation = copy.fComputation;
	192	fMotherPDG = copy.fMotherPDG;
	193	fMotherMass = copy.fMotherMass;
	194	fPairCuts = copy.fPairCuts;
	195	fOutputID = copy.fOutputID;
	196	fAxes = copy.fAxes;
	197	fComputed = copy.fComputed;
	198	fList = copy.fList;
	199
	200	Int_t i;
	201	for (i = 0; i < 2; i++) {
	202	fCutID[i] = copy.fCutID[i];
	203	fDaughter[i] = copy.fDaughter[i];
	204	fCharge[i] = copy.fCharge[i];
	205	}
	206
	207	fSel1.Set(0);
	208	fSel2.Set(0);
	209
03d23846	210	return (*this);
	211	}
	212
	213
	214	//__________________________________________________________________________________________________
	215	void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t min, Double_t max)
	216	{
	217	//
	218	// Create a new axis reference
	219	//
	220
61f275d1	221	Int_t size = fAxes.GetEntries();
03d23846	222	new (fAxes[size]) AliRsnMiniAxis(i, nbins, min, max);
	223	}
	224
	225	//__________________________________________________________________________________________________
	226	void AliRsnMiniOutput::AddAxis(Int_t i, Double_t min, Double_t max, Double_t step)
	227	{
	228	//
	229	// Create a new axis reference
	230	//
	231
61f275d1	232	Int_t size = fAxes.GetEntries();
03d23846	233	new (fAxes[size]) AliRsnMiniAxis(i, min, max, step);
	234	}
	235
	236	//__________________________________________________________________________________________________
	237	void AliRsnMiniOutput::AddAxis(Int_t i, Int_t nbins, Double_t *values)
	238	{
	239	//
	240	// Create a new axis reference
	241	//
	242
61f275d1	243	Int_t size = fAxes.GetEntries();
03d23846	244	new (fAxes[size]) AliRsnMiniAxis(i, nbins, values);
	245	}
	246
	247	//__________________________________________________________________________________________________
	248	Bool_t AliRsnMiniOutput::Init(const char prefix, TList list)
	249	{
	250	//
	251	// Initialize properly the histogram and add it to the argument list
	252	//
	253
	254	if (!list) {
	255	AliError("Required an output list");
	256	return kFALSE;
	257	}
61f275d1	258
03d23846	259	fList = list;
45aa62b9	260	Int_t size = fAxes.GetEntries();
	261	if (size < 1) {
	262	AliWarning(Form("[%s] Cannot initialize histogram with less than 1 axis", GetName()));
	263	return kFALSE;
	264	}
03d23846	265
	266	switch (fOutputType) {
	267	case kHistogram:
d573d2fb	268	if (size <= 3) {
	269	CreateHistogram(Form("%s_%s", prefix, GetName()));
	270	} else {
	271	AliInfo(Form("[%s] Added %d > 3 axes. Creating a sparse histogram", GetName(), size));
	272	fOutputType = kHistogramSparse;
	273	CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
	274	}
03d23846	275	return kTRUE;
03d23846	276	case kHistogramSparse:
6aaeb33c	277	CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
03d23846	278	return kTRUE;
	279	default:
	280	AliError("Wrong output histogram definition");
	281	return kFALSE;
	282	}
	283	}
	284
	285	//__________________________________________________________________________________________________
	286	void AliRsnMiniOutput::CreateHistogram(const char *name)
	287	{
	288	//
	289	// Initialize the 'default' TH1 output object.
	290	// In case one of the expected axes is NULL, the initialization fails.
	291	//
	292
	293	Int_t size = fAxes.GetEntries();
	294	AliInfo(Form("Histogram name = '%s', with %d axes", name, size));
	295
	296	// we expect to have maximum 3 axes in this case
	297	AliRsnMiniAxis xAxis = 0x0, yAxis = 0x0, *zAxis = 0x0;
61f275d1	298	if (size >= 1) xAxis = (AliRsnMiniAxis *)fAxes[0];
	299	if (size >= 2) yAxis = (AliRsnMiniAxis *)fAxes[1];
	300	if (size >= 3) zAxis = (AliRsnMiniAxis *)fAxes[2];
	301
03d23846	302	// create histogram depending on the number of axes
	303	TH1 *h1 = 0x0;
	304	if (xAxis && yAxis && zAxis) {
	305	h1 = new TH3F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray(), zAxis->NBins(), zAxis->BinArray());
	306	} else if (xAxis && yAxis) {
	307	h1 = new TH2F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray());
	308	} else if (xAxis) {
	309	h1 = new TH1F(name, "", xAxis->NBins(), xAxis->BinArray());
	310	} else {
	311	AliError("No axis was initialized");
	312	return;
	313	}
61f275d1	314
03d23846	315	// switch the correct computation of errors
	316	if (h1 && fList) {
	317	h1->Sumw2();
	318	fList->Add(h1);
	319	fOutputID = fList->IndexOf(h1);
	320	}
	321	}
	322
	323	//________________________________________________________________________________________
	324	void AliRsnMiniOutput::CreateHistogramSparse(const char *name)
	325	{
	326	//
	327	// Initialize the THnSparse output object.
	328	// In case one of the expected axes is NULL, the initialization fails.
	329	//
	330
d573d2fb	331	Int_t size = fAxes.GetEntries();
d573d2fb	332	AliInfo(Form("Sparse histogram name = '%s', with %d axes", name, size));
61f275d1	333
03d23846	334	// retrieve binnings and sizes of all axes
	335	// since the check for null values is done in Init(),
	336	// we assume that here they must all be well defined
03d23846	337	Int_t i, *nbins = new Int_t[size];
03d23846	338	for (i = 0; i < size; i++) {
61f275d1	339	AliRsnMiniAxis axis = (AliRsnMiniAxis )fAxes[i];
03d23846	340	nbins[i] = axis->NBins();
	341	}
	342
	343	// create fHSparseogram
	344	THnSparseF *h1 = new THnSparseF(name, "", size, nbins);
	345
	346	// update the various axes using the definitions given in the array of axes here
	347	for (i = 0; i < size; i++) {
61f275d1	348	AliRsnMiniAxis axis = (AliRsnMiniAxis )fAxes[i];
03d23846	349	h1->GetAxis(i)->Set(nbins[i], axis->BinArray());
	350	}
	351
	352	// clear heap
	353	delete [] nbins;
61f275d1	354
03d23846	355	// add to list
	356	if (h1 && fList) {
	357	h1->Sumw2();
	358	fList->Add(h1);
	359	fOutputID = fList->IndexOf(h1);
	360	}
	361	}
	362
	363	//________________________________________________________________________________________
45aa62b9	364	Bool_t AliRsnMiniOutput::FillEvent(AliRsnMiniEvent event, TClonesArray valueList)
03d23846	365	{
	366	//
	367	// Compute values for event-based computations (does not use the pair)
	368	//
	369
	370	// check computation type
	371	if (fComputation != kEventOnly) {
	372	AliError("This method can be called only for event-based computations");
	373	return kFALSE;
	374	}
	375
45aa62b9	376	// compute & fill
	377	ComputeValues(event, valueList);
	378	FillHistogram();
	379	return kTRUE;
03d23846	380	}
	381
	382	//________________________________________________________________________________________
45aa62b9	383	Bool_t AliRsnMiniOutput::FillMother(const AliRsnMiniPair pair, AliRsnMiniEvent event, TClonesArray *valueList)
03d23846	384	{
	385	//
	386	// Compute values for mother-based computations
	387	//
	388
	389	// check computation type
	390	if (fComputation != kMother) {
	391	AliError("This method can be called only for mother-based computations");
	392	return kFALSE;
	393	}
61f275d1	394
03d23846	395	// copy passed pair info
03d23846	396	fPair = (*pair);
61f275d1	397
03d23846	398	// check pair against cuts
	399	if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE;
	400
45aa62b9	401	// compute & fill
	402	ComputeValues(event, valueList);
	403	FillHistogram();
	404	return kTRUE;
03d23846	405	}
	406
	407	//________________________________________________________________________________________
d573d2fb	408	Int_t AliRsnMiniOutput::FillPair(AliRsnMiniEvent event1, AliRsnMiniEvent event2, TClonesArray *valueList, Bool_t refFirst)
03d23846	409	{
03d23846	410	//
45aa62b9	411	// Loops on the passed mini-event, and for each pair of particles
	412	// which satisfy the charge and cut requirements defined here, add an entry.
	413	// Returns the number of successful fillings.
d573d2fb	414	// Last argument tells if the reference event for event-based values is the first or the second.
03d23846	415	//
	416
	417	// check computation type
6aaeb33c	418	Bool_t okComp = kFALSE;
45aa62b9	419	if (fComputation == kTrackPair) okComp = kTRUE;
	420	if (fComputation == kTrackPairMix) okComp = kTRUE;
	421	if (fComputation == kTrackPairRotated1) okComp = kTRUE;
	422	if (fComputation == kTrackPairRotated2) okComp = kTRUE;
	423	if (fComputation == kTruePair) okComp = kTRUE;
6aaeb33c	424	if (!okComp) {
45aa62b9	425	AliError(Form("[%s] This method can be called only for pair-based computations", GetName()));
03d23846	426	return kFALSE;
03d23846	427	}
61f275d1	428
45aa62b9	429	// loop variables
45aa62b9	430	Int_t i1, i2, start, nadded = 0;
45aa62b9	431	AliRsnMiniParticle p1, p2;
61f275d1	432
45aa62b9	433	// it is necessary to know if criteria for the two daughters are the same
45aa62b9	434	// and if the two events are the same or not (mixing)
7196ee4f	435	//Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fCutID[0] == fCutID[1]));
7196ee4f	436	Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fDaughter[0] == fDaughter[1]));
45aa62b9	437	Bool_t sameEvent = (event1->ID() == event2->ID());
61f275d1	438
9e3a9020	439	TString selList1 = "";
9e3a9020	440	TString selList2 = "";
a2455d2a	441	Int_t n1 = event1->CountParticles(fSel1, fCharge[0], fCutID[0]);
	442	Int_t n2 = event2->CountParticles(fSel2, fCharge[1], fCutID[1]);
	443	for (i1 = 0; i1 < n1; i1++) selList1.Append(Form("%d ", fSel1[i1]));
	444	for (i2 = 0; i2 < n2; i2++) selList2.Append(Form("%d ", fSel2[i2]));
	445	AliDebugClass(1, Form("[%10s] Part #1: [%s] -- evID %6d -- charge = %c -- cut ID = %d --> %4d tracks (%s)", GetName(), (event1 == event2 ? "def" : "mix"), event1->ID(), fCharge[0], fCutID[0], n1, selList1.Data()));
	446	AliDebugClass(1, Form("[%10s] Part #2: [%s] -- evID %6d -- charge = %c -- cut ID = %d --> %4d tracks (%s)", GetName(), (event1 == event2 ? "def" : "mix"), event2->ID(), fCharge[1], fCutID[1], n2, selList2.Data()));
9e3a9020	447	if (!n1 \|\| !n2) {
d573d2fb	448	AliDebugClass(1, "No pairs to mix");
	449	return 0;
	450	}
61f275d1	451
45aa62b9	452	// external loop
45aa62b9	453	for (i1 = 0; i1 < n1; i1++) {
a2455d2a	454	p1 = event1->GetParticle(fSel1[i1]);
9e3a9020	455	//p1 = event1->GetParticle(i1);
	456	//if (p1->Charge() != fCharge[0]) continue;
	457	//if (!p1->HasCutBit(fCutID[0])) continue;
45aa62b9	458	// define starting point for inner loop
	459	// if daughter selection criteria (charge, cuts) are the same
	460	// and the two events coincide, internal loop must start from
	461	// the first track after current one;
	462	// otherwise it starts from the beginning
	463	start = ((sameEvent && sameCriteria) ? i1 + 1 : 0);
7196ee4f	464	AliDebugClass(2, Form("Start point = %d", start));
45aa62b9	465	// internal loop
45aa62b9	466	for (i2 = start; i2 < n2; i2++) {
a2455d2a	467	p2 = event2->GetParticle(fSel2[i2]);
9e3a9020	468	//p2 = event2->GetParticle(i2);
	469	//if (p2->Charge() != fCharge[1]) continue;
	470	//if (!p2->HasCutBit(fCutID[1])) continue;
45aa62b9	471	// avoid to mix a particle with itself
45aa62b9	472	if (sameEvent && (p1->Index() == p2->Index())) {
d573d2fb	473	AliDebugClass(2, "Skipping same index");
45aa62b9	474	continue;
	475	}
	476	// sum momenta
	477	fPair.Fill(p1, p2, GetMass(0), GetMass(1), fMotherMass);
	478	// do rotation if needed
	479	if (fComputation == kTrackPairRotated1) fPair.InvertP(kTRUE);
	480	if (fComputation == kTrackPairRotated2) fPair.InvertP(kFALSE);
	481	// if required, check that this is a true pair
	482	if (fComputation == kTruePair) {
	483	if (fPair.Mother() < 0) {
	484	continue;
	485	} else if (TMath::Abs(fPair.MotherPDG()) != fMotherPDG) {
	486	continue;
	487	}
	488	Bool_t decayMatch = kFALSE;
	489	if (p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
	490	decayMatch = kTRUE;
	491	if (p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
	492	decayMatch = kTRUE;
	493	if (!decayMatch) continue;
	494	}
	495	// check pair against cuts
	496	if (fPairCuts) {
9e7b94f5	497	if (!fPairCuts->IsSelected(&fPair)) continue;
45aa62b9	498	}
45aa62b9	499	// get computed values & fill histogram
9e3a9020	500	nadded++;
61f275d1	501	if (refFirst) ComputeValues(event1, valueList); else ComputeValues(event2, valueList);
45aa62b9	502	FillHistogram();
45aa62b9	503	} // end internal loop
45aa62b9	504	} // end external loop
61f275d1	505
d573d2fb	506	AliDebugClass(1, Form("Pairs added in total = %4d", nadded));
45aa62b9	507	return nadded;
03d23846	508	}
	509
	510	//________________________________________________________________________________________
45aa62b9	511	void AliRsnMiniOutput::ComputeValues(AliRsnMiniEvent event, TClonesArray valueList)
03d23846	512	{
	513	//
	514	// Using the arguments and the internal 'fPair' data member,
	515	// compute all values to be stored in the histogram
	516	//
	517
	518	// check size of computed array
	519	Int_t size = fAxes.GetEntries();
	520	if (fComputed.GetSize() != size) fComputed.Set(size);
	521
	522	Int_t i, ival, nval = valueList->GetEntries();
61f275d1	523
03d23846	524	for (i = 0; i < size; i++) {
03d23846	525	fComputed[i] = 1E20;
61f275d1	526	AliRsnMiniAxis axis = (AliRsnMiniAxis )fAxes[i];
03d23846	527	if (!axis) {
	528	AliError("Null axis");
	529	continue;
	530	}
	531	ival = axis->GetValueID();
	532	if (ival < 0 \|\| ival >= nval) {
	533	AliError(Form("Required value #%d, while maximum is %d", ival, nval));
	534	continue;
	535	}
61f275d1	536	AliRsnMiniValue val = (AliRsnMiniValue )valueList->At(ival);
03d23846	537	if (!val) {
	538	AliError(Form("Value in position #%d is NULL", ival));
	539	continue;
61f275d1	540	}
03d23846	541	// if none of the above exit points is taken, compute value
	542	fComputed[i] = val->Eval(&fPair, event);
	543	}
03d23846	544	}
	545
	546	//________________________________________________________________________________________
45aa62b9	547	void AliRsnMiniOutput::FillHistogram()
03d23846	548	{
	549	//
	550	// Fills the internal histogram using the current values stored in the
	551	// 'fComputed' array, in the order as they are stored, up to the max
	552	// dimension of the initialized histogram itself.
	553	//
	554
	555	// retrieve object from list
	556	if (!fList) {
	557	AliError("List pointer is NULL");
45aa62b9	558	return;
03d23846	559	}
	560	TObject *obj = fList->At(fOutputID);
	561
	562	if (obj->InheritsFrom(TH1F::Class())) {
61f275d1	563	((TH1F *)obj)->Fill(fComputed[0]);
03d23846	564	} else if (obj->InheritsFrom(TH2F::Class())) {
61f275d1	565	((TH2F *)obj)->Fill(fComputed[0], fComputed[1]);
03d23846	566	} else if (obj->InheritsFrom(TH3F::Class())) {
61f275d1	567	((TH3F *)obj)->Fill(fComputed[0], fComputed[1], fComputed[2]);
03d23846	568	} else if (obj->InheritsFrom(THnSparseF::Class())) {
61f275d1	569	((THnSparseF *)obj)->Fill(fComputed.GetArray());
03d23846	570	} else {
03d23846	571	AliError("No output initialized");
03d23846	572	}
03d23846	573	}