[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:
         CreateHistogram(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.
//

   // 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::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)
{
//
// 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.
//

   // 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 sameEvent = (event1->ID() == event2->ID());
   
   //Int_t nsel1 = event1->CountParticles(fCharge[0], fCutID[0]);
   //Int_t nsel2 = event2->CountParticles(fCharge[1], fCutID[1]);
   //cout << "Charge #1 = " << fCharge[0] << " cut bit #1 = " << fCutID[0] << ", tracks = " << nsel1 << endl;
   //cout << "Charge #2 = " << fCharge[1] << " cut bit #2 = " << fCutID[1] << ", tracks = " << nsel2 << endl;
   //if (!nsel1 || !nsel2) {
   //   cout << "Nothing to mix" << endl;
   //   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);
      // internal loop
      for (i2 = start; i2 < n2; i2++) {
         p2 = event2->GetParticle(i2);
         if (p2->Charge() != fCharge[1]) continue;
         if (!p2->HasCutBit(fCutID[1])) continue;
         //cout << "Matching: " << i1 << ' ' << i2 << endl;
         // avoid to mix a particle with itself
         if (sameEvent && (p1->Index() == p2->Index())) {
            //cout << "-- skipping same index" << endl;
            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)) {
               //cout << "-- cuts not passed: Y = " << TMath::Abs(fPair.Y(0)) << endl;
               continue;
            }
         }
         // get computed values & fill histogram
         ComputeValues(event1, valueList);
         FillHistogram();
         nadded++;
      } // end internal loop
   } // end external loop
   
   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:
0efb7042	256	CreateHistogram(Form("%s_%s", prefix, GetName()));
03d23846	257	return kTRUE;
03d23846	258	case kHistogramSparse:
6aaeb33c	259	CreateHistogramSparse(Form("%s_%s", prefix, GetName()));
03d23846	260	return kTRUE;
	261	default:
	262	AliError("Wrong output histogram definition");
	263	return kFALSE;
	264	}
	265	}
	266
	267	//__________________________________________________________________________________________________
	268	void AliRsnMiniOutput::CreateHistogram(const char *name)
	269	{
	270	//
	271	// Initialize the 'default' TH1 output object.
	272	// In case one of the expected axes is NULL, the initialization fails.
	273	//
	274
	275	Int_t size = fAxes.GetEntries();
	276	AliInfo(Form("Histogram name = '%s', with %d axes", name, size));
	277
	278	// we expect to have maximum 3 axes in this case
	279	AliRsnMiniAxis xAxis = 0x0, yAxis = 0x0, *zAxis = 0x0;
	280	if (size >= 1) xAxis = (AliRsnMiniAxis*)fAxes[0];
	281	if (size >= 2) yAxis = (AliRsnMiniAxis*)fAxes[1];
	282	if (size >= 3) zAxis = (AliRsnMiniAxis*)fAxes[2];
	283
	284	// create histogram depending on the number of axes
	285	TH1 *h1 = 0x0;
	286	if (xAxis && yAxis && zAxis) {
	287	h1 = new TH3F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray(), zAxis->NBins(), zAxis->BinArray());
	288	} else if (xAxis && yAxis) {
	289	h1 = new TH2F(name, "", xAxis->NBins(), xAxis->BinArray(), yAxis->NBins(), yAxis->BinArray());
	290	} else if (xAxis) {
	291	h1 = new TH1F(name, "", xAxis->NBins(), xAxis->BinArray());
	292	} else {
	293	AliError("No axis was initialized");
	294	return;
	295	}
	296
	297	// switch the correct computation of errors
	298	if (h1 && fList) {
	299	h1->Sumw2();
	300	fList->Add(h1);
	301	fOutputID = fList->IndexOf(h1);
	302	}
	303	}
	304
	305	//________________________________________________________________________________________
	306	void AliRsnMiniOutput::CreateHistogramSparse(const char *name)
	307	{
	308	//
	309	// Initialize the THnSparse output object.
	310	// In case one of the expected axes is NULL, the initialization fails.
	311	//
	312
	313	// retrieve binnings and sizes of all axes
	314	// since the check for null values is done in Init(),
	315	// we assume that here they must all be well defined
	316	Int_t size = fAxes.GetEntries();
	317	Int_t i, *nbins = new Int_t[size];
	318	for (i = 0; i < size; i++) {
	319	AliRsnMiniAxis axis = (AliRsnMiniAxis)fAxes[i];
	320	nbins[i] = axis->NBins();
	321	}
	322
	323	// create fHSparseogram
324	THnSparseF *h1 = new THnSparseF(name, "", size, nbins);
325
326	// update the various axes using the definitions given in the array of axes here
327	for (i = 0; i < size; i++) {
328	AliRsnMiniAxis axis = (AliRsnMiniAxis)fAxes[i];
329	h1->GetAxis(i)->Set(nbins[i], axis->BinArray());
330	}
331
332	// clear heap
333	delete [] nbins;
334
335	// add to list
336	if (h1 && fList) {
337	h1->Sumw2();
338	fList->Add(h1);
339	fOutputID = fList->IndexOf(h1);
340	}
341	}
342
343	//________________________________________________________________________________________
45aa62b9	344	Bool_t AliRsnMiniOutput::FillEvent(AliRsnMiniEvent event, TClonesArray valueList)
03d23846	345	{
	346	//
	347	// Compute values for event-based computations (does not use the pair)
	348	//
	349
	350	// check computation type
	351	if (fComputation != kEventOnly) {
	352	AliError("This method can be called only for event-based computations");
	353	return kFALSE;
	354	}
	355
45aa62b9	356	// compute & fill
	357	ComputeValues(event, valueList);
	358	FillHistogram();
	359	return kTRUE;
03d23846	360	}
	361
	362	//________________________________________________________________________________________
