[u/mrichter/AliRoot.git] / PWG2 / RESONANCES / AliRsnAnalysis.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/
 
//-------------------------------------------------------------------------
//                     Class AliRsnAnalysis
//            Reconstruction and analysis of a binary resonance
// ........................................
// ........................................
// ........................................
// ........................................
// ........................................
// 
// author: A. Pulvirenti             (email: alberto.pulvirenti@ct.infn.it)
//-------------------------------------------------------------------------

#include <Riostream.h>

#include <TH1.h>
#include <TTree.h>
#include <TDatabasePDG.h>

#include "AliRsnDaughter.h"
#include "AliRsnDaughterCut.h"
#include "AliRsnDaughterCutPair.h"
#include "AliRsnEvent.h"
#include "AliRsnAnalysis.h"

ClassImp(AliRsnAnalysis)

//--------------------------------------------------------------------------------------------------------
AliRsnAnalysis::AliRsnAnalysis() :
  TObject(),
  fRejectFakes(kFALSE),
  fNBins(0),
  fHistoMin(0.0),
  fHistoMax(0.0),
  fTrueMotherPDG(0),
  fMixPairDefs(0x0),
  fMixHistograms(0x0),
  fPairDefs(0x0),
  fHistograms(0x0),
  fPairCuts(0x0),
  fEventsTree(0x0)
{
//
// Constructor
// Initializes all pointers and collections to NULL.
//
	Int_t i;
	for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i] = 0;
}
//--------------------------------------------------------------------------------------------------------
AliRsnAnalysis::AliRsnAnalysis(const AliRsnAnalysis &copy) : 
  TObject(copy),
  fRejectFakes(copy.fRejectFakes),
  fNBins(copy.fNBins),
  fHistoMin(copy.fHistoMin),
  fHistoMax(copy.fHistoMax),
  fTrueMotherPDG(copy.fTrueMotherPDG),
  fMixPairDefs(0x0),
  fMixHistograms(0x0),
  fPairDefs(0x0),
  fHistograms(0x0),
  fPairCuts(0x0),
  fEventsTree(0x0)
{
//
// Copy constructor
// Initializes all pointers and collections to NULL anyway,
// but copies some settings from the argument.
//
	Int_t i;
	for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i] = 0;
}
//--------------------------------------------------------------------------------------------------------
void AliRsnAnalysis::AddCutPair(AliRsnDaughterCutPair *cut)
{
//
// Add a cut on pairs.
// This cut is global for all pairs.
//
	if (!fPairCuts) fPairCuts = new TObjArray(0);
	fPairCuts->AddLast(cut);
}
//--------------------------------------------------------------------------------------------------------
void AliRsnAnalysis::AddCutSingle(AliPID::EParticleType type, AliRsnDaughterCut *cut)
{
//
// Add a cut on single particles.
// This cut must be specified for each particle type.
//
	if (type >= AliPID::kElectron && type <= AliPID::kProton) {
		if (!fCuts[type]) fCuts[type] = new TObjArray(0);
		fCuts[type]->AddLast(cut);
	}
}
//--------------------------------------------------------------------------------------------------------
void AliRsnAnalysis::AddMixPairDef
(AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2)
{
//
// Adds a new pair definition to create a new histogram,
// for the event mixing step.
// If the pair defs array is NULL, it is initialized here.
//
	if (!fMixPairDefs) fMixPairDefs = new TObjArray(0);
	fMixPairDefs->AddLast( new AliPairDef(p1, sign1, p2, sign2, 0, kFALSE) );
}
//--------------------------------------------------------------------------------------------------------
void AliRsnAnalysis::AddPairDef
(AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2, Bool_t onlyTrue)
{
//
// Adds a new pair definition to create a new histogram,
// for the signal evaluation (same event) step.
// If the pair defs array is NULL, it is initialized here.
//
	if (!fPairDefs) fPairDefs = new TObjArray(0);
	fPairDefs->AddLast( new AliPairDef(p1, sign1, p2, sign2, fTrueMotherPDG, onlyTrue) );
}
//--------------------------------------------------------------------------------------------------------
void AliRsnAnalysis::Clear(Option_t* /* option */)
{
//
// Clear heap
//
	fHistograms->Clear("C");
	fPairDefs->Clear("C");
		
	fMixHistograms->Clear("C");
	fMixPairDefs->Clear("C");
	
	Int_t i;
	for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i]->Clear("C");
	fPairCuts->Clear("C");
}
//--------------------------------------------------------------------------------------------------------
Stat_t AliRsnAnalysis::EventMix(Int_t nmix, Int_t multDiffMax, Double_t vzDiffMax, Bool_t compareTotMult)
{
//
// Performs event mixing.
// It takes the array of fMixPairDefs and stores results in fMixHistograms.
// First parameter defines how many events must be mixed with each one.
// Events to be mixed are chosen using the other parameters:
//
// - multDiffMax defines the maximum allowed difference in particle multiplicity,
//   a) when last argument is kFALSE (default) the multiplicity comparison is made with respect
//      to the particles of 'type 2' in the pair
//   b) when last argument is kTRUE, the multiplicity comparison is made with respect of total
//      particle multiplicity in the events
//
// - vzDiffMax defines maximum allowed difference in Z coordinate of prim. vertex.
//
// If one wants to exchange the particle types, one has to add both combinations of particles
// as different pair defs.
//
// EXAMPLE:
// analysis->AddMixPairDef(AliRsnEvent::kPion, '+', AliRsnEvent::kKaon, '-');
// analysis->AddMixPairDef(AliRsnEvent::kKaon, '-', AliRsnEvent::kPion, '+');
//
	// allocate the histograms array
	Int_t i, npairdefs = (Int_t)fMixPairDefs->GetEntries();
	fMixHistograms = new TObjArray(npairdefs);
	TObjArray &histos = *fMixHistograms;
	AliPairDef *pd = 0;
	for (i = 0; i < npairdefs; i++) {
		pd = (AliPairDef*)fMixPairDefs->At(i);
		histos[i] = new TH1D(Form("hmix_%s", pd->GetName()), pd->GetTitle(), fNBins, fHistoMin, fHistoMax);
	}
	
	// Link events branch
	AliRsnEvent *event2 = new AliRsnEvent;
	fEventsTree->SetBranchAddress("events", &event2);
	
	// define variables to store data about particles
	Stat_t nPairs = 0;
	Int_t iev, ievmix, nEvents = (Int_t)fEventsTree->GetEntries();
	
	// loop on events
	Int_t nmixed, mult1, mult2, diffMult;
	Double_t vz1, vz2, diffVz;
	for (iev = 0; iev < nEvents; iev++) {
	
		// get event
		event2->Clear("DELETE");
		fEventsTree->GetEntry(iev);
		
		// copy this event into 'event 1'
		AliRsnEvent *event1 = new AliRsnEvent(*event2);
		
		// get Z position of primary vertex
		vz1 = event1->GetPrimaryVertexZ();
		
		// if total multiplicities must be used
		// it is computed here
		if (compareTotMult) {
			mult1 = event1->GetMultiplicity();
		}
		else {
			mult1 = 0;
		}
		
		// message
		if (iev % 10 == 0) cout << "\rMixing event " << iev << flush;
		
		// loop on pair definitions
		for (i = 0; i < npairdefs; i++) {
			
			// get pair definition
			pd = (AliPairDef*)fMixPairDefs->At(i);
			
			// get histogram reference
			TH1D *histogram = (TH1D*)histos[i];
			
			// get multiplicity of particle type '2' in event '1'
			if (!mult1) {
				mult1 = (Int_t)event1->GetTracks(pd->GetSign2(), pd->GetParticle2())->GetEntries();
			}
			
			// loop on events for mixing
			nmixed = 0;
			ievmix = iev;
			while (nmixed < nmix) {
				
				// get other event (it starts from iev + 1, and 
				// loops again to 0 if reachs end of tree)
				ievmix++;
				if (ievmix >= nEvents) ievmix -= nEvents;
				
				// if event-2-mix is equal to 'iev', then
				// a complete loop has been done, and no more
				// mixing can be done, even if they are too few
				// then the loop breaks anyway
				if (ievmix == iev) break;
				
				// get other event
				event2->Clear("DELETE");
				fEventsTree->GetEntry(ievmix);
				
				// get event stats related to event 2
				vz2 = event2->GetPrimaryVertexZ();
				if (compareTotMult) {
					mult2 = event2->GetMultiplicity();
				}
				else {
					mult2 = event2->GetTracks(pd->GetSign2(), pd->GetParticle2())->GetEntries();
				}
				
				// fill histogram if checks are passed
				diffMult = TMath::Abs(mult1 - mult2);
				diffVz = TMath::Abs(vz1 - vz2);
				if (diffVz <= vzDiffMax && diffMult <= multDiffMax) {
					//cout << ievmix << " " << flush;
					nPairs += Compute(pd, histogram, event1, event2);
					nmixed++;
				}
			}
		}
		
		event1->Clear("DELETE");
		delete event1;
		event1 = 0;
		
	} // end loop on events
	cout << endl;
	
	return nPairs;
}
//--------------------------------------------------------------------------------------------------------
Stat_t AliRsnAnalysis::Process()
{
//
// Reads the list 'fPairDefs', and builds an inv-mass histogram for each definition.
// For each event, particle of 'type 1' are combined with particles of 'type 2' as 
// defined in each pair definition specified in the list, taking the list of both types
// from the same event.
// This method is supposed to evaluate the 'signal' of the resonance, containing the peak.
// It can also be used when one wants to evaluate the 'background' with the 'wrong sign'
// particle pairs.
//
	// allocate the histograms array in order to contain 
	// as many objects as the number of pair definitionss
	Int_t i, npairdefs = (Int_t)fPairDefs->GetEntries();
	fHistograms = new TObjArray(npairdefs);
	TObjArray &histos = *fHistograms;
	
	// loop over pair defs in order to retrieve the particle species chosen
	// which are used to set an unique name for the output histogram
	// There will be a direct correspondence between the inder of pairdef in its array
	// and the corresponding histogram in the 'fHistograms' list.
	AliPairDef *pd = 0;
	for (i = 0; i < npairdefs; i++) {
		pd = (AliPairDef*)fPairDefs->At(i);
		histos[i] = new TH1D(Form("h_%s", pd->GetName()), pd->GetTitle(), fNBins, fHistoMin, fHistoMax);
	}
	
	// Link events branch
	AliRsnEvent *event = new AliRsnEvent;
	fEventsTree->SetBranchAddress("events", &event);
	
	// define variables to store data about particles
	Stat_t nPairs = 0;
	Int_t iev, nEvents = (Int_t)fEventsTree->GetEntries();
	
	// loop on events
	for (iev = 0; iev < nEvents; iev++) {
	
		// get event
		event->Clear("DELETE");
		fEventsTree->GetEntry(iev);
		if (iev % 1000 == 0) cout << "\rProcessing event " << iev << flush;
			
		// loop over the collection of pair defs
		for (i = 0; i < npairdefs; i++) {
			pd = (AliPairDef*)fPairDefs->At(i);
			TH1D *histogram = (TH1D*)histos[i];
			// invoke AliPairDef method to fill histogram
			nPairs += Compute(pd, histogram, event, event);
		}
		
	} // end loop on events
	cout << endl;
	
	return nPairs;
}
//--------------------------------------------------------------------------------------------------------
void AliRsnAnalysis::WriteHistograms() const
{
//
// Writes histograms in current directory
//
	TH1D *histogram;
	TObjArrayIter iter(fHistograms);
	while ( (histogram = (TH1D*)iter.Next()) ) {
		histogram->Write();
	}
	
	if (fMixHistograms) {
		TObjArrayIter iterMix(fMixHistograms);
		while ( (histogram = (TH1D*)iterMix.Next()) ) {
			histogram->Write();
		}
	}
}
//--------------------------------------------------------------------------------------------------------
Stat_t AliRsnAnalysis::Compute
(AliPairDef *pd, TH1D* &histogram, AliRsnEvent *event1, AliRsnEvent *event2)
{
//
// Adds to the specified histogram the invariant mass spectrum calculated taking
// particles of type 1 from event 1, and particles of type 2 from event 2.
// Events can be equal (signal) or different (background with event mixing).
//
	// define two 'cursor' objects
	AliRsnDaughter *track1 = 0, *track2 = 0;
	
	// define iterators for the two collections
	if (!event1 || !event1->GetTracks(pd->GetSign1(), pd->GetParticle1())) {
		Error("Compute", "Null pointer for particle 1 collection");
		return 0.0;
	}
	if (!event2 || !event2->GetTracks(pd->GetSign2(), pd->GetParticle2())) {
		Error("Compute", "Null pointer for particle 2 collection");
		return 0.0;
	}
	TObjArrayIter iter1(event1->GetTracks(pd->GetSign1(), pd->GetParticle1()));
	TObjArrayIter iter2(event2->GetTracks(pd->GetSign2(), pd->GetParticle2()));
	
	// define temporary variables for better code readability
	Stat_t nPairs = 0;
	
	// loop on particle of type 1 (in event 1)
	while ( (track1 = (AliRsnDaughter*)iter1.Next()) ) {

		if (fRejectFakes && (track1->GetLabel() < 0)) continue;
		if (!SingleCutCheck(pd->GetParticle1(), track1)) continue;
		
		iter2.Reset();
		
		// loop on particles of type 2 (in event 2)
		while ( (track2 = (AliRsnDaughter*)iter2.Next()) ) {
			
			if (fRejectFakes && (track2->GetLabel() < 0)) continue;
			if (event1 == event2 && track1->GetIndex() == track2->GetIndex()) continue;
			if (!SingleCutCheck(pd->GetParticle2(), track2)) continue;
			if (!PairCutCheck(track1, track2)) continue;
					
			/*
			// check
			Char_t sign1 = (track1->GetSign() > 0) ? '+' : '-';
			Char_t sign2 = (track2->GetSign() > 0) ? '+' : '-';
			Info("Compute", "Particle 1: PDG = %d, sign = %c --- Particle 2: PDG = %d, sign = %c", track1->GetTruePDG(), sign1, track2->GetTruePDG(), sign2);
			*/
				
			// total 4-momentum
			track1->SetMass(pd->GetMass1());
			track2->SetMass(pd->GetMass2());
			AliRsnDaughter sum = AliRsnDaughter::Sum(*track1, *track2);
				
			// if the choice to get ONLY true pairs is selected, a check is made on the mothers
			if (pd->GetOnlyTrue()) {
				// the 'sum' AliRsnDaughter object is initialized with
				// the PDG code of the common mother (if it is the case)
				if (TMath::Abs(sum.GetMotherPDG()) == TMath::Abs(fTrueMotherPDG)) {
					histogram->Fill(sum.GetMass());
					nPairs++;
				}
			}
			else {
				histogram->Fill(sum.GetMass());
				nPairs++;
			}	
				
		} // end loop 2
				
	} // end loop 1
	
	return nPairs;
}
//--------------------------------------------------------------------------------------------------------
Bool_t AliRsnAnalysis::SingleCutCheck(Int_t itype, AliRsnDaughter *track) const
{
//
// Checks a track against single particle cuts (if defined)
//
	if (!fCuts[itype]) return kTRUE;
	
	TObjArrayIter iter(fCuts[itype]);
	AliRsnDaughterCut *cut = 0;
	while ( (cut = (AliRsnDaughterCut*)iter.Next()) ) {
		if (!cut->Pass(track)) return kFALSE;
	}
	
	return kTRUE;
}
//--------------------------------------------------------------------------------------------------------
Bool_t AliRsnAnalysis::PairCutCheck(AliRsnDaughter *track1, AliRsnDaughter *track2) const
{
//
// Checks a pair against pair cuts (if defined)
//
	if (!fPairCuts) return kTRUE;
	
	TObjArrayIter iter(fPairCuts);
	AliRsnDaughterCutPair *cut = 0;
	while ( (cut = (AliRsnDaughterCutPair*)iter.Next()) ) {
		if (!cut->Pass(track1, track2)) return kFALSE;
	}
	
	return kTRUE;
}
//--------------------------------------------------------------------------------------------------------
AliRsnAnalysis::AliPairDef::AliPairDef
(AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2, Int_t pdgMother, Bool_t onlyTrue) :
  TNamed(),
  fOnlyTrue(onlyTrue),
  fTrueMotherPDG(pdgMother),
  fMass1(0.0),
  fSign1(sign1),
  fParticle1(p1),
  fMass2(0.0),
  fSign2(sign2),
  fParticle2(p2)
{
//
// Constructor for nested class
//
	if (!fOnlyTrue) fTrueMotherPDG = 0;
		
	// instance a PDG database to recovery true masses of particles
	TDatabasePDG *db = TDatabasePDG::Instance();
	Int_t pdg1 = AliPID::ParticleCode((Int_t)p1);
	Int_t pdg2 = AliPID::ParticleCode((Int_t)p2);
	fMass1 = db->GetParticle(pdg1)->Mass();
	fMass2 = db->GetParticle(pdg2)->Mass();
	
	// set name according to the chosen particles
	fName.Append(Form("%s(%c)_%s(%c)", ParticleName(p1), sign1, ParticleName(p2), sign2));
	fTitle.Append(Form("Inv. mass for %s(%c) and %s(%c)", ParticleName(p1), sign1, ParticleName(p2), sign2));
	if (onlyTrue) {
		fName.Append("_true");
		fTitle.Append(" (true pairs)");
	}
}
//--------------------------------------------------------------------------------------------------------
Text_t* AliRsnAnalysis::AliPairDef::ParticleName(AliPID::EParticleType part) const
{
//
// [PRIVATE]
// Returns the name of the particle in text format
//
	if (part == AliPID::kElectron) return ("E");
	else if (part == AliPID::kMuon) return ("Mu");
	else if (part == AliPID::kPion) return ("Pi");
	else if (part == AliPID::kKaon) return ("K");
	else if (part == AliPID::kProton) return ("P");
	else {
		Warning("ParticleName", "Unrecognized value of EParticle argument");
		return ("???");
	}
}
//--------------------------------------------------------------------------------------------------------
Commit	Line	Data
2f1637fb	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	//-------------------------------------------------------------------------
	17	// Class AliRsnAnalysis
	18	// Reconstruction and analysis of a binary resonance
