[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 <TRefArray.h>
#include <TParticle.h>
#include <TObjectTable.h>
#include <TDatabasePDG.h>

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

ClassImp(AliRsnAnalysis)

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