[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis2 / AliSPDMCTrackDensity.cxx

#include "AliSPDMCTrackDensity.h"
#include <AliMCEvent.h>
#include <AliTrackReference.h>
#include <AliStack.h>
#include <TMath.h>
#include <AliLog.h>
#include <TH3D.h>
#include <TH2D.h>
#include <TH1D.h>
#include <TList.h>
#include <TROOT.h>
#include <iostream>

//____________________________________________________________________
AliSPDMCTrackDensity::AliSPDMCTrackDensity()
  : TNamed(), 
    fUseOnlyPrimary(false), 
    fMinR(3.5), 
    fMaxR(4.5),
    fMinZ(-14.1), 
    fMaxZ(+14.1),
    fRZ(0), 
    fXYZ(0),
    fNRefs(0),
    fBinFlow(0)
{
  // Default constructor 
}

//____________________________________________________________________
AliSPDMCTrackDensity::AliSPDMCTrackDensity(const char*)
  : TNamed("spdMCTrackDensity","spdMCTrackDensity"), 
    fUseOnlyPrimary(false), 
    fMinR(3.5), 
    fMaxR(4.5),
    fMinZ(-14.1), 
    fMaxZ(+14.1),
    fRZ(0), 
    fXYZ(0),
    fNRefs(0),
    fBinFlow(0)
{
  // Normal constructor constructor 
}

//____________________________________________________________________
AliSPDMCTrackDensity::AliSPDMCTrackDensity(const AliSPDMCTrackDensity& o)
  : TNamed(o),
    fUseOnlyPrimary(o.fUseOnlyPrimary), 
    fMinR(o.fMinR), 
    fMaxR(o.fMaxR),
    fMinZ(o.fMinZ), 
    fMaxZ(o.fMaxZ),
    fRZ(o.fRZ), 
    fXYZ(o.fXYZ),
    fNRefs(o.fNRefs),
    fBinFlow(o.fBinFlow)
{
  // Normal constructor constructor 
}

//____________________________________________________________________
AliSPDMCTrackDensity&
AliSPDMCTrackDensity::operator=(const AliSPDMCTrackDensity& o)
{
  // Assignment operator 
  TNamed::operator=(o);
  fUseOnlyPrimary       = o.fUseOnlyPrimary;
  fMinR             = o.fMinR;
  fMaxR             = o.fMaxR;
  fMinZ             = o.fMinZ;
  fMaxZ             = o.fMaxZ;
  fRZ               = o.fRZ;
  fXYZ              = o.fXYZ;
  fNRefs            = o.fNRefs;
  fBinFlow          = o.fBinFlow;
  return *this;
}

//____________________________________________________________________
void
AliSPDMCTrackDensity::DefineOutput(TList* l)
{
  fRZ = new TH2D("rz", "(r,z) of used track references", 
		 50, 0, 5, 30, -15, 15);
  fRZ->SetXTitle("z [cm]");
  fRZ->SetYTitle("r [cm]");
  fRZ->SetDirectory(0);
  l->Add(fRZ);

  fXYZ = new TH3D("xyz", "(x,y,z) of used track references", 
		  100, -5., 5., 100, -5., 5., 30, -15., 15.);
  fXYZ->SetXTitle("x [cm]");
  fXYZ->SetYTitle("y [cm]");
  fXYZ->SetZTitle("z [cm]");
  fXYZ->SetDirectory(0);
  l->Add(fXYZ);

  fNRefs = new TH1D("nrefs", "Number of references used per track", 
		    11, -.5, 10.5); 
  fNRefs->SetXTitle("# of references per track");
  fNRefs->SetDirectory(0);
  l->Add(fNRefs);

  fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow", 
		      120, -3, 3, 40, -180, 180);
  fBinFlow->SetXTitle("#Delta#eta");
  fBinFlow->SetYTitle("#Delta#varphi");
  fBinFlow->SetDirectory(0);
  l->Add(fBinFlow);
}

//____________________________________________________________________
void
AliSPDMCTrackDensity::StoreParticle(AliMCParticle* particle, 
				    const AliMCParticle* mother, 
				    Int_t          refNo,
				    Double_t       vz, 
				    TH2D&     output) const
{
  // Store a particle. 
  if (refNo < 0) return;

  AliTrackReference* ref = particle->GetTrackReference(refNo);
  if (!ref) return;
    
  Double_t r = ref->R();
  Double_t x = ref->X();
  Double_t y = ref->Y();
  Double_t z = ref->Z();

  Double_t zr = z-vz;
  Double_t th = TMath::ATan2(r,zr);
  if (th < 0) th += 2*TMath::Pi();
  Double_t et = -TMath::Log(TMath::Tan(th/2));
  Double_t ph = TMath::ATan2(y,x);
  if (ph < 0) ph += 2*TMath::Pi();
  output.Fill(et,ph);

  const AliMCParticle* mp = (mother ? mother : particle);
  Double_t dEta = mp->Eta() - et;
  Double_t dPhi = (mp->Phi() - ph) * 180 / TMath::Pi();
  if (dPhi >  180) dPhi -= 360;
  if (dPhi < -180) dPhi += 360;
  fBinFlow->Fill(dEta, dPhi);
}