2f769150	19	// ........................................
	20	// ........................................
	21	// ........................................
	22	// ........................................
	23	// ........................................
2f1637fb	24	//
	25	// author: A. Pulvirenti (email: alberto.pulvirenti@ct.infn.it)
	26	//-------------------------------------------------------------------------
	27
	28	#include <Riostream.h>
	29
	30	#include <TH1.h>
	31	#include <TTree.h>
2f1637fb	32	#include <TDatabasePDG.h>
	33
	34	#include "AliRsnDaughter.h"
	35	#include "AliRsnDaughterCut.h"
c37c6481	36	#include "AliRsnDaughterCutPair.h"
2f1637fb	37	#include "AliRsnEvent.h"
	38	#include "AliRsnAnalysis.h"
	39
	40	ClassImp(AliRsnAnalysis)
	41
	42	//--------------------------------------------------------------------------------------------------------
c37c6481	43	AliRsnAnalysis::AliRsnAnalysis() :
	44	TObject(),
	45	fRejectFakes(kFALSE),
	46	fNBins(0),
	47	fHistoMin(0.0),
	48	fHistoMax(0.0),
	49	fTrueMotherPDG(0),
	50	fMixPairDefs(0x0),
	51	fMixHistograms(0x0),
	52	fPairDefs(0x0),
	53	fHistograms(0x0),
	54	fPairCuts(0x0),
	55	fEventsTree(0x0)
