[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliFMDMCTrackDensity.cxx

#include "AliFMDMCTrackDensity.h"
#include "AliESDFMD.h"
#include "AliTrackReference.h"
#include <TMath.h>
#include "AliFMDStripIndex.h"
#include "AliLog.h"
#include <TH2D.h>
#include <TH1D.h>
#include <TList.h>
#include <TROOT.h>
#include <iostream>

//____________________________________________________________________
void
AliFMDMCTrackDensity::State::Clear(Bool_t alsoCount)
{
  angle       = 0;
  oldDetector = 0;
  oldRing     = '\0';
  oldSector   = 1024;
  oldStrip    = 1024;
  startStrip  = 1024;
  nRefs       = 0;
  nStrips     = 0;
  longest     = 0x0;
  if (alsoCount) count = 0;
}

//____________________________________________________________________
AliFMDMCTrackDensity::State&
AliFMDMCTrackDensity::State::operator=(const State& o)
{
  if (&o == this) return *this;
  angle          = o.angle;
  oldDetector    = o.oldDetector;
  oldRing        = o.oldRing;
  oldSector      = o.oldSector;
  oldStrip       = o.oldStrip;
  startStrip     = o.startStrip;
  nRefs          = o.nRefs;
  nStrips        = o.nStrips;
  count          = o.count;
  longest        = o.longest;
  return *this;
}

//____________________________________________________________________
AliFMDMCTrackDensity::AliFMDMCTrackDensity()
  : AliBaseMCTrackDensity(), 
    fState(),
    fMaxConsequtiveStrips(3), 
    fNr(0), 
    fNt(0), 
    fNc(0),
    fNcr(0),
    fOutput(0)
{
  // Default constructor 
}

//____________________________________________________________________
AliFMDMCTrackDensity::AliFMDMCTrackDensity(const char*)
  : AliBaseMCTrackDensity("fmdMCTrackDensity"), 
    fState(),
    fMaxConsequtiveStrips(3), 
    fNr(0), 
    fNt(0), 
    fNc(0),
    fNcr(0),
    fOutput(0)
{
  // Normal constructor constructor 
}

//____________________________________________________________________
AliFMDMCTrackDensity::AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o)
  : AliBaseMCTrackDensity(o),
    fState(o.fState),
    fMaxConsequtiveStrips(o.fMaxConsequtiveStrips), 
    fNr(o.fNr), 
    fNt(o.fNt), 
    fNc(o.fNc),
    fNcr(o.fNcr),
    fOutput(o.fOutput)
{
  // Normal constructor constructor 
}

//____________________________________________________________________
AliFMDMCTrackDensity&
AliFMDMCTrackDensity::operator=(const AliFMDMCTrackDensity& o)
{
  // Assignment operator 
  if (&o == this) return *this; 
  AliBaseMCTrackDensity::operator=(o);
  fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
  fNr                   = o.fNr;
  fNt                   = o.fNt;
  fNc                   = o.fNc;
  fNcr                  = o.fNcr;
  fState                = o.fState;
  fOutput               = o.fOutput;

  return *this;
}