//____________________________________________________________________
const AliMCParticle*
AliSPDMCTrackDensity::GetMother(Int_t     iTr,
				const AliMCEvent& event) const
{
  // 
  // Track down primary mother 
  // 
  Int_t i  = iTr;
  do { 
    const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
    if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
    
    i = p->GetMother();
  } while (i > 0);

  return 0;
}  

//____________________________________________________________________
Bool_t
AliSPDMCTrackDensity::Calculate(const AliMCEvent& event, 
				Double_t          vz,
				TH2D&             output, 
				TH2D*             primary)
{
  // 
  // Filter the input kinematics and track references, using 
  // some of the ESD information
  // 
  // Parameters:
  //    input   Input ESD event
  //    event   Input MC event
  //    vz      Vertex position 
  //    output  Output ESD-like object
  //    primary Per-event histogram of primaries 
  //
  // Return:
  //    True on succes, false otherwise 
  //

  AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
  Int_t nTracks   = stack->GetNtrack();
  Int_t nPrim     = stack->GetNtrack();
  for (Int_t iTr = 0; iTr < nTracks; iTr++) { 
    AliMCParticle* particle = 
      static_cast<AliMCParticle*>(event.GetTrack(iTr));
    
    // Check the returned particle 
    if (!particle) continue;
    
    // Check if this charged and a primary 
    Bool_t isCharged = particle->Charge() != 0;
    if (!isCharged) continue;

    Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);

    // Fill 'dn/deta' histogram 
    if (isPrimary && iTr < nPrim) {
      if (primary) {
	Double_t eta = particle->Eta();
	Double_t phi = particle->Phi();
	primary->Fill(eta, phi);
      }
    }

    // Bail out if we're only processing primaries - perhaps we should
    // track back to the original primary?
    if (fUseOnlyPrimary && !isPrimary) continue;

    Int_t nTrRef  = particle->GetNumberOfTrackReferences();
    Int_t nRef    = 0;
    for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) { 
      AliTrackReference* ref = particle->GetTrackReference(iTrRef);
      
      // Check existence 
      if (!ref) continue;

      // Check that we hit an FMD element 
      if (ref->DetectorId() != AliTrackReference::kITS) 
	continue;
      
      // Get radius and z where the track reference was made 
      Double_t r = ref->R();
      Double_t x = ref->X();
      Double_t y = ref->Y();
      Double_t z = ref->Z();
      if (r > fMaxR || r < fMinR) continue;
      if (z > fMaxZ || z < fMinZ) continue;

      fRZ->Fill(z,r);
      fXYZ->Fill(x, y, z);
      nRef++;
      // Only fill first reference 
      if (nRef == 1) { 
	const AliMCParticle* mother = GetMother(iTrRef, event);
	StoreParticle(particle, mother, iTrRef, vz, output);
      }
    }
    fNRefs->Fill(nRef);
  }
  return kTRUE;
}
//____________________________________________________________________
void
AliSPDMCTrackDensity::Print(Option_t* /*option*/) const 
{
  char ind[gROOT->GetDirLevel()+1];
  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
  ind[gROOT->GetDirLevel()] = '\0';
  std::cout << ind << ClassName() << ": " << GetName() << '\n'
	    << std::boolalpha 
	    << ind << " Only primary tracks:    " << fUseOnlyPrimary << '\n'
	    << ind << " R range:                [" << fMinR << ',' << fMaxR
	    << "]\n"
	    << ind << " Z range:                [" << fMinZ << ',' << fMaxZ
	    << "]" << std::noboolalpha << std::endl;
  
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
f68d9069	1	#include "AliSPDMCTrackDensity.h"
	2	#include <AliMCEvent.h>
	3	#include <AliTrackReference.h>
	4	#include <AliStack.h>
	5	#include <TMath.h>
	6	#include <AliLog.h>
	7	#include <TH3D.h>
	8	#include <TH2D.h>
	9	#include <TH1D.h>
	10	#include <TList.h>
	11	#include <TROOT.h>
	12	#include <iostream>
	13
	14	//____________________________________________________________________
	15	AliSPDMCTrackDensity::AliSPDMCTrackDensity()
	16	: TNamed(),
	17	fUseOnlyPrimary(false),
	18	fMinR(3.5),
	19	fMaxR(4.5),
	20	fMinZ(-14.1),
	21	fMaxZ(+14.1),
	22	fRZ(0),
	23	fXYZ(0),
	24	fNRefs(0),
	25	fBinFlow(0)
	26	{
	27	// Default constructor
	28	}
	29
	30	//____________________________________________________________________
	31	AliSPDMCTrackDensity::AliSPDMCTrackDensity(const char*)
	32	: TNamed("spdMCTrackDensity","spdMCTrackDensity"),
	33	fUseOnlyPrimary(false),
	34	fMinR(3.5),
	35	fMaxR(4.5),
	36	fMinZ(-14.1),
	37	fMaxZ(+14.1),
	38	fRZ(0),
	39	fXYZ(0),
	40	fNRefs(0),
	41	fBinFlow(0)
	42	{
	43	// Normal constructor constructor
	44	}
	45
	46	//____________________________________________________________________
	47	AliSPDMCTrackDensity::AliSPDMCTrackDensity(const AliSPDMCTrackDensity& o)
	48	: TNamed(o),
	49	fUseOnlyPrimary(o.fUseOnlyPrimary),
	50	fMinR(o.fMinR),
	51	fMaxR(o.fMaxR),
	52	fMinZ(o.fMinZ),
	53	fMaxZ(o.fMaxZ),
	54	fRZ(o.fRZ),
	55	fXYZ(o.fXYZ),
	56	fNRefs(o.fNRefs),
	57	fBinFlow(o.fBinFlow)
	58	{
	59	// Normal constructor constructor
	60	}
	61
	62	//____________________________________________________________________
	63	AliSPDMCTrackDensity&
	64	AliSPDMCTrackDensity::operator=(const AliSPDMCTrackDensity& o)
65	{
66	// Assignment operator
67	TNamed::operator=(o);
68	fUseOnlyPrimary = o.fUseOnlyPrimary;
69	fMinR = o.fMinR;
70	fMaxR = o.fMaxR;
71	fMinZ = o.fMinZ;
72	fMaxZ = o.fMaxZ;
73	fRZ = o.fRZ;
74	fXYZ = o.fXYZ;
75	fNRefs = o.fNRefs;
76	fBinFlow = o.fBinFlow;
77	return *this;
78	}
79
80	//____________________________________________________________________
81	void
82	AliSPDMCTrackDensity::DefineOutput(TList* l)
83	{
84	fRZ = new TH2D("rz", "(r,z) of used track references",
85	50, 0, 5, 30, -15, 15);
86	fRZ->SetXTitle("z [cm]");
87	fRZ->SetYTitle("r [cm]");
88	fRZ->SetDirectory(0);
89	l->Add(fRZ);
90
91	fXYZ = new TH3D("xyz", "(x,y,z) of used track references",
92	100, -5., 5., 100, -5., 5., 30, -15., 15.);
93	fXYZ->SetXTitle("x [cm]");
94	fXYZ->SetYTitle("y [cm]");
95	fXYZ->SetZTitle("z [cm]");
96	fXYZ->SetDirectory(0);
97	l->Add(fXYZ);
98
99	fNRefs = new TH1D("nrefs", "Number of references used per track",
100	11, -.5, 10.5);
101	fNRefs->SetXTitle("# of references per track");
102	fNRefs->SetDirectory(0);
103	l->Add(fNRefs);
104
105	fBinFlow = new TH2D("binFlow", "#eta and #varphi bin flow",
106	120, -3, 3, 40, -180, 180);
107	fBinFlow->SetXTitle("#Delta#eta");
108	fBinFlow->SetYTitle("#Delta#varphi");
109	fBinFlow->SetDirectory(0);
110	l->Add(fBinFlow);
111	}
112
113	//____________________________________________________________________
114	void
115	AliSPDMCTrackDensity::StoreParticle(AliMCParticle* particle,
116	const AliMCParticle* mother,
117	Int_t refNo,
118	Double_t vz,
119	TH2D& output) const
120	{
121	// Store a particle.
122	if (refNo < 0) return;
123
124	AliTrackReference* ref = particle->GetTrackReference(refNo);
125	if (!ref) return;
126
127	Double_t r = ref->R();
128	Double_t x = ref->X();
129	Double_t y = ref->Y();
130	Double_t z = ref->Z();
131
132	Double_t zr = z-vz;
133	Double_t th = TMath::ATan2(r,zr);
134	if (th < 0) th += 2*TMath::Pi();
135	Double_t et = -TMath::Log(TMath::Tan(th/2));
136	Double_t ph = TMath::ATan2(y,x);
137	if (ph < 0) ph += 2*TMath::Pi();
138	output.Fill(et,ph);
139
140	const AliMCParticle* mp = (mother ? mother : particle);
141	Double_t dEta = mp->Eta() - et;
142	Double_t dPhi = (mp->Phi() - ph) * 180 / TMath::Pi();
143	if (dPhi > 180) dPhi -= 360;
144	if (dPhi < -180) dPhi += 360;
145	fBinFlow->Fill(dEta, dPhi);
146	}
147
148
149	//____________________________________________________________________
150	const AliMCParticle*
151	AliSPDMCTrackDensity::GetMother(Int_t iTr,
152	const AliMCEvent& event) const
153	{
154	//
155	// Track down primary mother
156	//
157	Int_t i = iTr;
158	do {
159	const AliMCParticle* p = static_cast<AliMCParticle*>(event.GetTrack(i));
160	if (const_cast<AliMCEvent&>(event).Stack()->IsPhysicalPrimary(i)) return p;
161
162	i = p->GetMother();
163	} while (i > 0);
164
165	return 0;
166	}
167
168	//____________________________________________________________________
169	Bool_t
170	AliSPDMCTrackDensity::Calculate(const AliMCEvent& event,
171	Double_t vz,
172	TH2D& output,
173	TH2D* primary)
174	{
175	//
176	// Filter the input kinematics and track references, using
177	// some of the ESD information
178	//
179	// Parameters:
180	// input Input ESD event
181	// event Input MC event
182	// vz Vertex position
183	// output Output ESD-like object
184	// primary Per-event histogram of primaries
185	//
186	// Return:
187	// True on succes, false otherwise
188	//
189
190	AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
191	Int_t nTracks = stack->GetNtrack();
192	Int_t nPrim = stack->GetNtrack();
193	for (Int_t iTr = 0; iTr < nTracks; iTr++) {
194	AliMCParticle* particle =
195	static_cast<AliMCParticle*>(event.GetTrack(iTr));
196
197	// Check the returned particle
198	if (!particle) continue;
199
200	// Check if this charged and a primary
201	Bool_t isCharged = particle->Charge() != 0;
202	if (!isCharged) continue;
203
204	Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
205
206	// Fill 'dn/deta' histogram
207	if (isPrimary && iTr < nPrim) {
208	if (primary) {
209	Double_t eta = particle->Eta();
210	Double_t phi = particle->Phi();
211	primary->Fill(eta, phi);
212	}
213	}
214
215	// Bail out if we're only processing primaries - perhaps we should
216	// track back to the original primary?
217	if (fUseOnlyPrimary && !isPrimary) continue;
218
219	Int_t nTrRef = particle->GetNumberOfTrackReferences();
220	Int_t nRef = 0;
221	for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
222	AliTrackReference* ref = particle->GetTrackReference(iTrRef);
223
224	// Check existence
225	if (!ref) continue;
226
227	// Check that we hit an FMD element
228	if (ref->DetectorId() != AliTrackReference::kITS)
229	continue;
230
231	// Get radius and z where the track reference was made
232	Double_t r = ref->R();
233	Double_t x = ref->X();
234	Double_t y = ref->Y();
235	Double_t z = ref->Z();
236	if (r > fMaxR \|\| r < fMinR) continue;
237	if (z > fMaxZ \|\| z < fMinZ) continue;
238
239	fRZ->Fill(z,r);
240	fXYZ->Fill(x, y, z);
241	nRef++;
242	// Only fill first reference
243	if (nRef == 1) {
244	const AliMCParticle* mother = GetMother(iTrRef, event);
245	StoreParticle(particle, mother, iTrRef, vz, output);
246	}
247	}
248	fNRefs->Fill(nRef);
249	}
250	return kTRUE;
251	}
252	//____________________________________________________________________
253	void
254	AliSPDMCTrackDensity::Print(Option_t* /option/) const
255	{
256	char ind[gROOT->GetDirLevel()+1];
257	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
258	ind[gROOT->GetDirLevel()] = '\0';
259	std::cout << ind << ClassName() << ": " << GetName() << '\n'
260	<< std::boolalpha
261	<< ind << " Only primary tracks: " << fUseOnlyPrimary << '\n'
262	<< ind << " R range: [" << fMinR << ',' << fMaxR
263	<< "]\n"
264	<< ind << " Z range: [" << fMinZ << ',' << fMaxZ
265	<< "]" << std::noboolalpha << std::endl;
266
267	}
268
269	//____________________________________________________________________
270	//
271	// EOF
272	//