2f769150	56	{
2f1637fb	57	//
	58	// Constructor
	59	// Initializes all pointers and collections to NULL.
	60	//
2f1637fb	61	Int_t i;
	62	for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i] = 0;
	63	}
	64	//--------------------------------------------------------------------------------------------------------
c37c6481	65	AliRsnAnalysis::AliRsnAnalysis(const AliRsnAnalysis &copy) :
	66	TObject(copy),
	67	fRejectFakes(copy.fRejectFakes),
	68	fNBins(copy.fNBins),
	69	fHistoMin(copy.fHistoMin),
	70	fHistoMax(copy.fHistoMax),
	71	fTrueMotherPDG(copy.fTrueMotherPDG),
	72	fMixPairDefs(0x0),
	73	fMixHistograms(0x0),
	74	fPairDefs(0x0),
	75	fHistograms(0x0),
	76	fPairCuts(0x0),
	77	fEventsTree(0x0)
	78	{
	79	//
	80	// Copy constructor
	81	// Initializes all pointers and collections to NULL anyway,
	82	// but copies some settings from the argument.
	83	//
	84	Int_t i;
	85	for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i] = 0;
	86	}
	87	//--------------------------------------------------------------------------------------------------------
	88	void AliRsnAnalysis::AddCutPair(AliRsnDaughterCutPair *cut)
