[u/mrichter/AliRoot.git] / PWG2 / UNICOR / AliDLoop.cxx

//Author: Dariusz Miskowiec 2007

//=============================================================================
// simple event loop manager
// First, events are classified according to centrality, reaction plane angle, 
// and z-vertex. Then, single particles and true pairs are processed. Finally, 
// event mixing is done within event classes. For this, the tree and the event 
// passed as arguments are cloned. 
// The class is, at present, deliberately left in a macro-like design, with all 
// the essential things happening inside the function AliDLoop::Run. This is 
// subject to future polishing if the class proves useful. 
//=============================================================================

#include <TROOT.h>
#include <TMath.h>
#include <TTree.h>
#include <TSystem.h>
#include <TStopwatch.h>
#include "AliDLoop.h"
#include "AliDAnalGlobal.h"
#include "AliDAnalSingle.h"
#include "AliDAnalCorrel.h"
#include "AliDAnalPtfluc.h"
#include "AliDEvent.h"

ClassImp(AliDLoop)

//=============================================================================
AliDLoop::AliDLoop(TTree *tr, AliDEvent *ev0, char *outfil): 
  TObject(), 
  fTree0(tr), fTree1(0x0), 
  fEv0(ev0), fEv1(0x0), 
  fOutputFilename(outfil)     
{
  // constructor
  // Clone the tree such that during the event mixing the two events get 
  // populated from two different trees. This costs some memory (37 MB for 
  // ceres) but allows to avoid switching one tree between two different 
  // events, especially difficult for alice because alice trees do not like 
  // to be attached more than once. 

  fTree1 = (TTree*) tr->Clone();
  fEv1 = (AliDEvent*) ev0->Clone();
  fEv0->AttachTree(fTree0);
  fEv1->AttachTree(fTree1);
}
//=============================================================================
AliDLoop::AliDLoop(const AliDLoop &loop) : 
  TObject(loop), fTree0(loop.fTree0), fTree1(loop.fTree1), fEv0(loop.fEv0), 
  fEv1(loop.fEv1), fOutputFilename(loop.fOutputFilename) 
{     
  //copy constructor, shallow copy, just to make the compiler happy
}
//=============================================================================
AliDLoop &AliDLoop::operator=(const AliDLoop &source) 
{
  // substitution operator, shallow copy, just to make the compiler happy

  TObject::operator=(source);
  fTree0 = source.fTree0;
  fTree1 = source.fTree1;
  fEv0 = source.fEv0;
  fEv1 = source.fEv1;
  fOutputFilename = source.fOutputFilename;
  return *this;
}
//=============================================================================
Int_t AliDLoop::Mem() const 
{
  // return memory in MB occupied by this process

  ProcInfo_t info;
  gSystem->GetProcInfo(&info);
  Int_t memmb = info.fMemVirtual/1024;
  return memmb;
}
//=============================================================================
void AliDLoop::Run(int nwish) const 
{
  // process nwish events from fTree0; nwish=-1 (default) means process all

  // control parameters


  //  const int kcen = 10; // number of centrality bins for mixing
  //  const int kphi = 8; // number of phi bins for mixing
  //  const int kzve = 13; // number of z-vertex bins for mixing
  const int kcen = 1; // number of centrality bins for mixing
  const int kphi = 1; // number of phi bins for mixing
  const int kzve = 1; // number of z-vertex bins for mixing
  int mixingFactor = 1; // number of wished event pairs/number of events
  int mixingStep[10] = {2,3,5,7,11,13,17,19,23,29};
  if (mixingFactor>=10) {printf("mixing factor too high\n"); exit(-1);}
  int minn = 20*mixingFactor; // minimum number of events per bin
  int pr = !gROOT->IsBatch(); // more  printing in interactive mode

  // initialize analysis

  Double_t etami = fEv0->Etamin();
  Double_t etama = fEv0->Etamax();
  AliDAnalGlobal *dag = new AliDAnalGlobal("dag");
  AliDAnalSingle *all = new AliDAnalSingle("all",etami,etama,0);
  AliDAnalSingle *pim = new AliDAnalSingle("pim",etami,etama,-211);
  AliDAnalSingle *pip = new AliDAnalSingle("pip",etami,etama, 211);
  AliDAnalCorrel *cnn = new AliDAnalCorrel("cnn",etami,etama,-211,-211);
  AliDAnalCorrel *cpp = new AliDAnalCorrel("cpp",etami,etama, 211, 211);
  AliDAnalPtfluc *ptf = new AliDAnalPtfluc("ptf",-211,-211);

  // initialize number of events etc

  printf("\ndetermining number of events in chain...");
  int nmax = fTree0->GetEntries();
  int nact = nmax;
  if (nwish>-1 && nwish<nmax) nact = nwish;
  printf("analyzing %d events of %d\n\n",nact,nmax);
  TStopwatch sto;

  sto.Reset();
  sto.Start();

  // find groups of of similar events
  // first allocate nact entries to each array; later reduce

  printf("%4d MB; booking   event index arrays\n",Mem());
  TArrayI *ar[kcen][kphi][kzve];
  int nar[kcen][kphi][kzve];
  for (int i=0; i<kcen; i++) 
    for (int j=0; j<kphi; j++)
      for (int k=0; k<kzve; k++) {
	ar[i][j][k] = new TArrayI(nact);
	nar[i][j][k] = 0;
      }
  printf("%4d MB; filling   event index arrays\n",Mem());
  for (int i=0; i<nact; i++) {
    fTree0->GetEvent(i);
    if (!fEv0->Good()) continue;
    int icen = (int) (kcen*fEv0->Centrality());
    int iphi = (int) (kphi*(fEv0->RPphi()/TMath::Pi()+1)/2);
    int izve = (int) (kzve*(fEv0->Zver()+1.0)/2.0);
    if (pr) printf("\revent %5d   %2d%2d%2d",i,icen,iphi,izve);
    if (icen<0 || icen>=kcen) continue;
    if (iphi<0 || iphi>=kphi) continue;
    if (izve<0 || izve>=kzve) continue;
    ar[icen][iphi][izve]->AddAt(i,nar[icen][iphi][izve]);
    nar[icen][iphi][izve]++;
  }
  if (pr) printf("\n");

  printf("%4d MB; shrinking event index arrays\n",Mem());
  int nall = 0;
  int nsup = 0;
  for (int i=0; i<kcen; i++) 
    for (int j=0; j<kphi; j++)
      for (int k=0; k<kzve; k++) {
	nall += nar[i][j][k];
	// suppress scarse bins
	if (nar[i][j][k]<minn) {
	  nsup += nar[i][j][k];
	  nar[i][j][k]=0;
	}
	ar[i][j][k]->Set(nar[i][j][k]);
      }
  printf("%4d MB; %d events sorted into bins\n",Mem(),nall);
  printf("%d events suppressed because bin occupancy was below %d\n",nsup,minn);
  printf("sorting events took %.1f s CPU %.1f s real %.1f ms real per event\n\n",
	 sto.CpuTime(),sto.RealTime(),1000*sto.RealTime()/nall);
  if ((nall-=nsup) < 1) return;

  // actual loop

  printf("%4d MB; starting the main loop\n",Mem());
  sto.Reset();
  sto.Start();
  for (int i=0; i<kcen; i++) 
    for (int j=0; j<kphi; j++)
      for (int k=0; k<kzve; k++) {

	// true pairs

	int nevents = ar[i][j][k]->GetSize(); // number of events in this bin
	for (int l=0; l<nevents; l++) {
	  if (pr) printf("\revent %6d of %6d  %3d%3d%3d",l,nevents,i,j,k);
	  int m = ar[i][j][k]->At(l); 
	  fTree0->GetEvent(m);
	  dag->Process(fEv0);
	  all->Process(fEv0);
	  pim->Process(fEv0);
	  pip->Process(fEv0);
	  cnn->Process(0,fEv0,fEv0,0);
	  cnn->Process(2,fEv0,fEv0,TMath::DegToRad()*180);	  
	  cpp->Process(0,fEv0,fEv0,0);
	  cpp->Process(2,fEv0,fEv0,TMath::DegToRad()*180);	  
	  ptf->Process(0,fEv0,fEv0);	  
	}
	if (pr && nevents) printf("\n");

	// mixed pairs
	// loop somewhat optimized to reduce jumping back and forth

	for (int imix=0; imix<mixingFactor; imix++) {
	  int step = mixingStep[imix];
	  for (int phase=0; phase<step; phase++) for (int l=phase; l<nevents; l+=step) {
	    int m0 = ar[i][j][k]->At(l); 
	    int next = (l+step)%nevents;
	    int m1 = ar[i][j][k]->At(next); 
	    if (pr) printf("\rmixing %5d and %5d",m0,m1);
	    fTree0->GetEvent(m0);
	    fTree1->GetEvent(m1);
	    //printf("   event0 %f   event1 %f\n",fEv0->ParticleP(0),fEv1->ParticleP(0));
	    cnn->Process(1,fEv0,fEv1,0);
	    cpp->Process(1,fEv0,fEv1,0);
	    ptf->Process(1,fEv0,fEv1);	  
	  }
	if (pr && nevents) printf("\r");
	}
      }

  sto.Stop();
  printf("event loop took %.1f s CPU %.1f s real %.1f ms real per event\n\n",
	 sto.CpuTime(),sto.RealTime(),1000*sto.RealTime()/nall);

  // save analysis results

  dag->Save(fOutputFilename.Data(),"recreate");
  all->Save(fOutputFilename.Data());
  pim->Save(fOutputFilename.Data());
  pip->Save(fOutputFilename.Data());
  cnn->Save(fOutputFilename.Data());
  cpp->Save(fOutputFilename.Data());
  ptf->Save(fOutputFilename.Data());
}
//=============================================================================
Commit	Line	Data
7148817a	1	//Author: Dariusz Miskowiec 2007
	2
	3	//=============================================================================
	4	// simple event loop manager
	5	// First, events are classified according to centrality, reaction plane angle,
	6	// and z-vertex. Then, single particles and true pairs are processed. Finally,
	7	// event mixing is done within event classes. For this, the tree and the event
	8	// passed as arguments are cloned.
	9	// The class is, at present, deliberately left in a macro-like design, with all
	10	// the essential things happening inside the function AliDLoop::Run. This is
	11	// subject to future polishing if the class proves useful.
	12	//=============================================================================
	13
	14	#include <TROOT.h>