45aa62b9	363	Bool_t AliRsnMiniOutput::FillMother(const AliRsnMiniPair pair, AliRsnMiniEvent event, TClonesArray *valueList)
03d23846	364	{
	365	//
	366	// Compute values for mother-based computations
	367	//
	368
	369	// check computation type
	370	if (fComputation != kMother) {
	371	AliError("This method can be called only for mother-based computations");
	372	return kFALSE;
	373	}
	374
	375	// copy passed pair info
	376	fPair = (*pair);
	377
	378	// check pair against cuts
	379	if (fPairCuts) if (!fPairCuts->IsSelected(&fPair)) return kFALSE;
	380
45aa62b9	381	// compute & fill
	382	ComputeValues(event, valueList);
	383	FillHistogram();
	384	return kTRUE;
03d23846	385	}
	386
	387	//________________________________________________________________________________________
45aa62b9	388	Int_t AliRsnMiniOutput::FillPair(AliRsnMiniEvent event1, AliRsnMiniEvent event2, TClonesArray *valueList)
03d23846	389	{
03d23846	390	//
45aa62b9	391	// Loops on the passed mini-event, and for each pair of particles
	392	// which satisfy the charge and cut requirements defined here, add an entry.
	393	// Returns the number of successful fillings.
03d23846	394	//
	395
	396	// check computation type
6aaeb33c	397	Bool_t okComp = kFALSE;
45aa62b9	398	if (fComputation == kTrackPair) okComp = kTRUE;
	399	if (fComputation == kTrackPairMix) okComp = kTRUE;
	400	if (fComputation == kTrackPairRotated1) okComp = kTRUE;
	401	if (fComputation == kTrackPairRotated2) okComp = kTRUE;
	402	if (fComputation == kTruePair) okComp = kTRUE;
6aaeb33c	403	if (!okComp) {
45aa62b9	404	AliError(Form("[%s] This method can be called only for pair-based computations", GetName()));
03d23846	405	return kFALSE;
	406	}
	407
45aa62b9	408	// loop variables
	409	Int_t i1, i2, start, nadded = 0;
	410	Int_t n1 = event1->Particles().GetEntriesFast();
	411	Int_t n2 = event2->Particles().GetEntriesFast();
	412	AliRsnMiniParticle p1, p2;
03d23846	413
45aa62b9	414	// it is necessary to know if criteria for the two daughters are the same
	415	// and if the two events are the same or not (mixing)
	416	Bool_t sameCriteria = ((fCharge[0] == fCharge[1]) && (fCutID[0] == fCutID[1]));
	417	Bool_t sameEvent = (event1->ID() == event2->ID());
6aaeb33c	418
45aa62b9	419	//Int_t nsel1 = event1->CountParticles(fCharge[0], fCutID[0]);
	420	//Int_t nsel2 = event2->CountParticles(fCharge[1], fCutID[1]);
	421	//cout << "Charge #1 = " << fCharge[0] << " cut bit #1 = " << fCutID[0] << ", tracks = " << nsel1 << endl;
	422	//cout << "Charge #2 = " << fCharge[1] << " cut bit #2 = " << fCutID[1] << ", tracks = " << nsel2 << endl;
	423	//if (!nsel1 \|\| !nsel2) {
	424	// cout << "Nothing to mix" << endl;
	425	// return 0;
	426	//}
03d23846	427
45aa62b9	428	// external loop
	429	for (i1 = 0; i1 < n1; i1++) {
	430	p1 = event1->GetParticle(i1);
	431	if (p1->Charge() != fCharge[0]) continue;
	432	if (!p1->HasCutBit(fCutID[0])) continue;
	433	// define starting point for inner loop
	434	// if daughter selection criteria (charge, cuts) are the same
	435	// and the two events coincide, internal loop must start from
	436	// the first track after current one;
	437	// otherwise it starts from the beginning
	438	start = ((sameEvent && sameCriteria) ? i1 + 1 : 0);
	439	// internal loop
	440	for (i2 = start; i2 < n2; i2++) {
	441	p2 = event2->GetParticle(i2);
	442	if (p2->Charge() != fCharge[1]) continue;
	443	if (!p2->HasCutBit(fCutID[1])) continue;
	444	//cout << "Matching: " << i1 << ' ' << i2 << endl;
	445	// avoid to mix a particle with itself
	446	if (sameEvent && (p1->Index() == p2->Index())) {
	447	//cout << "-- skipping same index" << endl;
	448	continue;
	449	}
	450	// sum momenta
	451	fPair.Fill(p1, p2, GetMass(0), GetMass(1), fMotherMass);
	452	// do rotation if needed
	453	if (fComputation == kTrackPairRotated1) fPair.InvertP(kTRUE);
	454	if (fComputation == kTrackPairRotated2) fPair.InvertP(kFALSE);
	455	// if required, check that this is a true pair
	456	if (fComputation == kTruePair) {
	457	if (fPair.Mother() < 0) {
	458	continue;
	459	} else if (TMath::Abs(fPair.MotherPDG()) != fMotherPDG) {
	460	continue;
	461	}
	462	Bool_t decayMatch = kFALSE;
	463	if (p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
	464	decayMatch = kTRUE;
	465	if (p2->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[0]) && p1->PDGAbs() == AliRsnDaughter::SpeciesPDG(fDaughter[1]))
	466	decayMatch = kTRUE;
	467	if (!decayMatch) continue;
	468	}
	469	// check pair against cuts
	470	if (fPairCuts) {
	471	if (!fPairCuts->IsSelected(&fPair)) {
	472	//cout << "-- cuts not passed: Y = " << TMath::Abs(fPair.Y(0)) << endl;
	473	continue;
	474	}
	475	}
	476	// get computed values & fill histogram
	477	ComputeValues(event1, valueList);
	478	FillHistogram();
	479	nadded++;
	480	} // end internal loop
	481	} // end external loop
03d23846	482
45aa62b9	483	return nadded;
03d23846	484	}
	485
	486	//________________________________________________________________________________________
