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

#include "AliFMDMCTrackELoss.h"
#include "AliESDFMD.h"
#include "AliTrackReference.h"
#include <TMath.h>
#include "AliFMDStripIndex.h"
#include "AliFMDEncodedEdx.h"
#include "AliLog.h"
#include "AliMCParticle.h"
#include <TH2D.h>
#include <TH1D.h>
#include <TList.h>
#include <TROOT.h>
#include <TTree.h>
#include <TLorentzVector.h>
#include <TParticle.h>
#include <TClonesArray.h>
#include <iostream>

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

//____________________________________________________________________
AliFMDMCTrackELoss::State&
AliFMDMCTrackELoss::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;
  de             = o.de;
  dx             = o.dx;      
  return *this;
}

//====================================================================
AliFMDMCTrackELoss::Hit::Hit() 
  : fGamma(0), fBeta(0), fEta(0), fDe(0), fDx(0), fDetId(0), fPdg(0), 
    fPrimary(true), fDetector(0), fRing('\0'), 
    fSector(0xFFFF), fStrip(0xFFFF)
{}
//____________________________________________________________________
void 
AliFMDMCTrackELoss::Hit::Decode() const
{ 
  if (fDetector > 0) return;
  AliFMDStripIndex::Unpack(fDetId, fDetector, fRing, fSector, fStrip);
}
//____________________________________________________________________
UInt_t
AliFMDMCTrackELoss::Hit::AbsPdg() const 
{ 
  return TMath::Abs(fPdg); 
}

//====================================================================
AliFMDMCTrackELoss::AliFMDMCTrackELoss()
  : AliBaseMCTrackDensity(), 
    fState(),
    fEvent(),
    fMaxConsequtiveStrips(3), 
    fUseTree(false), 
    fHits(0),
    fTree(0),
    fNr(0), 
    fNt(0), 
    fNc(0),
    fNcr(0),
    fBetaGammadEdx(0),
    fBetaGammaEta(0),
    fDEdxEta(0),
    fPrimaries(0),
    fSecondaries(0),
    fAll(0),
    fEta(0)
{
  // Default constructor 
}

//____________________________________________________________________
AliFMDMCTrackELoss::AliFMDMCTrackELoss(const char*)
  : AliBaseMCTrackDensity("fmdMCTrackELoss"), 
    fState(),
    fEvent(),
    fMaxConsequtiveStrips(3), 
    fUseTree(false), 
    fHits(0),
    fTree(0),
    fNr(0), 
    fNt(0), 
    fNc(0),
    fNcr(0),
    fBetaGammadEdx(0),
    fBetaGammaEta(0),
    fDEdxEta(0),
    fPrimaries(0),
    fSecondaries(0),
    fAll(0),
    fEta(0)
{
  // Normal constructor constructor 
  SetTitle("mcTrackELoss");
  fHits = new TClonesArray("AliFMDMCTrackELoss::Hit");
}

#if 0
//____________________________________________________________________
AliFMDMCTrackELoss::AliFMDMCTrackELoss(const AliFMDMCTrackELoss& 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 
}

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

  if (o.fTree) 
    fTree               = static_cast<TTree*>(o.fTree->Clone());

  return *this;
}
#endif

//____________________________________________________________________
void
AliFMDMCTrackELoss::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);
  

  TArrayD edepArray;
  TArrayD bgArray(601);
  AliFMDEncodedEdx::GetdEAxis().FillBinArray(edepArray);
  AliForwardUtil::MakeLogScale(600, -2, 5, bgArray);
  for (Int_t i = 0; i < edepArray.GetSize(); i++) 
    edepArray[i] *= 10;

  fBetaGammadEdx = new TH2D("betaGammadEdx", "Energy loss", 
			    bgArray.GetSize()-1, bgArray.GetArray(),
			    edepArray.GetSize()-1, edepArray.GetArray());
  fBetaGammadEdx->SetDirectory(0);
  fBetaGammadEdx->SetXTitle("#it{#beta#gamma}");
  fBetaGammadEdx->SetYTitle("d#it{#Delta}/d#it{x}");
  fBetaGammadEdx->Sumw2();
  ll->Add(fBetaGammadEdx);
  
  fBetaGammaEta = new TH2D("betaGammaEta", "#beta#gamma", 
			   200, -4, 6, bgArray.GetSize()-1, bgArray.GetArray());
  fBetaGammaEta->SetXTitle("#eta");
  fBetaGammaEta->SetYTitle("#it{#beta#gamma}");
  fBetaGammaEta->Sumw2();
  ll->Add(fBetaGammaEta);

  fDEdxEta = new TH2D("dEdxEta", "d#it{#Delta}/d#it{x}", 
		      200, -4, 6, edepArray.GetSize()-1, edepArray.GetArray());
  fDEdxEta->SetXTitle("#eta");
  fDEdxEta->SetYTitle("d#it{#Delta}/d#it{x}");
  fDEdxEta->Sumw2();
  ll->Add(fDEdxEta);

  if (!fUseTree) return;

  fTree  = new TTree("tree", "Tree of hits");
  fTree->Branch("event", &fEvent, "ipZ/D:cent");
  fTree->Branch("hits", &fHits);
}

