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