2f769150	89	{
2f1637fb	90	//
	91	// Add a cut on pairs.
	92	// This cut is global for all pairs.
	93	//
2f1637fb	94	if (!fPairCuts) fPairCuts = new TObjArray(0);
2f1637fb	95	fPairCuts->AddLast(cut);
	96	}
	97	//--------------------------------------------------------------------------------------------------------
	98	void AliRsnAnalysis::AddCutSingle(AliPID::EParticleType type, AliRsnDaughterCut *cut)
2f769150	99	{
2f1637fb	100	//
	101	// Add a cut on single particles.
	102	// This cut must be specified for each particle type.
	103	//
2f1637fb	104	if (type >= AliPID::kElectron && type <= AliPID::kProton) {
	105	if (!fCuts[type]) fCuts[type] = new TObjArray(0);
	106	fCuts[type]->AddLast(cut);
	107	}
	108	}
	109	//--------------------------------------------------------------------------------------------------------
	110	void AliRsnAnalysis::AddMixPairDef
	111	(AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2)
2f769150	112	{
2f1637fb	113	//
	114	// Adds a new pair definition to create a new histogram,
	115	// for the event mixing step.
	116	// If the pair defs array is NULL, it is initialized here.
	117	//
2f1637fb	118	if (!fMixPairDefs) fMixPairDefs = new TObjArray(0);
	119	fMixPairDefs->AddLast( new AliPairDef(p1, sign1, p2, sign2, 0, kFALSE) );
	120	}
	121	//--------------------------------------------------------------------------------------------------------
	122	void AliRsnAnalysis::AddPairDef
	123	(AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2, Bool_t onlyTrue)
