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

//
// Class to do the sharing correction of FMD ESD data
//
#include "AliFMDMCSharingFilter.h"
#include <AliESDFMD.h>
#include <AliMCEvent.h>
#include <AliTrackReference.h>
#include <AliStack.h>
#include <TAxis.h>
#include <TList.h>
#include <TH1.h>
#include <TMath.h>
#include "AliForwardCorrectionManager.h"
// #include "AliFMDAnaParameters.h"
#include "AliFMDStripIndex.h"
#include <AliLog.h>
#include <TROOT.h>
#include <iostream>
#include <iomanip>

ClassImp(AliFMDMCSharingFilter)
#if 0
; // This is for Emacs
#endif 


//____________________________________________________________________
AliFMDMCSharingFilter::AliFMDMCSharingFilter()
  : AliFMDSharingFilter(), 
    fFMD1i(0),
    fFMD2i(0),
    fFMD2o(0),
    fFMD3i(0),
    fFMD3o(0), 
    fSumEta(0)
{}

//____________________________________________________________________
AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title)
  : AliFMDSharingFilter(title), 
    fFMD1i(0),
    fFMD2i(0),
    fFMD2o(0),
    fFMD3i(0),
    fFMD3o(0),
    fSumEta(0)
{
  fFMD1i = new TH2D("FMD1i_corr", "Merged vs MC", 21, -.5, 20.5, 100, 0, 20);
  fFMD2i = new TH2D("FMD2i_corr", "Merged vs MC", 21, -.5, 20.5, 100, 0, 20);
  fFMD2o = new TH2D("FMD2o_corr", "Merged vs MC", 21, -.5, 20.5, 100, 0, 20);
  fFMD3i = new TH2D("FMD3i_corr", "Merged vs MC", 21, -.5, 20.5, 100, 0, 20);
  fFMD3o = new TH2D("FMD3o_corr", "Merged vs MC", 21, -.5, 20.5, 100, 0, 20);
  fFMD1i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
  fFMD1i->SetXTitle("Hits (MC)");
  fFMD2i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
  fFMD2i->SetXTitle("Hits (MC)");
  fFMD2o->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
  fFMD2o->SetXTitle("Hits (MC)");
  fFMD3i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
  fFMD3i->SetXTitle("Hits (MC)");
  fFMD3o->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
  fFMD3o->SetXTitle("Hits (MC)");
  fFMD1i->SetDirectory(0);
  fFMD2i->SetDirectory(0);
  fFMD2o->SetDirectory(0);
  fFMD3i->SetDirectory(0);
  fFMD3o->SetDirectory(0);
  fSumEta = new TH1D("mcSumEta", "MC INEL Truth", 200, -4, 6);
  fSumEta->SetXTitle("#eta");
  fSumEta->SetYTitle("dN_{ch}/d#eta");
  fSumEta->SetDirectory(0);
  fSumEta->Sumw2();
  fSumEta->SetMarkerColor(kOrange+2);
  fSumEta->SetMarkerStyle(22);
  fSumEta->SetFillColor(0);
  fSumEta->SetFillStyle(0);
  
}

//____________________________________________________________________
AliFMDMCSharingFilter::AliFMDMCSharingFilter(const AliFMDMCSharingFilter& o)
  : AliFMDSharingFilter(o), 
    fFMD1i(o.fFMD1i),
    fFMD2i(o.fFMD2i),
    fFMD2o(o.fFMD2o),
    fFMD3i(o.fFMD3i),
    fFMD3o(o.fFMD3o),
    fSumEta(o.fSumEta)
{
}

//____________________________________________________________________
AliFMDMCSharingFilter::~AliFMDMCSharingFilter()
{
  if (fFMD1i)  delete fFMD1i;
  if (fFMD2i)  delete fFMD2i;
  if (fFMD2o)  delete fFMD2o;
  if (fFMD3i)  delete fFMD3i;
  if (fFMD3o)  delete fFMD3o;
  if (fSumEta) delete fSumEta;
}

//____________________________________________________________________
AliFMDMCSharingFilter&
AliFMDMCSharingFilter::operator=(const AliFMDMCSharingFilter& o)
{
  AliFMDSharingFilter::operator=(o);
  return *this;
}