5a6d201c	15	#include <TMath.h>
7148817a	16	#include <TTree.h>
	17	#include <TSystem.h>
	18	#include <TStopwatch.h>
	19	#include "AliDLoop.h"
	20	#include "AliDAnalGlobal.h"
	21	#include "AliDAnalSingle.h"
	22	#include "AliDAnalCorrel.h"
	23	#include "AliDAnalPtfluc.h"
	24	#include "AliDEvent.h"
	25
	26	ClassImp(AliDLoop)
	27
	28	//=============================================================================
	29	AliDLoop::AliDLoop(TTree tr, AliDEvent ev0, char *outfil):
	30	TObject(),
	31	fTree0(tr), fTree1(0x0),
	32	fEv0(ev0), fEv1(0x0),
	33	fOutputFilename(outfil)
	34	{
	35	// constructor
	36	// Clone the tree such that during the event mixing the two events get
	37	// populated from two different trees. This costs some memory (37 MB for
	38	// ceres) but allows to avoid switching one tree between two different
	39	// events, especially difficult for alice because alice trees do not like
	40	// to be attached more than once.
	41
	42	fTree1 = (TTree*) tr->Clone();
	43	fEv1 = (AliDEvent*) ev0->Clone();
	44	fEv0->AttachTree(fTree0);
	45	fEv1->AttachTree(fTree1);
	46	}
	47	//=============================================================================
	48	AliDLoop::AliDLoop(const AliDLoop &loop) :
	49	TObject(loop), fTree0(loop.fTree0), fTree1(loop.fTree1), fEv0(loop.fEv0),
	50	fEv1(loop.fEv1), fOutputFilename(loop.fOutputFilename)
	51	{
	52	//copy constructor, shallow copy, just to make the compiler happy
	53	}
	54	//=============================================================================
	55	AliDLoop &AliDLoop::operator=(const AliDLoop &source)
	56	{
	57	// substitution operator, shallow copy, just to make the compiler happy
	58
	59	TObject::operator=(source);
	60	fTree0 = source.fTree0;
	61	fTree1 = source.fTree1;
	62	fEv0 = source.fEv0;
	63	fEv1 = source.fEv1;
	64	fOutputFilename = source.fOutputFilename;
	65	return *this;
	66	}
	67	//=============================================================================
	68	Int_t AliDLoop::Mem() const
	69	{
	70	// return memory in MB occupied by this process
	71
	72	ProcInfo_t info;
	73	gSystem->GetProcInfo(&info);
	74	Int_t memmb = info.fMemVirtual/1024;
	75	return memmb;
	76	}
	77	//=============================================================================