//____________________________________________________________________
void
AliFMDMCTrackDensity::CreateOutputObjects(TList* l)
{
  AliBaseMCTrackDensity::CreateOutputObjects(l);
  TList* ll = static_cast<TList*>(l->FindObject(GetTitle()));
  if (!ll) ll = l;

  fNr = new TH1D("clusterRefs", "# track references per cluster",
		 21, -.5, 20.5);
  fNr->SetXTitle("# of track references/cluster");
  fNr->SetFillColor(kRed+1);
  fNr->SetFillStyle(3001);
  fNr->SetDirectory(0);
  ll->Add(fNr);

  fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
  fNt->SetXTitle("Cluster size [strips]");
  fNt->SetFillColor(kBlue+1);
  fNt->SetFillStyle(3001);
  fNt->SetDirectory(0);
  ll->Add(fNt);

  fNc = new TH1D("nClusters", "# clusters per track", 21, -.5, 20.5);
  fNc->SetXTitle("# clusters");
  fNc->SetFillColor(kGreen+1);
  fNc->SetFillStyle(3001);
  fNc->SetDirectory(0);
  ll->Add(fNc);

  fNcr = new TH2D("clusterVsRefs", "# clusters vs. # refs", 
		  21, -.5, 20.5, 21, -.5, 20.5);
  fNcr->SetXTitle("# References");
  fNcr->SetYTitle("# Clusters");
  fNcr->SetOption("COLZ");
  fNcr->SetDirectory(0);
  ll->Add(fNcr);
  
		  
}
//____________________________________________________________________
Int_t
AliFMDMCTrackDensity::GetDetectorId() const
{
  return AliTrackReference::kFMD;
}

//____________________________________________________________________
void
AliFMDMCTrackDensity::BeginTrackRefs()
{
  fState.Clear(true);
}

//____________________________________________________________________
void
AliFMDMCTrackDensity::EndTrackRefs(Int_t nRefs)
{
  fNc->Fill(fState.count);
  fNcr->Fill(nRefs, fState.count);
  fState.Clear(true);
}
  
//____________________________________________________________________
AliTrackReference*
AliFMDMCTrackDensity::ProcessRef(AliMCParticle*       particle,
				 const AliMCParticle* mother,
				 AliTrackReference*   ref)
{
  // Get the detector coordinates 
  UShort_t d, s, t;
  Char_t r;
  AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
    
  // Calculate distance of previous reference to base of cluster 
  UShort_t nT = TMath::Abs(t - fState.startStrip) + 1;

  // Now check if we should flush to output 
  Bool_t used = false;

  // If this is a new detector/ring, then reset the other one 
  // Check if we have a valid old detectorm ring, and sector 
  if (fState.oldDetector >  0 && 
      fState.oldRing     != '\0' && 
      fState.oldSector   != 1024) {
    // New detector, new ring, or new sector 
    if (d != fState.oldDetector   || 
	r != fState.oldRing       || 
	s != fState.oldSector) {
      if (fDebug) Info("Process", "New because new sector");
      used = true;
    }
    else if (nT > fMaxConsequtiveStrips) {
      if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)", 
		       fState.nStrips, t, fState.oldStrip, fState.startStrip);
      used = true;
    }
  }
  if (used) {
    if (fDebug) 
      Info("Process", "I=%p L=%p D=%d (was %d), R=%c (was %c), "
	   "S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
	   ref, fState.longest, 
	   d, fState.oldDetector, 
	   r, fState.oldRing, 
	   s, fState.oldSector, 
	   t, fState.oldStrip, 
	   fState.nStrips, fMaxConsequtiveStrips);
    // Int_t nnT   = TMath::Abs(fState.oldStrip - fState.startStrip) + 1;
    StoreParticle(particle, mother, fState.longest);
    fState.Clear(false);
  }

  // If base of cluster not set, set it here. 
  if (fState.startStrip == 1024) fState.startStrip = t;
  
  // Calculate distance of previous reference to base of cluster 
  fState.nStrips = TMath::Abs(t - fState.startStrip) + 1;

  // Count number of track refs in this sector 
  fState.nRefs++;

  fState.oldDetector = d;
  fState.oldRing     = r;
  fState.oldSector   = s;
  fState.oldStrip    = t;

  // Debug output 
  if (fDebug) {
    if (t == fState.startStrip) 
      Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]", 
	   d, r, s, t);
    else 
      Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
	   "length=%3d (now in %3d, previous %3d)", 
	   d, r, s, fState.startStrip, fState.nStrips, t, fState.oldStrip);
  }
    
  // The longest passage is determined through the angle 
  Double_t ang  = GetTrackRefTheta(ref);
  if (ang > fState.angle) {
    fState.longest = ref;
    fState.angle   = ang;
  }

  return fState.longest;
}