//____________________________________________________________________
void
AliFMDMCSharingFilter::StoreParticle(UShort_t   d, 
				     Char_t     r,
				     UShort_t   s, 
				     UShort_t   t, 
				     AliESDFMD& output) const
{
  Double_t old = output.Multiplicity(d,r,s,t);
  if (old == AliESDFMD::kInvalidMult) old = 0;
  output.SetMultiplicity(d,r,s,t,old+1);
}

//____________________________________________________________________
Bool_t
AliFMDMCSharingFilter::FilterMC(const AliESDFMD&  input, 
				const AliMCEvent& event,
				Double_t          vz,
				AliESDFMD&        output)
{
  output.Clear();

  // 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));
    }
  }
  AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
  Int_t nTracks   = event.GetNumberOfTracks();
  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) fSumEta->Fill(particle->Eta());

    Int_t    nTrRef  = particle->GetNumberOfTrackReferences();
    Int_t    longest = -1;
    Double_t angle   = 0;
    UShort_t oD = 0, oS = 1024, oT = 1024;
    Char_t   oR = '\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::kFMD) 
	continue;

      // Get the detector coordinates 
      UShort_t d, s, t;
      Char_t r;
      AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
      if (oD > 0 && oR != '\0' && (d != oD || r != oR)) {
	longest = -1;
	StoreParticle(oD, oR, oS, oT, output);
      }

      oD = d;
      oR = r;
      oS = s;
      oT = t;

      // The longest passage is determined through the angle 
      Double_t x    = ref->X();
      Double_t y    = ref->Y();
      Double_t z    = ref->Z()-vz;
      Double_t rr   = TMath::Sqrt(x*x+y*y);
      Double_t theta= TMath::ATan2(rr,z);
      Double_t ang  = TMath::Abs(TMath::Pi()-theta);
      if (ang > angle) {
	longest = iTrRef;
	angle   = ang;
      }
    } // Loop over track references
    if (longest < 0) continue;

    // Get the reference corresponding to the longest path through the detector
    AliTrackReference* ref = particle->GetTrackReference(longest);

    // Get the detector coordinates 
    UShort_t d, s, t;
    Char_t r;
    AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
    
    StoreParticle(d,r,s,t,output);
  } // Loop over tracks
  return kTRUE;
}

//____________________________________________________________________
void
AliFMDMCSharingFilter::CompareResults(const AliESDFMD&  esd, 
				      const AliESDFMD&  mc)
{
  // Copy eta values to output 
  for (UShort_t d = 1; d <= 3; d++) { 
    UShort_t nq = (d == 1 ? 1 : 2);
    for (UShort_t q = 0; q < nq; q++) {
      Char_t   r  = (q == 0 ? 'I' : 'O');
      UShort_t ns = (q == 0 ?  20 :  40);
      UShort_t nt = (q == 0 ? 512 : 256);
      TH2*     co = 0;
      switch (d) { 
      case 1: co = fFMD1i; break;
      case 2: co = (q == 0 ? fFMD2i : fFMD2o); break;
      case 3: co = (q == 0 ? fFMD3i : fFMD3o); break;
      }

      for (UShort_t s = 0; s < ns; s++) {
	for (UShort_t t = 0; t < nt; t++) { 
	  Float_t mEsd = esd.Multiplicity(d, r, s, t);
	  Float_t mMc  = mc.Multiplicity(d, r, s, t);

	  co->Fill(mMc, mEsd);
	} 
      }
    }
  }
}
  
//____________________________________________________________________
void
AliFMDMCSharingFilter::DefineOutput(TList* dir)
{
  AliFMDSharingFilter::DefineOutput(dir);
  TList* d = static_cast<TList*>(dir->FindObject(GetName()));
  TList* cd = new TList;
  cd->SetName("esd_mc_comparion");
  d->Add(cd);
  cd->Add(fFMD1i);
  cd->Add(fFMD2i);
  cd->Add(fFMD2o);
  cd->Add(fFMD3i);
  cd->Add(fFMD3o);
  dir->Add(fSumEta);
}

