[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 "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 &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;

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