//____________________________________________________________________
Double_t
AliFMDMCTrackDensity::StoreParticle(AliMCParticle*       particle, 
				    const AliMCParticle* mother, 
				    AliTrackReference*   ref) const
{
  Double_t w = 
    AliBaseMCTrackDensity::StoreParticle(particle, mother, ref);
  if (w <= 0) return w;

  // Get the detector coordinates 
  UShort_t d, s, t;
  Char_t r;
  AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);

  // Check if we have value already 
  Double_t old = fOutput->Multiplicity(d,r,s,t);

  // If invalid, force it valid 
  if (old == AliESDFMD::kInvalidMult) old = 0;

  // Increment count 
  fOutput->SetMultiplicity(d,r,s,t,old+w);

  // Fill histograms 
  fNr->Fill(fState.nRefs);
  fNt->Fill(fState.nStrips);

  fState.count++;

  return w;
}  

//____________________________________________________________________
Bool_t
AliFMDMCTrackDensity::Calculate(const AliESDFMD&  input, 
				const AliMCEvent& event,
				Double_t          vz,
				AliESDFMD&        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 
  //
  output.Clear();
  fOutput = &output;

  // Copy eta values to output 
  for (UShort_t ed = 1; ed <= 3; ed++) { 
    UShort_t nq = (ed == 1 ? 1 : 2);
    for (UShort_t eq = 0; eq < nq; eq++) {
      Char_t   er = (eq == 0 ? 'I' : 'O');
      UShort_t ns = (eq == 0 ?  20 :  40);
      UShort_t nt = (eq == 0 ? 512 : 256);
      for (UShort_t es = 0; es < ns; es++) 
	for (UShort_t et = 0; et < nt; et++) 
	  output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
    }
  }

  return ProcessTracks(event, vz, primary);
}
//____________________________________________________________________
void
AliFMDMCTrackDensity::Print(Option_t* option) const 
{
  AliBaseMCTrackDensity::Print(option);
  char ind[gROOT->GetDirLevel()+1];
  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
  ind[gROOT->GetDirLevel()] = '\0';
  std::cout << ind << " Max cluster size:       " << fMaxConsequtiveStrips
	    << std::endl;
  
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
f68d9069	1	#include "AliFMDMCTrackDensity.h"
f53fb4f6	2	#include "AliESDFMD.h"
f53fb4f6	3	#include "AliTrackReference.h"
f68d9069	4	#include <TMath.h>
f68d9069	5	#include "AliFMDStripIndex.h"
f53fb4f6	6	#include "AliLog.h"
f68d9069	7	#include <TH2D.h>
	8	#include <TH1D.h>
	9	#include <TList.h>
	10	#include <TROOT.h>
	11	#include <iostream>
4bcdcbc1	12
	13	//____________________________________________________________________
	14	void
	15	AliFMDMCTrackDensity::State::Clear(Bool_t alsoCount)
	16	{
	17	angle = 0;
	18	oldDetector = 0;
	19	oldRing = '\0';
	20	oldSector = 1024;
	21	oldStrip = 1024;
	22	startStrip = 1024;
	23	nRefs = 0;
	24	nStrips = 0;
	25	longest = 0x0;
	26	if (alsoCount) count = 0;
	27	}
	28
	29	//____________________________________________________________________
	30	AliFMDMCTrackDensity::State&
	31	AliFMDMCTrackDensity::State::operator=(const State& o)
	32	{
6ab100ec	33	if (&o == this) return *this;
4bcdcbc1	34	angle = o.angle;
	35	oldDetector = o.oldDetector;
	36	oldRing = o.oldRing;
	37	oldSector = o.oldSector;
	38	oldStrip = o.oldStrip;
	39	startStrip = o.startStrip;
	40	nRefs = o.nRefs;
	41	nStrips = o.nStrips;
	42	count = o.count;
	43	longest = o.longest;
	44	return *this;
	45	}
f68d9069	46
	47	//____________________________________________________________________
	48	AliFMDMCTrackDensity::AliFMDMCTrackDensity()
4bcdcbc1	49	: AliBaseMCTrackDensity(),
4bcdcbc1	50	fState(),
f68d9069	51	fMaxConsequtiveStrips(3),
f68d9069	52	fNr(0),
4bcdcbc1	53	fNt(0),
	54	fNc(0),
	55	fNcr(0),
	56	fOutput(0)
f68d9069	57	{
	58	// Default constructor
	59	}
	60
	61	//____________________________________________________________________
	62	AliFMDMCTrackDensity::AliFMDMCTrackDensity(const char*)
4bcdcbc1	63	: AliBaseMCTrackDensity("fmdMCTrackDensity"),
4bcdcbc1	64	fState(),
f68d9069	65	fMaxConsequtiveStrips(3),
f68d9069	66	fNr(0),
4bcdcbc1	67	fNt(0),
	68	fNc(0),
	69	fNcr(0),
	70	fOutput(0)
f68d9069	71	{
	72	// Normal constructor constructor
	73	}
	74
	75	//____________________________________________________________________
	76	AliFMDMCTrackDensity::AliFMDMCTrackDensity(const AliFMDMCTrackDensity& o)
4bcdcbc1	77	: AliBaseMCTrackDensity(o),
4bcdcbc1	78	fState(o.fState),
f68d9069	79	fMaxConsequtiveStrips(o.fMaxConsequtiveStrips),
f68d9069	80	fNr(o.fNr),
4bcdcbc1	81	fNt(o.fNt),
	82	fNc(o.fNc),
	83	fNcr(o.fNcr),
	84	fOutput(o.fOutput)
f68d9069	85	{
	86	// Normal constructor constructor
	87	}
	88
	89	//____________________________________________________________________
	90	AliFMDMCTrackDensity&
	91	AliFMDMCTrackDensity::operator=(const AliFMDMCTrackDensity& o)
	92	{
	93	// Assignment operator
d015ecfe	94	if (&o == this) return *this;
4bcdcbc1	95	AliBaseMCTrackDensity::operator=(o);
f68d9069	96	fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
	97	fNr = o.fNr;
	98	fNt = o.fNt;
4bcdcbc1	99	fNc = o.fNc;
	100	fNcr = o.fNcr;
	101	fState = o.fState;
	102	fOutput = o.fOutput;
	103
f68d9069	104	return *this;
	105	}
	106
	107	//____________________________________________________________________
	108	void
5934a3e3	109	AliFMDMCTrackDensity::CreateOutputObjects(TList* l)
f68d9069	110	{
5934a3e3	111	AliBaseMCTrackDensity::CreateOutputObjects(l);
4bcdcbc1	112	TList* ll = static_cast<TList*>(l->FindObject(GetTitle()));
	113	if (!ll) ll = l;
	114
f68d9069	115	fNr = new TH1D("clusterRefs", "# track references per cluster",
	116	21, -.5, 20.5);
	117	fNr->SetXTitle("# of track references/cluster");
4bcdcbc1	118	fNr->SetFillColor(kRed+1);
4bcdcbc1	119	fNr->SetFillStyle(3001);
f68d9069	120	fNr->SetDirectory(0);
4bcdcbc1	121	ll->Add(fNr);
f68d9069	122
	123	fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
	124	fNt->SetXTitle("Cluster size [strips]");
4bcdcbc1	125	fNt->SetFillColor(kBlue+1);
4bcdcbc1	126	fNt->SetFillStyle(3001);
f68d9069	127	fNt->SetDirectory(0);
4bcdcbc1	128	ll->Add(fNt);
	129
	130	fNc = new TH1D("nClusters", "# clusters per track", 21, -.5, 20.5);
	131	fNc->SetXTitle("# clusters");
	132	fNc->SetFillColor(kGreen+1);
	133	fNc->SetFillStyle(3001);
	134	fNc->SetDirectory(0);
	135	ll->Add(fNc);
	136
	137	fNcr = new TH2D("clusterVsRefs", "# clusters vs. # refs",
	138	21, -.5, 20.5, 21, -.5, 20.5);
	139	fNcr->SetXTitle("# References");
	140	fNcr->SetYTitle("# Clusters");
	141	fNcr->SetOption("COLZ");
	142	fNcr->SetDirectory(0);
	143	ll->Add(fNcr);
	144
	145
	146	}
	147	//____________________________________________________________________
	148	Int_t
	149	AliFMDMCTrackDensity::GetDetectorId() const
	150	{
	151	return AliTrackReference::kFMD;
	152	}
	153
	154	//____________________________________________________________________
	155	void
	156	AliFMDMCTrackDensity::BeginTrackRefs()
	157	{
	158	fState.Clear(true);
f68d9069	159	}
	160
	161	//____________________________________________________________________
	162	void
4bcdcbc1	163	AliFMDMCTrackDensity::EndTrackRefs(Int_t nRefs)
	164	{
	165	fNc->Fill(fState.count);
	166	fNcr->Fill(nRefs, fState.count);
	167	fState.Clear(true);
	168	}
	169
	170	//____________________________________________________________________
	171	AliTrackReference*
	172	AliFMDMCTrackDensity::ProcessRef(AliMCParticle* particle,
	173	const AliMCParticle* mother,
	174	AliTrackReference* ref)
	175	{
	176	// Get the detector coordinates
	177	UShort_t d, s, t;
	178	Char_t r;
	179	AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
	180
	181	// Calculate distance of previous reference to base of cluster
	182	UShort_t nT = TMath::Abs(t - fState.startStrip) + 1;
	183
	184	// Now check if we should flush to output
	185	Bool_t used = false;
	186
	187	// If this is a new detector/ring, then reset the other one
	188	// Check if we have a valid old detectorm ring, and sector
	189	if (fState.oldDetector > 0 &&
	190	fState.oldRing != '\0' &&
	191	fState.oldSector != 1024) {
	192	// New detector, new ring, or new sector
	193	if (d != fState.oldDetector \|\|
	194	r != fState.oldRing \|\|
	195	s != fState.oldSector) {
	196	if (fDebug) Info("Process", "New because new sector");
	197	used = true;
	198	}
	199	else if (nT > fMaxConsequtiveStrips) {
	200	if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)",
	201	fState.nStrips, t, fState.oldStrip, fState.startStrip);
	202	used = true;
	203	}
	204	}
	205	if (used) {
	206	if (fDebug)
	207	Info("Process", "I=%p L=%p D=%d (was %d), R=%c (was %c), "
	208	"S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
	209	ref, fState.longest,
	210	d, fState.oldDetector,
	211	r, fState.oldRing,
	212	s, fState.oldSector,
	213	t, fState.oldStrip,
	214	fState.nStrips, fMaxConsequtiveStrips);
	215	// Int_t nnT = TMath::Abs(fState.oldStrip - fState.startStrip) + 1;
	216	StoreParticle(particle, mother, fState.longest);
	217	fState.Clear(false);
	218	}
	219
	220	// If base of cluster not set, set it here.
	221	if (fState.startStrip == 1024) fState.startStrip = t;
	222
	223	// Calculate distance of previous reference to base of cluster
	224	fState.nStrips = TMath::Abs(t - fState.startStrip) + 1;
	225
	226	// Count number of track refs in this sector