//____________________________________________________________________
void
AliFMDMCSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
{
  AliFMDSharingFilter::ScaleHistograms(dir, nEvents);
  TH1D* sumEta = static_cast<TH1D*>(dir->FindObject("mcSumEta"));
  if (!sumEta) { 
    AliWarning(Form("No mcSumEta histogram found in output list"));
    return;
  }
  sumEta->Scale(1. / nEvents, "width");
}

//____________________________________________________________________
void
AliFMDMCSharingFilter::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';
  AliFMDSharingFilter::Print(option);
}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
dc0b1641	1	//
	2	// Class to do the sharing correction of FMD ESD data
	3	//
	4	#include "AliFMDMCSharingFilter.h"
	5	#include <AliESDFMD.h>
	6	#include <AliMCEvent.h>
	7	#include <AliTrackReference.h>
	8	#include <AliStack.h>
	9	#include <TAxis.h>
	10	#include <TList.h>
	11	#include <TH1.h>
	12	#include <TMath.h>
	13	#include "AliForwardCorrectionManager.h"
	14	// #include "AliFMDAnaParameters.h"
	15	#include "AliFMDStripIndex.h"
	16	#include <AliLog.h>
	17	#include <TROOT.h>
	18	#include <iostream>
	19	#include <iomanip>
	20
	21	ClassImp(AliFMDMCSharingFilter)
	22	#if 0
	23	; // This is for Emacs
	24	#endif
	25
	26
	27	//____________________________________________________________________
	28	AliFMDMCSharingFilter::AliFMDMCSharingFilter()
	29	: AliFMDSharingFilter(),
	30	fFMD1i(0),
	31	fFMD2i(0),
	32	fFMD2o(0),
	33	fFMD3i(0),
	34	fFMD3o(0),
	35	fSumEta(0)
	36	{}
	37
	38	//____________________________________________________________________
	39	AliFMDMCSharingFilter::AliFMDMCSharingFilter(const char* title)
	40	: AliFMDSharingFilter(title),
	41	fFMD1i(0),
	42	fFMD2i(0),
	43	fFMD2o(0),
	44	fFMD3i(0),
	45	fFMD3o(0),
	46	fSumEta(0)
	47	{
	48	fFMD1i = new TH2D("FMD1i_corr", "Merged vs MC", 21, -.5, 20.5, 100, 0, 20);
	49	fFMD2i = new TH2D("FMD2i_corr", "Merged vs MC", 21, -.5, 20.5, 100, 0, 20);
	50	fFMD2o = new TH2D("FMD2o_corr", "Merged vs MC", 21, -.5, 20.5, 100, 0, 20);
	51	fFMD3i = new TH2D("FMD3i_corr", "Merged vs MC", 21, -.5, 20.5, 100, 0, 20);
	52	fFMD3o = new TH2D("FMD3o_corr", "Merged vs MC", 21, -.5, 20.5, 100, 0, 20);
	53	fFMD1i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
	54	fFMD1i->SetXTitle("Hits (MC)");
	55	fFMD2i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
	56	fFMD2i->SetXTitle("Hits (MC)");
	57	fFMD2o->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
	58	fFMD2o->SetXTitle("Hits (MC)");
	59	fFMD3i->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
	60	fFMD3i->SetXTitle("Hits (MC)");
	61	fFMD3o->SetYTitle("#Delta E/#Delta_{mip} (ESD)");
	62	fFMD3o->SetXTitle("Hits (MC)");
	63	fFMD1i->SetDirectory(0);
	64	fFMD2i->SetDirectory(0);
65	fFMD2o->SetDirectory(0);
66	fFMD3i->SetDirectory(0);
67	fFMD3o->SetDirectory(0);
68	fSumEta = new TH1D("mcSumEta", "MC INEL Truth", 200, -4, 6);
69	fSumEta->SetXTitle("#eta");
70	fSumEta->SetYTitle("dN_{ch}/d#eta");
71	fSumEta->SetDirectory(0);
72	fSumEta->Sumw2();
73	fSumEta->SetMarkerColor(kOrange+2);
74	fSumEta->SetMarkerStyle(22);
75	fSumEta->SetFillColor(0);
76	fSumEta->SetFillStyle(0);
77
78	}
79
80	//____________________________________________________________________
81	AliFMDMCSharingFilter::AliFMDMCSharingFilter(const AliFMDMCSharingFilter& o)
82	: AliFMDSharingFilter(o),
83	fFMD1i(o.fFMD1i),
84	fFMD2i(o.fFMD2i),
85	fFMD2o(o.fFMD2o),
86	fFMD3i(o.fFMD3i),
87	fFMD3o(o.fFMD3o),
88	fSumEta(o.fSumEta)
89	{
90	}
91
92	//____________________________________________________________________
93	AliFMDMCSharingFilter::~AliFMDMCSharingFilter()
94	{
95	if (fFMD1i) delete fFMD1i;
96	if (fFMD2i) delete fFMD2i;
97	if (fFMD2o) delete fFMD2o;
98	if (fFMD3i) delete fFMD3i;
99	if (fFMD3o) delete fFMD3o;
100	if (fSumEta) delete fSumEta;
101	}
102
103	//____________________________________________________________________
104	AliFMDMCSharingFilter&
105	AliFMDMCSharingFilter::operator=(const AliFMDMCSharingFilter& o)
106	{
107	AliFMDSharingFilter::operator=(o);
108	return *this;
109	}
110
111	//____________________________________________________________________
112	void
113	AliFMDMCSharingFilter::StoreParticle(UShort_t d,
114	Char_t r,
115	UShort_t s,
116	UShort_t t,
117	AliESDFMD& output) const
118	{
119	Double_t old = output.Multiplicity(d,r,s,t);
120	if (old == AliESDFMD::kInvalidMult) old = 0;
121	output.SetMultiplicity(d,r,s,t,old+1);
122	}
123
124	//____________________________________________________________________
125	Bool_t
126	AliFMDMCSharingFilter::FilterMC(const AliESDFMD& input,
127	const AliMCEvent& event,
128	Double_t vz,
129	AliESDFMD& output)
130	{
131	output.Clear();
132
133	// Copy eta values to output
134	for (UShort_t ed = 1; ed <= 3; ed++) {
135	UShort_t nq = (ed == 1 ? 1 : 2);
136	for (UShort_t eq = 0; eq < nq; eq++) {
137	Char_t er = (eq == 0 ? 'I' : 'O');
138	UShort_t ns = (eq == 0 ? 20 : 40);
139	UShort_t nt = (eq == 0 ? 512 : 256);
140	for (UShort_t es = 0; es < ns; es++)
141	for (UShort_t et = 0; et < nt; et++)
142	output.SetEta(ed, er, es, et, input.Eta(ed, er, es, et));
143	}
144	}
145	AliStack* stack = const_cast<AliMCEvent&>(event).Stack();
146	Int_t nTracks = event.GetNumberOfTracks();
147	for (Int_t iTr = 0; iTr < nTracks; iTr++) {
148	AliMCParticle* particle =
149	static_cast<AliMCParticle*>(event.GetTrack(iTr));
150
151	// Check the returned particle
152	if (!particle) continue;
153
154	// Check if this charged and a primary
155	Bool_t isCharged = particle->Charge() != 0;
156	if (!isCharged) continue;
157	Bool_t isPrimary = stack->IsPhysicalPrimary(iTr);
158
159	// Fill 'dn/deta' histogram
160	if (isPrimary) fSumEta->Fill(particle->Eta());
161
162	Int_t nTrRef = particle->GetNumberOfTrackReferences();
163	Int_t longest = -1;
164	Double_t angle = 0;
165	UShort_t oD = 0, oS = 1024, oT = 1024;
166	Char_t oR = '\0';
167	for (Int_t iTrRef = 0; iTrRef < nTrRef; iTrRef++) {
168	AliTrackReference* ref = particle->GetTrackReference(iTrRef);
169
170	// Check existence
171	if (!ref) continue;
172
173	// Check that we hit an FMD element
174	if (ref->DetectorId() != AliTrackReference::kFMD)
175	continue;
176
177	// Get the detector coordinates
178	UShort_t d, s, t;
179	Char_t r;
180	AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
181	if (oD > 0 && oR != '\0' && (d != oD \|\| r != oR)) {
182	longest = -1;
183	StoreParticle(oD, oR, oS, oT, output);
184	}
185
186	oD = d;
187	oR = r;
188	oS = s;
189	oT = t;
190
191	// The longest passage is determined through the angle
192	Double_t x = ref->X();
193	Double_t y = ref->Y();
194	Double_t z = ref->Z()-vz;
195	Double_t rr = TMath::Sqrt(xx+yy);
196	Double_t theta= TMath::ATan2(rr,z);
197	Double_t ang = TMath::Abs(TMath::Pi()-theta);
198	if (ang > angle) {
199	longest = iTrRef;
200	angle = ang;
201	}
202	} // Loop over track references
203	if (longest < 0) continue;
204
205	// Get the reference corresponding to the longest path through the detector
206	AliTrackReference* ref = particle->GetTrackReference(longest);
207
208	// Get the detector coordinates
209	UShort_t d, s, t;
210	Char_t r;
211	AliFMDStripIndex::Unpack(ref->UserId(), d, r, s, t);
212
213	StoreParticle(d,r,s,t,output);
214	} // Loop over tracks
215	return kTRUE;
216	}
217
218	//____________________________________________________________________
219	void
220	AliFMDMCSharingFilter::CompareResults(const AliESDFMD& esd,
221	const AliESDFMD& mc)
222	{
223	// Copy eta values to output
224	for (UShort_t d = 1; d <= 3; d++) {
225	UShort_t nq = (d == 1 ? 1 : 2);
226	for (UShort_t q = 0; q < nq; q++) {
227	Char_t r = (q == 0 ? 'I' : 'O');
228	UShort_t ns = (q == 0 ? 20 : 40);
229	UShort_t nt = (q == 0 ? 512 : 256);
230	TH2* co = 0;
231	switch (d) {
232	case 1: co = fFMD1i; break;
233	case 2: co = (q == 0 ? fFMD2i : fFMD2o); break;
234	case 3: co = (q == 0 ? fFMD3i : fFMD3o); break;
235	}
236
237	for (UShort_t s = 0; s < ns; s++) {
238	for (UShort_t t = 0; t < nt; t++) {
239	Float_t mEsd = esd.Multiplicity(d, r, s, t);
240	Float_t mMc = mc.Multiplicity(d, r, s, t);
241
242	co->Fill(mMc, mEsd);
243	}
244	}
245	}
246	}
247	}
248
249	//____________________________________________________________________
250	void
251	AliFMDMCSharingFilter::DefineOutput(TList* dir)
252	{
253	AliFMDSharingFilter::DefineOutput(dir);
254	TList* d = static_cast<TList*>(dir->FindObject(GetName()));
255	TList* cd = new TList;
256	cd->SetName("esd_mc_comparion");
257	d->Add(cd);
258	cd->Add(fFMD1i);
259	cd->Add(fFMD2i);
260	cd->Add(fFMD2o);
261	cd->Add(fFMD3i);
262	cd->Add(fFMD3o);
263	dir->Add(fSumEta);
264	}
265
266	//____________________________________________________________________
267	void
268	AliFMDMCSharingFilter::ScaleHistograms(TList* dir, Int_t nEvents)
269	{
270	AliFMDSharingFilter::ScaleHistograms(dir, nEvents);
271	TH1D* sumEta = static_cast<TH1D*>(dir->FindObject("mcSumEta"));
272	if (!sumEta) {
273	AliWarning(Form("No mcSumEta histogram found in output list"));
274	return;
275	}
276	sumEta->Scale(1. / nEvents, "width");
277	}
278
279	//____________________________________________________________________
280	void
281	AliFMDMCSharingFilter::Print(Option_t* option) const
282	{
283	char ind[gROOT->GetDirLevel()+1];
284	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
285	ind[gROOT->GetDirLevel()] = '\0';
286	AliFMDSharingFilter::Print(option);
287	}
288
289	//____________________________________________________________________
290	//
291	// EOF
292	//