2f769150	124	{
2f1637fb	125	//
	126	// Adds a new pair definition to create a new histogram,
	127	// for the signal evaluation (same event) step.
	128	// If the pair defs array is NULL, it is initialized here.
	129	//
2f1637fb	130	if (!fPairDefs) fPairDefs = new TObjArray(0);
	131	fPairDefs->AddLast( new AliPairDef(p1, sign1, p2, sign2, fTrueMotherPDG, onlyTrue) );
	132	}
	133	//--------------------------------------------------------------------------------------------------------
	134	void AliRsnAnalysis::Clear(Option_t* /* option */)
2f769150	135	{
2f1637fb	136	//
	137	// Clear heap
	138	//
2f1637fb	139	fHistograms->Clear("C");
	140	fPairDefs->Clear("C");
	141
	142	fMixHistograms->Clear("C");
	143	fMixPairDefs->Clear("C");
	144
	145	Int_t i;
	146	for (i = 0; i < AliPID::kSPECIES; i++) fCuts[i]->Clear("C");
	147	fPairCuts->Clear("C");
	148	}
	149	//--------------------------------------------------------------------------------------------------------
	150	Stat_t AliRsnAnalysis::EventMix(Int_t nmix, Int_t multDiffMax, Double_t vzDiffMax, Bool_t compareTotMult)