//____________________________________________________________________
void
AliFMDMCTrackELoss::Clear(Option_t*)
{
  fPrimaries.Reset(0);
  fSecondaries.Reset(0);
  fAll.Reset(0);
  fEta.Reset(AliESDFMD::kInvalidEta);
  fHits->Clear();
}

//____________________________________________________________________
Int_t
AliFMDMCTrackELoss::GetDetectorId() const
{
  return AliTrackReference::kFMD;
}

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

//____________________________________________________________________
void
AliFMDMCTrackELoss::EndTrackRefs(Int_t nRefs)
{
  fNc->Fill(fState.count);
  fNcr->Fill(nRefs, fState.count);
  fState.Clear(true);
}
  
//____________________________________________________________________
AliTrackReference*
AliFMDMCTrackELoss::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);
  Double_t edep, length, dEdep, dLength;
  AliFMDEncodedEdx::Decode((ref->UserId() >> 19), edep, length, dEdep, dLength);

		 
  // 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;
  fState.de          += edep;
  fState.dx          += length;

  // 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
AliFMDMCTrackELoss::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);
  
  Double_t de   = fState.de;
  Double_t dx   = fState.dx;
  Bool_t   prim = particle == mother;
  fAll(d,r,s,t) += de;
  if (prim) fPrimaries(d,r,s,t)   += de;
  else      fSecondaries(d,r,s,t) += de;

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

  fState.count++;

  if (de <= 0 || dx <= 0) return w;
  
  TLorentzVector v;  
  particle->Particle()->Momentum(v);
  if (v.E() <= 0) return w;
  if (v.Beta() < 0 || v.Beta() > 1) return w;

  Double_t eta   = fEta(d,r,s,t);
  Double_t beta  = v.Beta();
  Double_t gamma = v.Gamma();
  Double_t dEdx  = de/dx;


  fBetaGammadEdx->Fill(beta*gamma, dEdx);
  fBetaGammaEta->Fill(eta, beta*gamma);
  fDEdxEta->Fill(eta, dEdx);

  Int_t nHits = fHits->GetEntries();
  Hit* hit = new((*fHits)[nHits]) Hit;
  hit->fPrimary   = prim;
  hit->fEta       = eta;
  hit->fDe        = de;
  hit->fDx        = dx;
  hit->fDetId     = ref->UserId() & AliFMDStripIndex::kIdMask;
  hit->fBeta      = beta;
  hit->fGamma     = gamma;
  hit->fPdg       = particle->PdgCode();

  return w;
}  

//____________________________________________________________________
Bool_t
AliFMDMCTrackELoss::Calculate(const AliESDFMD&  input, 
			      const AliMCEvent& event,
			      Double_t          vz,
			      Double_t          cent)
{
  // 
  // 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 
  //
  Clear();
  fEvent.fCent = cent;
  fEvent.fIpZ  = vz;

  // 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++) 
	  fEta(ed, er, es, et) = input.Eta(ed, er, es, et);
    }
  }

  Bool_t ret = ProcessTracks(event, vz, 0);

  if (fTree) fTree->Fill();

  return ret;
}

#define PFV(N,VALUE)					\
  do {							\
    AliForwardUtil::PrintName(N);			\
    std::cout << (VALUE) << std::endl; } while(false)
