[u/mrichter/AliRoot.git] / PWG2 / 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)
{
//
// 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)
//    -- "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)
{
//
// 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
//

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