227	fState.nRefs++;
228
229	fState.oldDetector = d;
230	fState.oldRing = r;
231	fState.oldSector = s;
232	fState.oldStrip = t;
233
234	// Debug output
235	if (fDebug) {
236	if (t == fState.startStrip)
237	Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]",
238	d, r, s, t);
239	else
240	Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
241	"length=%3d (now in %3d, previous %3d)",
242	d, r, s, fState.startStrip, fState.nStrips, t, fState.oldStrip);
243	}
244
245	// The longest passage is determined through the angle
246	Double_t ang = GetTrackRefTheta(ref);
247	if (ang > fState.angle) {
248	fState.longest = ref;
249	fState.angle = ang;
250	}
251
252	return fState.longest;
253	}
254
255	//____________________________________________________________________
256	Double_t
f68d9069	257	AliFMDMCTrackDensity::StoreParticle(AliMCParticle* particle,
f68d9069	258	const AliMCParticle* mother,
4bcdcbc1	259	AliTrackReference* ref) const
f68d9069	260	{
4bcdcbc1	261	Double_t w =
	262	AliBaseMCTrackDensity::StoreParticle(particle, mother, ref);
	263	if (w <= 0) return w;
f68d9069	264
f68d9069	265	// Get the detector coordinates
	266	UShort_t d, s, t;
	267	Char_t r;
	268	AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
	269
	270	// Check if we have value already
4bcdcbc1	271	Double_t old = fOutput->Multiplicity(d,r,s,t);
f68d9069	272
	273	// If invalid, force it valid
	274	if (old == AliESDFMD::kInvalidMult) old = 0;
	275
	276	// Increment count
4bcdcbc1	277	fOutput->SetMultiplicity(d,r,s,t,old+w);
f68d9069	278
f68d9069	279	// Fill histograms
4bcdcbc1	280	fNr->Fill(fState.nRefs);
4bcdcbc1	281	fNt->Fill(fState.nStrips);
f68d9069	282
4bcdcbc1	283	fState.count++;
f68d9069	284
4bcdcbc1	285	return w;
f68d9069	286	}
	287
	288	//____________________________________________________________________
	289	Bool_t
	290	AliFMDMCTrackDensity::Calculate(const AliESDFMD& input,
	291	const AliMCEvent& event,
	292	Double_t vz,
	293	AliESDFMD& output,
	294	TH2D* primary)
	295	{
	296	//
	297	// Filter the input kinematics and track references, using
	298	// some of the ESD information
	299	//
	300	// Parameters:
	301	// input Input ESD event
	302	// event Input MC event
	303	// vz Vertex position
	304	// output Output ESD-like object
	305	// primary Per-event histogram of primaries
	306	//
	307	// Return:
	308	// True on succes, false otherwise
	309	//
	310	output.Clear();
4bcdcbc1	311	fOutput = &output;
f68d9069	312
	313	// Copy eta values to output
	314	for (UShort_t ed = 1; ed <= 3; ed++) {
	315	UShort_t nq = (ed == 1 ? 1 : 2);
	316	for (UShort_t eq = 0; eq < nq; eq++) {
	317	Char_t er = (eq == 0 ? 'I' : 'O');
	318	UShort_t ns = (eq == 0 ? 20 : 40);
	319	UShort_t nt = (eq == 0 ? 512 : 256);
	320	for (UShort_t es = 0; es < ns; es++)
	321	for (UShort_t et = 0; et < nt; et++)
	322	output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
	323	}
	324	}
2b556440	325
4bcdcbc1	326	return ProcessTracks(event, vz, primary);
2b556440	327	}
f68d9069	328	//____________________________________________________________________
f68d9069	329	void
4bcdcbc1	330	AliFMDMCTrackDensity::Print(Option_t* option) const
f68d9069	331	{
4bcdcbc1	332	AliBaseMCTrackDensity::Print(option);
f68d9069	333	char ind[gROOT->GetDirLevel()+1];
	334	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
	335	ind[gROOT->GetDirLevel()] = '\0';
4bcdcbc1	336	std::cout << ind << " Max cluster size: " << fMaxConsequtiveStrips
4bcdcbc1	337	<< std::endl;
f68d9069	338
	339	}
	340
	341	//____________________________________________________________________
	342	//
	343	// EOF
	344	//