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