2f769150	151	{
2f1637fb	152	//
	153	// Performs event mixing.
	154	// It takes the array of fMixPairDefs and stores results in fMixHistograms.
	155	// First parameter defines how many events must be mixed with each one.
	156	// Events to be mixed are chosen using the other parameters:
	157	//
	158	// - multDiffMax defines the maximum allowed difference in particle multiplicity,
	159	// a) when last argument is kFALSE (default) the multiplicity comparison is made with respect
	160	// to the particles of 'type 2' in the pair
	161	// b) when last argument is kTRUE, the multiplicity comparison is made with respect of total
	162	// particle multiplicity in the events
	163	//
	164	// - vzDiffMax defines maximum allowed difference in Z coordinate of prim. vertex.
	165	//
	166	// If one wants to exchange the particle types, one has to add both combinations of particles
	167	// as different pair defs.
	168	//
	169	// EXAMPLE:
	170	// analysis->AddMixPairDef(AliRsnEvent::kPion, '+', AliRsnEvent::kKaon, '-');
	171	// analysis->AddMixPairDef(AliRsnEvent::kKaon, '-', AliRsnEvent::kPion, '+');
	172	//
2f1637fb	173	// allocate the histograms array
	174	Int_t i, npairdefs = (Int_t)fMixPairDefs->GetEntries();
	175	fMixHistograms = new TObjArray(npairdefs);
	176	TObjArray &histos = *fMixHistograms;
	177	AliPairDef *pd = 0;
	178	for (i = 0; i < npairdefs; i++) {
	179	pd = (AliPairDef*)fMixPairDefs->At(i);
	180	histos[i] = new TH1D(Form("hmix_%s", pd->GetName()), pd->GetTitle(), fNBins, fHistoMin, fHistoMax);
	181	}
	182
	183	// Link events branch
	184	AliRsnEvent *event2 = new AliRsnEvent;
	185	fEventsTree->SetBranchAddress("events", &event2);
	186
	187	// define variables to store data about particles
	188	Stat_t nPairs = 0;
	189	Int_t iev, ievmix, nEvents = (Int_t)fEventsTree->GetEntries();
	190
	191	// loop on events
	192	Int_t nmixed, mult1, mult2, diffMult;
	193	Double_t vz1, vz2, diffVz;
	194	for (iev = 0; iev < nEvents; iev++) {
	195
	196	// get event
	197	event2->Clear("DELETE");
	198	fEventsTree->GetEntry(iev);
	199
	200	// copy this event into 'event 1'
	201	AliRsnEvent event1 = new AliRsnEvent(event2);
	202
	203	// get Z position of primary vertex
	204	vz1 = event1->GetPrimaryVertexZ();
	205
	206	// if total multiplicities must be used
	207	// it is computed here
	208	if (compareTotMult) {
a8a7b929	209	mult1 = event1->GetMultiplicity();
2f1637fb	210	}
	211	else {
	212	mult1 = 0;
	213	}
	214
	215	// message
	216	if (iev % 10 == 0) cout << "\rMixing event " << iev << flush;
	217
	218	// loop on pair definitions
	219	for (i = 0; i < npairdefs; i++) {
	220
	221	// get pair definition
	222	pd = (AliPairDef*)fMixPairDefs->At(i);
	223
	224	// get histogram reference
	225	TH1D histogram = (TH1D)histos[i];
	226
	227	// get multiplicity of particle type '2' in event '1'
	228	if (!mult1) {
	229	mult1 = (Int_t)event1->GetTracks(pd->GetSign2(), pd->GetParticle2())->GetEntries();
	230	}
	231
	232	// loop on events for mixing
	233	nmixed = 0;
	234	ievmix = iev;
	235	while (nmixed < nmix) {
	236
	237	// get other event (it starts from iev + 1, and
	238	// loops again to 0 if reachs end of tree)
	239	ievmix++;
	240	if (ievmix >= nEvents) ievmix -= nEvents;
	241
	242	// if event-2-mix is equal to 'iev', then
	243	// a complete loop has been done, and no more
	244	// mixing can be done, even if they are too few
	245	// then the loop breaks anyway
	246	if (ievmix == iev) break;
	247
	248	// get other event
	249	event2->Clear("DELETE");
	250	fEventsTree->GetEntry(ievmix);
	251
	252	// get event stats related to event 2
	253	vz2 = event2->GetPrimaryVertexZ();
	254	if (compareTotMult) {
a8a7b929	255	mult2 = event2->GetMultiplicity();
2f1637fb	256	}
	257	else {
	258	mult2 = event2->GetTracks(pd->GetSign2(), pd->GetParticle2())->GetEntries();
	259	}
	260
	261	// fill histogram if checks are passed
	262	diffMult = TMath::Abs(mult1 - mult2);
	263	diffVz = TMath::Abs(vz1 - vz2);
	264	if (diffVz <= vzDiffMax && diffMult <= multDiffMax) {
	265	//cout << ievmix << " " << flush;
	266	nPairs += Compute(pd, histogram, event1, event2);
	267	nmixed++;
	268	}
	269	}
	270	}
	271
	272	event1->Clear("DELETE");
	273	delete event1;
	274	event1 = 0;
	275
	276	} // end loop on events
	277	cout << endl;
	278
	279	return nPairs;
	280	}
	281	//--------------------------------------------------------------------------------------------------------
	282	Stat_t AliRsnAnalysis::Process()