45aa62b9	487	void AliRsnMiniOutput::ComputeValues(AliRsnMiniEvent event, TClonesArray valueList)
03d23846	488	{
	489	//
	490	// Using the arguments and the internal 'fPair' data member,
	491	// compute all values to be stored in the histogram
	492	//
	493
	494	// check size of computed array
	495	Int_t size = fAxes.GetEntries();
	496	if (fComputed.GetSize() != size) fComputed.Set(size);
	497
	498	Int_t i, ival, nval = valueList->GetEntries();
	499
	500	for (i = 0; i < size; i++) {
	501	fComputed[i] = 1E20;
	502	AliRsnMiniAxis axis = (AliRsnMiniAxis)fAxes[i];
	503	if (!axis) {
	504	AliError("Null axis");
	505	continue;
	506	}
	507	ival = axis->GetValueID();
	508	if (ival < 0 \|\| ival >= nval) {
	509	AliError(Form("Required value #%d, while maximum is %d", ival, nval));
	510	continue;
	511	}
	512	AliRsnMiniValue val = (AliRsnMiniValue)valueList->At(ival);
	513	if (!val) {
	514	AliError(Form("Value in position #%d is NULL", ival));
	515	continue;
	516	}
	517	// if none of the above exit points is taken, compute value
	518	fComputed[i] = val->Eval(&fPair, event);
	519	}
03d23846	520	}
	521
	522	//________________________________________________________________________________________
45aa62b9	523	void AliRsnMiniOutput::FillHistogram()
03d23846	524	{
	525	//
	526	// Fills the internal histogram using the current values stored in the
	527	// 'fComputed' array, in the order as they are stored, up to the max
	528	// dimension of the initialized histogram itself.
	529	//
	530
	531	// retrieve object from list
	532	if (!fList) {
	533	AliError("List pointer is NULL");
45aa62b9	534	return;
03d23846	535	}
	536	TObject *obj = fList->At(fOutputID);
	537
	538	if (obj->InheritsFrom(TH1F::Class())) {
	539	((TH1F*)obj)->Fill(fComputed[0]);
	540	} else if (obj->InheritsFrom(TH2F::Class())) {
	541	((TH2F*)obj)->Fill(fComputed[0], fComputed[1]);
	542	} else if (obj->InheritsFrom(TH3F::Class())) {
	543	((TH3F*)obj)->Fill(fComputed[0], fComputed[1], fComputed[2]);
	544	} else if (obj->InheritsFrom(THnSparseF::Class())) {
	545	((THnSparseF*)obj)->Fill(fComputed.GetArray());
	546	} else {
	547	AliError("No output initialized");
03d23846	548	}
03d23846	549	}