5a6d201c	78	void AliDLoop::Run(int nwish) const
7148817a	79	{
	80	// process nwish events from fTree0; nwish=-1 (default) means process all
	81
	82	// control parameters
	83
	84
	85	// const int kcen = 10; // number of centrality bins for mixing
	86	// const int kphi = 8; // number of phi bins for mixing
	87	// const int kzve = 13; // number of z-vertex bins for mixing
	88	const int kcen = 1; // number of centrality bins for mixing
	89	const int kphi = 1; // number of phi bins for mixing
	90	const int kzve = 1; // number of z-vertex bins for mixing
	91	int mixingFactor = 1; // number of wished event pairs/number of events
	92	int mixingStep[10] = {2,3,5,7,11,13,17,19,23,29};
	93	if (mixingFactor>=10) {printf("mixing factor too high\n"); exit(-1);}
	94	int minn = 20*mixingFactor; // minimum number of events per bin
	95	int pr = !gROOT->IsBatch(); // more printing in interactive mode
	96
	97	// initialize analysis
	98
	99	Double_t etami = fEv0->Etamin();
	100	Double_t etama = fEv0->Etamax();
	101	AliDAnalGlobal *dag = new AliDAnalGlobal("dag");
	102	AliDAnalSingle *all = new AliDAnalSingle("all",etami,etama,0);
	103	AliDAnalSingle *pim = new AliDAnalSingle("pim",etami,etama,-211);
	104	AliDAnalSingle *pip = new AliDAnalSingle("pip",etami,etama, 211);
	105	AliDAnalCorrel *cnn = new AliDAnalCorrel("cnn",etami,etama,-211,-211);
	106	AliDAnalCorrel *cpp = new AliDAnalCorrel("cpp",etami,etama, 211, 211);
	107	AliDAnalPtfluc *ptf = new AliDAnalPtfluc("ptf",-211,-211);
	108
	109	// initialize number of events etc
	110
	111	printf("\ndetermining number of events in chain...");
	112	int nmax = fTree0->GetEntries();
	113	int nact = nmax;
	114	if (nwish>-1 && nwish<nmax) nact = nwish;
	115	printf("analyzing %d events of %d\n\n",nact,nmax);
	116	TStopwatch sto;
	117
	118	sto.Reset();
	119	sto.Start();
	120
	121	// find groups of of similar events
	122	// first allocate nact entries to each array; later reduce
	123
	124	printf("%4d MB; booking event index arrays\n",Mem());
	125	TArrayI *ar[kcen][kphi][kzve];
	126	int nar[kcen][kphi][kzve];
	127	for (int i=0; i<kcen; i++)
	128	for (int j=0; j<kphi; j++)
	129	for (int k=0; k<kzve; k++) {
	130	ar[i][j][k] = new TArrayI(nact);
	131	nar[i][j][k] = 0;
	132	}
	133	printf("%4d MB; filling event index arrays\n",Mem());
	134	for (int i=0; i<nact; i++) {
	135	fTree0->GetEvent(i);
	136	if (!fEv0->Good()) continue;
	137	int icen = (int) (kcen*fEv0->Centrality());
	138	int iphi = (int) (kphi*(fEv0->RPphi()/TMath::Pi()+1)/2);
	139	int izve = (int) (kzve*(fEv0->Zver()+1.0)/2.0);
	140	if (pr) printf("\revent %5d %2d%2d%2d",i,icen,iphi,izve);
	141	if (icen<0 \|\| icen>=kcen) continue;
	142	if (iphi<0 \|\| iphi>=kphi) continue;
143	if (izve<0 \|\| izve>=kzve) continue;
144	ar[icen][iphi][izve]->AddAt(i,nar[icen][iphi][izve]);
145	nar[icen][iphi][izve]++;
146	}
147	if (pr) printf("\n");
148
149	printf("%4d MB; shrinking event index arrays\n",Mem());
150	int nall = 0;
151	int nsup = 0;
152	for (int i=0; i<kcen; i++)
153	for (int j=0; j<kphi; j++)
154	for (int k=0; k<kzve; k++) {
155	nall += nar[i][j][k];
156	// suppress scarse bins
157	if (nar[i][j][k]<minn) {
158	nsup += nar[i][j][k];
159	nar[i][j][k]=0;
160	}
161	ar[i][j][k]->Set(nar[i][j][k]);
162	}
163	printf("%4d MB; %d events sorted into bins\n",Mem(),nall);
164	printf("%d events suppressed because bin occupancy was below %d\n",nsup,minn);
165	printf("sorting events took %.1f s CPU %.1f s real %.1f ms real per event\n\n",
166	sto.CpuTime(),sto.RealTime(),1000*sto.RealTime()/nall);
167	if ((nall-=nsup) < 1) return;
168
169	// actual loop
170
171	printf("%4d MB; starting the main loop\n",Mem());
172	sto.Reset();
173	sto.Start();
174	for (int i=0; i<kcen; i++)
175	for (int j=0; j<kphi; j++)
176	for (int k=0; k<kzve; k++) {
177
178	// true pairs
179
180	int nevents = ar[i][j][k]->GetSize(); // number of events in this bin
181	for (int l=0; l<nevents; l++) {
182	if (pr) printf("\revent %6d of %6d %3d%3d%3d",l,nevents,i,j,k);
183	int m = ar[i][j][k]->At(l);
184	fTree0->GetEvent(m);
185	dag->Process(fEv0);
186	all->Process(fEv0);
187	pim->Process(fEv0);
188	pip->Process(fEv0);
189	cnn->Process(0,fEv0,fEv0,0);
190	cnn->Process(2,fEv0,fEv0,TMath::DegToRad()*180);
191	cpp->Process(0,fEv0,fEv0,0);
192	cpp->Process(2,fEv0,fEv0,TMath::DegToRad()*180);
193	ptf->Process(0,fEv0,fEv0);
194	}
195	if (pr && nevents) printf("\n");
196
197	// mixed pairs
198	// loop somewhat optimized to reduce jumping back and forth
199
200	for (int imix=0; imix<mixingFactor; imix++) {
201	int step = mixingStep[imix];
202	for (int phase=0; phase<step; phase++) for (int l=phase; l<nevents; l+=step) {
203	int m0 = ar[i][j][k]->At(l);
204	int next = (l+step)%nevents;
205	int m1 = ar[i][j][k]->At(next);
206	if (pr) printf("\rmixing %5d and %5d",m0,m1);
207	fTree0->GetEvent(m0);
208	fTree1->GetEvent(m1);
209	//printf(" event0 %f event1 %f\n",fEv0->ParticleP(0),fEv1->ParticleP(0));
210	cnn->Process(1,fEv0,fEv1,0);
211	cpp->Process(1,fEv0,fEv1,0);
212	ptf->Process(1,fEv0,fEv1);
213	}
214	if (pr && nevents) printf("\r");
215	}
216	}
217
218	sto.Stop();
219	printf("event loop took %.1f s CPU %.1f s real %.1f ms real per event\n\n",
220	sto.CpuTime(),sto.RealTime(),1000*sto.RealTime()/nall);
221
222	// save analysis results
223
224	dag->Save(fOutputFilename.Data(),"recreate");
225	all->Save(fOutputFilename.Data());
226	pim->Save(fOutputFilename.Data());
227	pip->Save(fOutputFilename.Data());
228	cnn->Save(fOutputFilename.Data());
229	cpp->Save(fOutputFilename.Data());
230	ptf->Save(fOutputFilename.Data());
231	}
232	//=============================================================================