2f769150	283	{
2f1637fb	284	//
	285	// Reads the list 'fPairDefs', and builds an inv-mass histogram for each definition.
	286	// For each event, particle of 'type 1' are combined with particles of 'type 2' as
	287	// defined in each pair definition specified in the list, taking the list of both types
	288	// from the same event.
	289	// This method is supposed to evaluate the 'signal' of the resonance, containing the peak.
	290	// It can also be used when one wants to evaluate the 'background' with the 'wrong sign'
	291	// particle pairs.
	292	//
2f1637fb	293	// allocate the histograms array in order to contain
	294	// as many objects as the number of pair definitionss
	295	Int_t i, npairdefs = (Int_t)fPairDefs->GetEntries();
	296	fHistograms = new TObjArray(npairdefs);
	297	TObjArray &histos = *fHistograms;
	298
	299	// loop over pair defs in order to retrieve the particle species chosen
	300	// which are used to set an unique name for the output histogram
	301	// There will be a direct correspondence between the inder of pairdef in its array
	302	// and the corresponding histogram in the 'fHistograms' list.
	303	AliPairDef *pd = 0;
	304	for (i = 0; i < npairdefs; i++) {
	305	pd = (AliPairDef*)fPairDefs->At(i);
	306	histos[i] = new TH1D(Form("h_%s", pd->GetName()), pd->GetTitle(), fNBins, fHistoMin, fHistoMax);
	307	}
	308
	309	// Link events branch
	310	AliRsnEvent *event = new AliRsnEvent;
	311	fEventsTree->SetBranchAddress("events", &event);
	312
	313	// define variables to store data about particles
	314	Stat_t nPairs = 0;
	315	Int_t iev, nEvents = (Int_t)fEventsTree->GetEntries();
	316
	317	// loop on events
	318	for (iev = 0; iev < nEvents; iev++) {
	319
	320	// get event
	321	event->Clear("DELETE");
	322	fEventsTree->GetEntry(iev);
	323	if (iev % 1000 == 0) cout << "\rProcessing event " << iev << flush;
	324
	325	// loop over the collection of pair defs
	326	for (i = 0; i < npairdefs; i++) {
	327	pd = (AliPairDef*)fPairDefs->At(i);
	328	TH1D histogram = (TH1D)histos[i];
	329	// invoke AliPairDef method to fill histogram
	330	nPairs += Compute(pd, histogram, event, event);
	331	}
	332
	333	} // end loop on events
	334	cout << endl;
	335
	336	return nPairs;
	337	}
	338	//--------------------------------------------------------------------------------------------------------
2f769150	339	void AliRsnAnalysis::WriteHistograms() const
2f769150	340	{
2f1637fb	341	//
	342	// Writes histograms in current directory
	343	//
2f1637fb	344	TH1D *histogram;
	345	TObjArrayIter iter(fHistograms);
	346	while ( (histogram = (TH1D*)iter.Next()) ) {
	347	histogram->Write();
	348	}
	349
	350	if (fMixHistograms) {
	351	TObjArrayIter iterMix(fMixHistograms);
	352	while ( (histogram = (TH1D*)iterMix.Next()) ) {
	353	histogram->Write();
	354	}
	355	}
	356	}
	357	//--------------------------------------------------------------------------------------------------------
	358	Stat_t AliRsnAnalysis::Compute
	359	(AliPairDef pd, TH1D &histogram, AliRsnEvent event1, AliRsnEvent event2)
2f769150	360	{
2f1637fb	361	//
	362	// Adds to the specified histogram the invariant mass spectrum calculated taking
	363	// particles of type 1 from event 1, and particles of type 2 from event 2.
	364	// Events can be equal (signal) or different (background with event mixing).
	365	//
2f1637fb	366	// define two 'cursor' objects
	367	AliRsnDaughter track1 = 0, track2 = 0;
	368
	369	// define iterators for the two collections
	370	if (!event1 \|\| !event1->GetTracks(pd->GetSign1(), pd->GetParticle1())) {
	371	Error("Compute", "Null pointer for particle 1 collection");
	372	return 0.0;
	373	}
	374	if (!event2 \|\| !event2->GetTracks(pd->GetSign2(), pd->GetParticle2())) {
	375	Error("Compute", "Null pointer for particle 2 collection");
	376	return 0.0;
	377	}
	378	TObjArrayIter iter1(event1->GetTracks(pd->GetSign1(), pd->GetParticle1()));
	379	TObjArrayIter iter2(event2->GetTracks(pd->GetSign2(), pd->GetParticle2()));
	380
	381	// define temporary variables for better code readability
	382	Stat_t nPairs = 0;
	383
	384	// loop on particle of type 1 (in event 1)
	385	while ( (track1 = (AliRsnDaughter*)iter1.Next()) ) {
	386
	387	if (fRejectFakes && (track1->GetLabel() < 0)) continue;
	388	if (!SingleCutCheck(pd->GetParticle1(), track1)) continue;
	389
	390	iter2.Reset();
	391
	392	// loop on particles of type 2 (in event 2)
	393	while ( (track2 = (AliRsnDaughter*)iter2.Next()) ) {
	394
	395	if (fRejectFakes && (track2->GetLabel() < 0)) continue;
	396	if (event1 == event2 && track1->GetIndex() == track2->GetIndex()) continue;
	397	if (!SingleCutCheck(pd->GetParticle2(), track2)) continue;
	398	if (!PairCutCheck(track1, track2)) continue;
	399
	400	/*
	401	// check
	402	Char_t sign1 = (track1->GetSign() > 0) ? '+' : '-';
	403	Char_t sign2 = (track2->GetSign() > 0) ? '+' : '-';
	404	Info("Compute", "Particle 1: PDG = %d, sign = %c --- Particle 2: PDG = %d, sign = %c", track1->GetTruePDG(), sign1, track2->GetTruePDG(), sign2);
	405	*/
	406
	407	// total 4-momentum
	408	track1->SetMass(pd->GetMass1());
	409	track2->SetMass(pd->GetMass2());
2f769150	410	AliRsnDaughter sum = AliRsnDaughter::Sum(track1, track2);
2f1637fb	411
	412	// if the choice to get ONLY true pairs is selected, a check is made on the mothers
	413	if (pd->GetOnlyTrue()) {
	414	// the 'sum' AliRsnDaughter object is initialized with
	415	// the PDG code of the common mother (if it is the case)
	416	if (TMath::Abs(sum.GetMotherPDG()) == TMath::Abs(fTrueMotherPDG)) {
	417	histogram->Fill(sum.GetMass());
	418	nPairs++;
	419	}
	420	}
	421	else {
	422	histogram->Fill(sum.GetMass());
	423	nPairs++;
	424	}
	425
	426	} // end loop 2
	427
	428	} // end loop 1
	429
	430	return nPairs;
	431	}
	432	//--------------------------------------------------------------------------------------------------------
