[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemtoUser / AliFemtoCutMonitorParticleEtCorr.cxx

////////////////////////////////////////////////////////////////////////////////
//                                                                            //
// AliFemtoCutMonitorParticleEtCorr - the cut monitor for particles           //
// which saves particles' et histogram and makes the bin-by-bin correlation   //
//                                                                            //
// Author: Adam.Kisiel@cern.ch                                                //
//                                                                            //
////////////////////////////////////////////////////////////////////////////////
#include "AliFemtoCutMonitorParticleEtCorr.h"
#include "AliFemtoModelHiddenInfo.h"
#include <TH1D.h>
#include <TH2D.h>
#include <TList.h>
#include <TMath.h>

AliFemtoCutMonitorParticleEtCorr::AliFemtoCutMonitorParticleEtCorr():
  AliFemtoCutMonitor(),
  fPhiBins(60),
  fPtPerPhi(0),     
  fPtCovPerPhi(0),
  fPtMultPerPhi(0),
  fNEventsProcessed(0)
{
  // Default constructor
}

AliFemtoCutMonitorParticleEtCorr::AliFemtoCutMonitorParticleEtCorr(const char *aName, int aPhiBins):
  AliFemtoCutMonitor(),
  fPhiBins(aPhiBins),
  fPtPerPhi(0),     
  fPtCovPerPhi(0),
  fPtMultPerPhi(0),
  fNEventsProcessed(0)
{
  // Normal constructor
  char name[200];
  snprintf(name, 200, "EtCorrAvgPt%s", aName);
  fPtPerPhi = new TH1D(name, "Average Pt Per Phi", aPhiBins, -0.5, aPhiBins-0.5);
  snprintf(name, 200, "EtCorrMult%s", aName);
  fPtMultPerPhi = new TH2D(name, "Multiplicity Per Phi", aPhiBins, -0.5, aPhiBins-0.5, aPhiBins, -0.5, aPhiBins-0.5);
  snprintf(name, 200, "EtCorrAvgPtCov%s", aName);
  fPtCovPerPhi = new TH2D(name, "Covariance of Average Pt Per Phi", aPhiBins, -0.5, aPhiBins-0.5, aPhiBins, -0.5, aPhiBins-0.5);

  fPtPerPhi->Sumw2();
  fPtCovPerPhi->Sumw2();
  fPtMultPerPhi->Sumw2();
  fPhiBins = aPhiBins;
}

AliFemtoCutMonitorParticleEtCorr::AliFemtoCutMonitorParticleEtCorr(const AliFemtoCutMonitorParticleEtCorr &aCut):
  AliFemtoCutMonitor(),
  fPhiBins(0),
  fPtPerPhi(0),     
  fPtCovPerPhi(0),
  fPtMultPerPhi(0),
  fNEventsProcessed(0)
{
  // copy constructor
  if (fPtCovPerPhi) delete fPtCovPerPhi;
  fPtCovPerPhi = new TH2D(*aCut.fPtCovPerPhi);
  if (fPtPerPhi) delete fPtPerPhi;
  fPtPerPhi = new TH1D(*aCut.fPtPerPhi);
  if (fPtMultPerPhi) delete fPtMultPerPhi;
  fPtMultPerPhi = new TH2D(*aCut.fPtMultPerPhi);
  fPhiBins = aCut.fPhiBins;
  fNEventsProcessed = aCut.fNEventsProcessed;
}

AliFemtoCutMonitorParticleEtCorr::~AliFemtoCutMonitorParticleEtCorr()
{
  // Destructor
  delete fPtPerPhi;
  delete fPtMultPerPhi;
  delete fPtCovPerPhi;
}

AliFemtoCutMonitorParticleEtCorr& AliFemtoCutMonitorParticleEtCorr::operator=(const AliFemtoCutMonitorParticleEtCorr& aCut)
{
  // assignment operator
  if (this == &aCut) 
    return *this;

  if (fPtCovPerPhi) delete fPtCovPerPhi;
  fPtCovPerPhi = new TH2D(*aCut.fPtCovPerPhi);
  if (fPtPerPhi) delete fPtPerPhi;
  fPtPerPhi = new TH1D(*aCut.fPtPerPhi);
  if (fPtMultPerPhi) delete fPtMultPerPhi;
  fPtMultPerPhi = new TH2D(*aCut.fPtMultPerPhi);
  fPhiBins = aCut.fPhiBins;
  fNEventsProcessed = aCut.fNEventsProcessed;
  
  return *this;
}

AliFemtoString AliFemtoCutMonitorParticleEtCorr::Report(){ 
  // Prepare report from the execution
  string stemp = "*** AliFemtoCutMonitorParticleEtCorr report"; 
  AliFemtoString returnThis = stemp;
  return returnThis; 
}

void AliFemtoCutMonitorParticleEtCorr::Fill(const AliFemtoTrack* aTrack)
{
  // Fill in the monitor histograms with the values from the current track
  //  float tEnergy = ::sqrt(aTrack->P().mag2()+fMass*fMass);
  //  float tRapidity = 0.5*::log((tEnergy+aTrack->P().z())/(tEnergy-aTrack->P().z()));
  float tPt = ::sqrt((aTrack->P().x())*(aTrack->P().x())+(aTrack->P().y())*(aTrack->P().y()));
  //  float tEta = -TMath::Log(TMath::Tan(aTrack->P().theta()/2.0));
  float tPhi = aTrack->P().phi();
  Double_t tPiTwo = TMath::Pi()*2;

  while (tPhi > tPiTwo) tPhi -= tPiTwo;
  while (tPhi < 0) tPhi += tPiTwo;

  int nbin = (int) floor(tPhi * fPhiBins / tPiTwo);
  fPtSumEvent[nbin] += tPt;
  fMultSumEvent[nbin] += 1;
}

void AliFemtoCutMonitorParticleEtCorr::Write()
{
  // Write out the relevant histograms
  fPtPerPhi->Write();
  fPtCovPerPhi->Write();
  fPtMultPerPhi->Write();
}

TList *AliFemtoCutMonitorParticleEtCorr::GetOutputList()
{
  TList *tOutputList = new TList();
  tOutputList->Add(fPtPerPhi);
  tOutputList->Add(fPtCovPerPhi);
  tOutputList->Add(fPtMultPerPhi);

  return tOutputList;
}

void AliFemtoCutMonitorParticleEtCorr::EventBegin(const AliFemtoEvent* aEvent)
{
  for (int iter=0; iter<fPhiBins; iter++) {
    fPtSumEvent[iter] = 0;
    fMultSumEvent[iter] = 0;
  }
}

void AliFemtoCutMonitorParticleEtCorr::EventEnd(const AliFemtoEvent* aEvent)
{
  for (int ispt=0; ispt<fPhiBins; ispt++) {
    fPtPerPhi->Fill(ispt, fPtSumEvent[ispt]);
    for (int ispt2=0; ispt2<fPhiBins; ispt2++) {
      fPtCovPerPhi->Fill(ispt, ispt2, fPtSumEvent[ispt]*fPtSumEvent[ispt2]);
      fPtMultPerPhi->Fill(ispt, ispt2, fMultSumEvent[ispt]*fMultSumEvent[ispt2]);
    }
  }
  fNEventsProcessed++;
}
Commit	Line	Data
95150c4b	1	////////////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoCutMonitorParticleEtCorr - the cut monitor for particles //
	4	// which saves particles' et histogram and makes the bin-by-bin correlation //
	5	// //
	6	// Author: Adam.Kisiel@cern.ch //
	7	// //
	8	////////////////////////////////////////////////////////////////////////////////
	9	#include "AliFemtoCutMonitorParticleEtCorr.h"
	10	#include "AliFemtoModelHiddenInfo.h"
	11	#include <TH1D.h>
	12	#include <TH2D.h>
	13	#include <TList.h>
	14	#include <TMath.h>
	15
	16	AliFemtoCutMonitorParticleEtCorr::AliFemtoCutMonitorParticleEtCorr():
	17	AliFemtoCutMonitor(),
	18	fPhiBins(60),
	19	fPtPerPhi(0),
	20	fPtCovPerPhi(0),
	21	fPtMultPerPhi(0),
	22	fNEventsProcessed(0)
	23	{
	24	// Default constructor
	25	}
	26
	27	AliFemtoCutMonitorParticleEtCorr::AliFemtoCutMonitorParticleEtCorr(const char *aName, int aPhiBins):
	28	AliFemtoCutMonitor(),
	29	fPhiBins(aPhiBins),
	30	fPtPerPhi(0),
	31	fPtCovPerPhi(0),
	32	fPtMultPerPhi(0),
	33	fNEventsProcessed(0)
	34	{
	35	// Normal constructor
	36	char name[200];
	37	snprintf(name, 200, "EtCorrAvgPt%s", aName);
	38	fPtPerPhi = new TH1D(name, "Average Pt Per Phi", aPhiBins, -0.5, aPhiBins-0.5);
	39	snprintf(name, 200, "EtCorrMult%s", aName);
	40	fPtMultPerPhi = new TH2D(name, "Multiplicity Per Phi", aPhiBins, -0.5, aPhiBins-0.5, aPhiBins, -0.5, aPhiBins-0.5);
	41	snprintf(name, 200, "EtCorrAvgPtCov%s", aName);
	42	fPtCovPerPhi = new TH2D(name, "Covariance of Average Pt Per Phi", aPhiBins, -0.5, aPhiBins-0.5, aPhiBins, -0.5, aPhiBins-0.5);
	43
	44	fPtPerPhi->Sumw2();
	45	fPtCovPerPhi->Sumw2();
	46	fPtMultPerPhi->Sumw2();
	47	fPhiBins = aPhiBins;
	48	}
	49
	50	AliFemtoCutMonitorParticleEtCorr::AliFemtoCutMonitorParticleEtCorr(const AliFemtoCutMonitorParticleEtCorr &aCut):
	51	AliFemtoCutMonitor(),
	52	fPhiBins(0),
	53	fPtPerPhi(0),
	54	fPtCovPerPhi(0),
	55	fPtMultPerPhi(0),
	56	fNEventsProcessed(0)
	57	{
	58	// copy constructor
	59	if (fPtCovPerPhi) delete fPtCovPerPhi;
	60	fPtCovPerPhi = new TH2D(*aCut.fPtCovPerPhi);
	61	if (fPtPerPhi) delete fPtPerPhi;
	62	fPtPerPhi = new TH1D(*aCut.fPtPerPhi);
	63	if (fPtMultPerPhi) delete fPtMultPerPhi;
	64	fPtMultPerPhi = new TH2D(*aCut.fPtMultPerPhi);
65	fPhiBins = aCut.fPhiBins;
66	fNEventsProcessed = aCut.fNEventsProcessed;
67	}
68
69	AliFemtoCutMonitorParticleEtCorr::~AliFemtoCutMonitorParticleEtCorr()
70	{
71	// Destructor
72	delete fPtPerPhi;
73	delete fPtMultPerPhi;
74	delete fPtCovPerPhi;
75	}
76
77	AliFemtoCutMonitorParticleEtCorr& AliFemtoCutMonitorParticleEtCorr::operator=(const AliFemtoCutMonitorParticleEtCorr& aCut)
78	{
79	// assignment operator
80	if (this == &aCut)
81	return *this;
82
83	if (fPtCovPerPhi) delete fPtCovPerPhi;
84	fPtCovPerPhi = new TH2D(*aCut.fPtCovPerPhi);
85	if (fPtPerPhi) delete fPtPerPhi;
86	fPtPerPhi = new TH1D(*aCut.fPtPerPhi);
87	if (fPtMultPerPhi) delete fPtMultPerPhi;
88	fPtMultPerPhi = new TH2D(*aCut.fPtMultPerPhi);
89	fPhiBins = aCut.fPhiBins;
90	fNEventsProcessed = aCut.fNEventsProcessed;
91
92	return *this;
93	}
94
95	AliFemtoString AliFemtoCutMonitorParticleEtCorr::Report(){
96	// Prepare report from the execution
97	string stemp = "*** AliFemtoCutMonitorParticleEtCorr report";
98	AliFemtoString returnThis = stemp;
99	return returnThis;
100	}
101
102	void AliFemtoCutMonitorParticleEtCorr::Fill(const AliFemtoTrack* aTrack)
103	{
104	// Fill in the monitor histograms with the values from the current track
105	// float tEnergy = ::sqrt(aTrack->P().mag2()+fMass*fMass);
106	// float tRapidity = 0.5*::log((tEnergy+aTrack->P().z())/(tEnergy-aTrack->P().z()));
107	float tPt = ::sqrt((aTrack->P().x())(aTrack->P().x())+(aTrack->P().y())(aTrack->P().y()));
108	// float tEta = -TMath::Log(TMath::Tan(aTrack->P().theta()/2.0));
109	float tPhi = aTrack->P().phi();
110	Double_t tPiTwo = TMath::Pi()*2;
111
112	while (tPhi > tPiTwo) tPhi -= tPiTwo;
113	while (tPhi < 0) tPhi += tPiTwo;
114
115	int nbin = (int) floor(tPhi * fPhiBins / tPiTwo);
116	fPtSumEvent[nbin] += tPt;
117	fMultSumEvent[nbin] += 1;
118	}
119
120	void AliFemtoCutMonitorParticleEtCorr::Write()
121	{
122	// Write out the relevant histograms
123	fPtPerPhi->Write();
124	fPtCovPerPhi->Write();
125	fPtMultPerPhi->Write();
126	}
127
128	TList *AliFemtoCutMonitorParticleEtCorr::GetOutputList()
129	{
130	TList *tOutputList = new TList();
131	tOutputList->Add(fPtPerPhi);
132	tOutputList->Add(fPtCovPerPhi);
133	tOutputList->Add(fPtMultPerPhi);
134
135	return tOutputList;
136	}
137
138	void AliFemtoCutMonitorParticleEtCorr::EventBegin(const AliFemtoEvent* aEvent)
139	{
140	for (int iter=0; iter<fPhiBins; iter++) {
141	fPtSumEvent[iter] = 0;
142	fMultSumEvent[iter] = 0;
143	}
144	}
145
146	void AliFemtoCutMonitorParticleEtCorr::EventEnd(const AliFemtoEvent* aEvent)
147	{
148	for (int ispt=0; ispt<fPhiBins; ispt++) {
149	fPtPerPhi->Fill(ispt, fPtSumEvent[ispt]);
150	for (int ispt2=0; ispt2<fPhiBins; ispt2++) {
151	fPtCovPerPhi->Fill(ispt, ispt2, fPtSumEvent[ispt]*fPtSumEvent[ispt2]);
152	fPtMultPerPhi->Fill(ispt, ispt2, fMultSumEvent[ispt]*fMultSumEvent[ispt2]);
153	}
154	}
155	fNEventsProcessed++;
156	}
157