//____________________________________________________________________
void
AliFMDMCTrackELoss::Print(Option_t* option) const 
{
  AliBaseMCTrackDensity::Print(option);
  gROOT->IncreaseDirLevel();
  PFV("Max cluster size", fMaxConsequtiveStrips);
  gROOT->DecreaseDirLevel();
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
7a3e34f4	1	#include "AliFMDMCTrackELoss.h"
	2	#include "AliESDFMD.h"
	3	#include "AliTrackReference.h"
	4	#include <TMath.h>
	5	#include "AliFMDStripIndex.h"
	6	#include "AliFMDEncodedEdx.h"
	7	#include "AliLog.h"
	8	#include "AliMCParticle.h"
	9	#include <TH2D.h>
	10	#include <TH1D.h>
	11	#include <TList.h>
	12	#include <TROOT.h>
	13	#include <TTree.h>
	14	#include <TLorentzVector.h>
	15	#include <TParticle.h>
	16	#include <TClonesArray.h>
	17	#include <iostream>
	18
	19	//====================================================================
	20	void
	21	AliFMDMCTrackELoss::State::Clear(Bool_t alsoCount)
	22	{
	23	angle = 0;
	24	oldDetector = 0;
	25	oldRing = '\0';
	26	oldSector = 1024;
	27	oldStrip = 1024;
	28	startStrip = 1024;
	29	nRefs = 0;
	30	nStrips = 0;
	31	longest = 0x0;
	32	de = 0;
	33	dx = 0;
	34	if (alsoCount) count = 0;
	35	}
	36
	37	//____________________________________________________________________
	38	AliFMDMCTrackELoss::State&
	39	AliFMDMCTrackELoss::State::operator=(const State& o)
	40	{
	41	if (&o == this) return *this;
	42	angle = o.angle;
	43	oldDetector = o.oldDetector;
	44	oldRing = o.oldRing;
	45	oldSector = o.oldSector;
	46	oldStrip = o.oldStrip;
	47	startStrip = o.startStrip;
	48	nRefs = o.nRefs;
	49	nStrips = o.nStrips;
	50	count = o.count;
	51	longest = o.longest;
	52	de = o.de;
	53	dx = o.dx;
	54	return *this;
	55	}
	56
	57	//====================================================================
	58	AliFMDMCTrackELoss::Hit::Hit()
	59	: fGamma(0), fBeta(0), fEta(0), fDe(0), fDx(0), fDetId(0), fPdg(0),
	60	fPrimary(true), fDetector(0), fRing('\0'),
	61	fSector(0xFFFF), fStrip(0xFFFF)
	62	{}
	63	//____________________________________________________________________
	64	void
65	AliFMDMCTrackELoss::Hit::Decode() const
66	{
67	if (fDetector > 0) return;
68	AliFMDStripIndex::Unpack(fDetId, fDetector, fRing, fSector, fStrip);
69	}
70	//____________________________________________________________________
71	UInt_t
72	AliFMDMCTrackELoss::Hit::AbsPdg() const
73	{
74	return TMath::Abs(fPdg);
75	}
76
77	//====================================================================
78	AliFMDMCTrackELoss::AliFMDMCTrackELoss()
79	: AliBaseMCTrackDensity(),
80	fState(),
81	fEvent(),
82	fMaxConsequtiveStrips(3),
83	fUseTree(false),
84	fHits(0),
85	fTree(0),
86	fNr(0),
87	fNt(0),
88	fNc(0),
89	fNcr(0),
90	fBetaGammadEdx(0),
91	fBetaGammaEta(0),
92	fDEdxEta(0),
93	fPrimaries(0),
94	fSecondaries(0),
95	fAll(0),
96	fEta(0)
97	{
98	// Default constructor
99	}
100
101	//____________________________________________________________________
102	AliFMDMCTrackELoss::AliFMDMCTrackELoss(const char*)
103	: AliBaseMCTrackDensity("fmdMCTrackELoss"),
104	fState(),
105	fEvent(),
106	fMaxConsequtiveStrips(3),
107	fUseTree(false),
108	fHits(0),
109	fTree(0),
110	fNr(0),
111	fNt(0),
112	fNc(0),
113	fNcr(0),
114	fBetaGammadEdx(0),
115	fBetaGammaEta(0),
116	fDEdxEta(0),
117	fPrimaries(0),
118	fSecondaries(0),
119	fAll(0),
120	fEta(0)
121	{
122	// Normal constructor constructor
123	SetTitle("mcTrackELoss");
124	fHits = new TClonesArray("AliFMDMCTrackELoss::Hit");
125	}
126
127	#if 0
128	//____________________________________________________________________
129	AliFMDMCTrackELoss::AliFMDMCTrackELoss(const AliFMDMCTrackELoss& o)
130	: AliBaseMCTrackDensity(o),
131	fState(o.fState),
132	fMaxConsequtiveStrips(o.fMaxConsequtiveStrips),
133	fNr(o.fNr),
134	fNt(o.fNt),
135	fNc(o.fNc),
136	fNcr(o.fNcr),
137	fOutput(o.fOutput)
138	{
139	// Normal constructor constructor
140	}
141
142	//____________________________________________________________________
143	AliFMDMCTrackELoss&
144	AliFMDMCTrackELoss::operator=(const AliFMDMCTrackELoss& o)
145	{
146	// Assignment operator
147	if (&o == this) return *this;
148	AliBaseMCTrackDensity::operator=(o);
149	fMaxConsequtiveStrips = o.fMaxConsequtiveStrips;
150	fUseTree = o.fUseTree;
151	fNr = o.fNr;
152	fNt = o.fNt;
153	fNc = o.fNc;
154	fNcr = o.fNcr;
155	fState = o.fState;
156
157	if (o.fTree)
158	fTree = static_cast<TTree*>(o.fTree->Clone());
159
160	return *this;
161	}
162	#endif
163
164	//____________________________________________________________________
165	void
166	AliFMDMCTrackELoss::CreateOutputObjects(TList* l)
167	{
168	AliBaseMCTrackDensity::CreateOutputObjects(l);
169	TList* ll = static_cast<TList*>(l->FindObject(GetTitle()));
170	if (!ll) ll = l;
171
172	fNr = new TH1D("clusterRefs", "# track references per cluster",
173	21, -.5, 20.5);
174	fNr->SetXTitle("# of track references/cluster");
175	fNr->SetFillColor(kRed+1);
176	fNr->SetFillStyle(3001);
177	fNr->SetDirectory(0);
178	ll->Add(fNr);
179
180	fNt = new TH1D("clusterSize", "cluster length in strips", 21, -.5, 20.5);
181	fNt->SetXTitle("Cluster size [strips]");
182	fNt->SetFillColor(kBlue+1);
183	fNt->SetFillStyle(3001);
184	fNt->SetDirectory(0);
185	ll->Add(fNt);
186
187	fNc = new TH1D("nClusters", "# clusters per track", 21, -.5, 20.5);
188	fNc->SetXTitle("# clusters");
189	fNc->SetFillColor(kGreen+1);
190	fNc->SetFillStyle(3001);
191	fNc->SetDirectory(0);
192	ll->Add(fNc);
193
194	fNcr = new TH2D("clusterVsRefs", "# clusters vs. # refs",
195	21, -.5, 20.5, 21, -.5, 20.5);
196	fNcr->SetXTitle("# References");
197	fNcr->SetYTitle("# Clusters");
198	fNcr->SetOption("COLZ");
199	fNcr->SetDirectory(0);
200	ll->Add(fNcr);
201
202
203	TArrayD edepArray;
204	TArrayD bgArray(601);
205	AliFMDEncodedEdx::GetdEAxis().FillBinArray(edepArray);
206	AliForwardUtil::MakeLogScale(600, -2, 5, bgArray);
207	for (Int_t i = 0; i < edepArray.GetSize(); i++)
208	edepArray[i] *= 10;
209
210	fBetaGammadEdx = new TH2D("betaGammadEdx", "Energy loss",
211	bgArray.GetSize()-1, bgArray.GetArray(),
212	edepArray.GetSize()-1, edepArray.GetArray());
213	fBetaGammadEdx->SetDirectory(0);
214	fBetaGammadEdx->SetXTitle("#it{#beta#gamma}");
215	fBetaGammadEdx->SetYTitle("d#it{#Delta}/d#it{x}");
216	fBetaGammadEdx->Sumw2();
217	ll->Add(fBetaGammadEdx);
218
219	fBetaGammaEta = new TH2D("betaGammaEta", "#beta#gamma",
220	200, -4, 6, bgArray.GetSize()-1, bgArray.GetArray());
221	fBetaGammaEta->SetXTitle("#eta");
222	fBetaGammaEta->SetYTitle("#it{#beta#gamma}");
223	fBetaGammaEta->Sumw2();
224	ll->Add(fBetaGammaEta);
225
226	fDEdxEta = new TH2D("dEdxEta", "d#it{#Delta}/d#it{x}",
227	200, -4, 6, edepArray.GetSize()-1, edepArray.GetArray());
228	fDEdxEta->SetXTitle("#eta");
229	fDEdxEta->SetYTitle("d#it{#Delta}/d#it{x}");
230	fDEdxEta->Sumw2();
231	ll->Add(fDEdxEta);
232
233	if (!fUseTree) return;
234
235	fTree = new TTree("tree", "Tree of hits");
236	fTree->Branch("event", &fEvent, "ipZ/D:cent");
237	fTree->Branch("hits", &fHits);
238	}
239
240	//____________________________________________________________________
241	void
242	AliFMDMCTrackELoss::Clear(Option_t*)
243	{
244	fPrimaries.Reset(0);
245	fSecondaries.Reset(0);
246	fAll.Reset(0);
247	fEta.Reset(AliESDFMD::kInvalidEta);
248	fHits->Clear();
249	}
250
251	//____________________________________________________________________
252	Int_t
253	AliFMDMCTrackELoss::GetDetectorId() const
254	{
255	return AliTrackReference::kFMD;
256	}
257
258	//____________________________________________________________________
259	void
260	AliFMDMCTrackELoss::BeginTrackRefs()
261	{
262	fState.Clear(true);
263	}
264
265	//____________________________________________________________________
266	void
267	AliFMDMCTrackELoss::EndTrackRefs(Int_t nRefs)
268	{
269	fNc->Fill(fState.count);
270	fNcr->Fill(nRefs, fState.count);
271	fState.Clear(true);
272	}
273
274	//____________________________________________________________________
275	AliTrackReference*
276	AliFMDMCTrackELoss::ProcessRef(AliMCParticle* particle,
277	const AliMCParticle* mother,
278	AliTrackReference* ref)
279	{
280	// Get the detector coordinates
281	UShort_t d, s, t;
282	Char_t r;
283	AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
284	Double_t edep, length, dEdep, dLength;
285	AliFMDEncodedEdx::Decode((ref->UserId() >> 19), edep, length, dEdep, dLength);
286
287
288
289
290	// Calculate distance of previous reference to base of cluster
291	UShort_t nT = TMath::Abs(t - fState.startStrip) + 1;
292
293	// Now check if we should flush to output
294	Bool_t used = false;
295
296	// If this is a new detector/ring, then reset the other one
297	// Check if we have a valid old detectorm ring, and sector
298	if (fState.oldDetector > 0 &&
299	fState.oldRing != '\0' &&
300	fState.oldSector != 1024) {
301	// New detector, new ring, or new sector
302	if (d != fState.oldDetector \|\|
303	r != fState.oldRing \|\|
304	s != fState.oldSector) {
305	if (fDebug) Info("Process", "New because new sector");
306	used = true;
307	}
308	else if (nT > fMaxConsequtiveStrips) {
309	if (fDebug) Info("Process", "New because too long: %d (%d,%d,%d)",
310	fState.nStrips, t, fState.oldStrip, fState.startStrip);
311	used = true;
312	}
313	}
314	if (used) {
315	if (fDebug)
316	Info("Process", "I=%p L=%p D=%d (was %d), R=%c (was %c), "
317	"S=%2d (was %2d) t=%3d (was %3d) nT=%3d/%4d",
318	ref, fState.longest,
319	d, fState.oldDetector,
320	r, fState.oldRing,
321	s, fState.oldSector,
322	t, fState.oldStrip,
323	fState.nStrips, fMaxConsequtiveStrips);
324	// Int_t nnT = TMath::Abs(fState.oldStrip - fState.startStrip) + 1;
325	StoreParticle(particle, mother, fState.longest);
326	fState.Clear(false);
327	}
328
329	// If base of cluster not set, set it here.
330	if (fState.startStrip == 1024) fState.startStrip = t;
331
332	// Calculate distance of previous reference to base of cluster
333	fState.nStrips = TMath::Abs(t - fState.startStrip) + 1;
334
335	// Count number of track refs in this sector
336	fState.nRefs++;
337
338	fState.oldDetector = d;
339	fState.oldRing = r;
340	fState.oldSector = s;
341	fState.oldStrip = t;
342	fState.de += edep;
343	fState.dx += length;
344
345	// Debug output
346	if (fDebug) {
347	if (t == fState.startStrip)
348	Info("Process", "New cluster starting at FMD%d%c[%2d,%3d]",
349	d, r, s, t);
350	else
351	Info("Process", "Adding to cluster starting at FMD%d%c[%2d,%3d], "
352	"length=%3d (now in %3d, previous %3d)",
353	d, r, s, fState.startStrip, fState.nStrips, t, fState.oldStrip);
354	}
355
356	// The longest passage is determined through the angle
357	Double_t ang = GetTrackRefTheta(ref);
358	if (ang > fState.angle) {
359	fState.longest = ref;
360	fState.angle = ang;
361	}
362
363	return fState.longest;
364	}
365
366	//____________________________________________________________________
367	Double_t
368	AliFMDMCTrackELoss::StoreParticle(AliMCParticle* particle,
369	const AliMCParticle* mother,
370	AliTrackReference* ref) const
371	{
372	Double_t w =
373	AliBaseMCTrackDensity::StoreParticle(particle, mother, ref);
374	if (w <= 0) return w;
375
376	// Get the detector coordinates
377	UShort_t d, s, t;
378	Char_t r;
379	AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
380
381	Double_t de = fState.de;
382	Double_t dx = fState.dx;
383	Bool_t prim = particle == mother;
384	fAll(d,r,s,t) += de;
385	if (prim) fPrimaries(d,r,s,t) += de;
386	else fSecondaries(d,r,s,t) += de;
387
388	// Fill histograms
389	fNr->Fill(fState.nRefs);
390	fNt->Fill(fState.nStrips);
391
392	fState.count++;
393
394	if (de <= 0 \|\| dx <= 0) return w;
395
396	TLorentzVector v;
397	particle->Particle()->Momentum(v);
398	if (v.E() <= 0) return w;
399	if (v.Beta() < 0 \|\| v.Beta() > 1) return w;
400
401	Double_t eta = fEta(d,r,s,t);
402	Double_t beta = v.Beta();
403	Double_t gamma = v.Gamma();
404	Double_t dEdx = de/dx;
405
406
407	fBetaGammadEdx->Fill(beta*gamma, dEdx);
408	fBetaGammaEta->Fill(eta, beta*gamma);
409	fDEdxEta->Fill(eta, dEdx);
410
411	Int_t nHits = fHits->GetEntries();
412	Hit* hit = new((*fHits)[nHits]) Hit;
413	hit->fPrimary = prim;
414	hit->fEta = eta;
415	hit->fDe = de;
416	hit->fDx = dx;
417	hit->fDetId = ref->UserId() & AliFMDStripIndex::kIdMask;
418	hit->fBeta = beta;
419	hit->fGamma = gamma;
420	hit->fPdg = particle->PdgCode();
421
422	return w;
423	}
424
425	//____________________________________________________________________
426	Bool_t
427	AliFMDMCTrackELoss::Calculate(const AliESDFMD& input,
428	const AliMCEvent& event,
429	Double_t vz,
430	Double_t cent)
431	{
432	//
433	// Filter the input kinematics and track references, using
434	// some of the ESD information
435	//
436	// Parameters:
437	// input Input ESD event
438	// event Input MC event
439	// vz Vertex position
440	// output Output ESD-like object
441	// primary Per-event histogram of primaries
442	//
443	// Return:
444	// True on succes, false otherwise
445	//
446	Clear();
447	fEvent.fCent = cent;
448	fEvent.fIpZ = vz;
449
450	// Copy eta values to output
451	for (UShort_t ed = 1; ed <= 3; ed++) {
452	UShort_t nq = (ed == 1 ? 1 : 2);
453	for (UShort_t eq = 0; eq < nq; eq++) {
454	Char_t er = (eq == 0 ? 'I' : 'O');
455	UShort_t ns = (eq == 0 ? 20 : 40);
456	UShort_t nt = (eq == 0 ? 512 : 256);
457	for (UShort_t es = 0; es < ns; es++)
458	for (UShort_t et = 0; et < nt; et++)
459	fEta(ed, er, es, et) = input.Eta(ed, er, es, et);
460	}
461	}
462
463	Bool_t ret = ProcessTracks(event, vz, 0);
464
465	if (fTree) fTree->Fill();
466
467	return ret;
468	}
469
470	#define PFV(N,VALUE) \
471	do { \
472	AliForwardUtil::PrintName(N); \
473	std::cout << (VALUE) << std::endl; } while(false)
474	//____________________________________________________________________
475	void
476	AliFMDMCTrackELoss::Print(Option_t* option) const
477	{
478	AliBaseMCTrackDensity::Print(option);
479	gROOT->IncreaseDirLevel();
480	PFV("Max cluster size", fMaxConsequtiveStrips);
481	gROOT->DecreaseDirLevel();
482	}
483
484	//____________________________________________________________________
485	//
486	// EOF
487	//