2f769150	433	Bool_t AliRsnAnalysis::SingleCutCheck(Int_t itype, AliRsnDaughter *track) const
2f769150	434	{
2f1637fb	435	//
	436	// Checks a track against single particle cuts (if defined)
	437	//
2f1637fb	438	if (!fCuts[itype]) return kTRUE;
	439
	440	TObjArrayIter iter(fCuts[itype]);
	441	AliRsnDaughterCut *cut = 0;
	442	while ( (cut = (AliRsnDaughterCut*)iter.Next()) ) {
	443	if (!cut->Pass(track)) return kFALSE;
	444	}
	445
	446	return kTRUE;
	447	}
	448	//--------------------------------------------------------------------------------------------------------
2f769150	449	Bool_t AliRsnAnalysis::PairCutCheck(AliRsnDaughter track1, AliRsnDaughter track2) const
2f769150	450	{
2f1637fb	451	//
	452	// Checks a pair against pair cuts (if defined)
	453	//
2f1637fb	454	if (!fPairCuts) return kTRUE;
	455
	456	TObjArrayIter iter(fPairCuts);
c37c6481	457	AliRsnDaughterCutPair *cut = 0;
c37c6481	458	while ( (cut = (AliRsnDaughterCutPair*)iter.Next()) ) {
2f1637fb	459	if (!cut->Pass(track1, track2)) return kFALSE;
	460	}
	461
	462	return kTRUE;
	463	}
	464	//--------------------------------------------------------------------------------------------------------
	465	AliRsnAnalysis::AliPairDef::AliPairDef
c37c6481	466	(AliPID::EParticleType p1, Char_t sign1, AliPID::EParticleType p2, Char_t sign2, Int_t pdgMother, Bool_t onlyTrue) :
	467	TNamed(),
	468	fOnlyTrue(onlyTrue),
	469	fTrueMotherPDG(pdgMother),
	470	fMass1(0.0),
	471	fSign1(sign1),
	472	fParticle1(p1),
	473	fMass2(0.0),
	474	fSign2(sign2),
	475	fParticle2(p2)
2f769150	476	{
2f1637fb	477	//
	478	// Constructor for nested class
	479	//
c37c6481	480	if (!fOnlyTrue) fTrueMotherPDG = 0;
c37c6481	481
2f1637fb	482	// instance a PDG database to recovery true masses of particles
	483	TDatabasePDG *db = TDatabasePDG::Instance();
	484	Int_t pdg1 = AliPID::ParticleCode((Int_t)p1);
	485	Int_t pdg2 = AliPID::ParticleCode((Int_t)p2);
	486	fMass1 = db->GetParticle(pdg1)->Mass();
	487	fMass2 = db->GetParticle(pdg2)->Mass();
	488
	489	// set name according to the chosen particles
	490	fName.Append(Form("%s(%c)_%s(%c)", ParticleName(p1), sign1, ParticleName(p2), sign2));
	491	fTitle.Append(Form("Inv. mass for %s(%c) and %s(%c)", ParticleName(p1), sign1, ParticleName(p2), sign2));
	492	if (onlyTrue) {
	493	fName.Append("_true");
	494	fTitle.Append(" (true pairs)");
	495	}
	496	}
	497	//--------------------------------------------------------------------------------------------------------
2f769150	498	Text_t* AliRsnAnalysis::AliPairDef::ParticleName(AliPID::EParticleType part) const
2f769150	499	{
2f1637fb	500	//
	501	// [PRIVATE]
	502	// Returns the name of the particle in text format
	503	//
2f1637fb	504	if (part == AliPID::kElectron) return ("E");
	505	else if (part == AliPID::kMuon) return ("Mu");
	506	else if (part == AliPID::kPion) return ("Pi");
	507	else if (part == AliPID::kKaon) return ("K");
	508	else if (part == AliPID::kProton) return ("P");
	509	else {
	510	Warning("ParticleName", "Unrecognized value of EParticle argument");
	511	return ("???");
	512	}
	513	}
	514	//--------------------------------------------------------------------------------------------------------