StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StHbtThreeParticleAnalysis.cxx
1 /***************************************************************************
2  *
3  * $Id: StHbtThreeParticleAnalysis.cxx,v 1.8 2003/09/02 17:58:32 perev Exp $
4  *
5  * Author: Robert Willson, Ohio State, willson@bnl.gov
6  ***************************************************************************
7  *
8  * Description: part of STAR HBT Framework: StHbtMaker package.
9  * This is the derived class for three-particle analysis objects.
10  * Each of the simultaneous analyses should have one of derived
11  * analysis classes running these instantiated. They link
12  * into the aanager in an analysis collection.
13  *
14  ***************************************************************************
15  *
16  * $Log: StHbtThreeParticleAnalysis.cxx,v $
17  * Revision 1.8 2003/09/02 17:58:32 perev
18  * gcc 3.2 updates + WarnOff
19  *
20  * Revision 1.7 2002/08/22 13:14:25 willson
21  * Removed some warnings
22  *
23  * Revision 1.6 2001/06/03 20:56:40 willson
24  * Sectoring, Cos(phi) calculation added
25  *
26  * Revision 1.5 2000/08/11 16:35:41 rcwells
27  * Added number of events processed to each HBT analysis
28  *
29  * Revision 1.4 2000/07/25 03:26:52 willson
30  * Error with small event collections fixed.
31  *
32  * Revision 1.3 2000/05/11 21:18:56 willson
33  * Removed StHbtThreeParticleCorrFctn's...put methods in StHbtCorrFctn
34  * Some methods in derived analysis classes moved to base analysis class
35  *
36  * Revision 1.2 2000/04/12 01:54:20 willson
37  * Initial Installation - Comments Added
38  *
39  *
40  ***************************************************************************/
41 
42 #include "StHbtMaker/Infrastructure/StHbtThreeParticleAnalysis.h"
43 #include "StHbtMaker/Infrastructure/StHbtParticleCollection.hh"
44 #include "StHbtMaker/Base/StHbtTrackCut.h"
45 #include "StHbtMaker/Base/StHbtV0Cut.h"
46 #include "TFile.h"
47 
48 #ifdef __ROOT__
49 ClassImp(StHbtThreeParticleAnalysis)
50 #endif
51 
52 // Factorial Function
53 int Factorial(int arg) {
54  int fact=1;
55  for (int i=arg; i>0; i--) fact*=i;
56  return fact;
57 }
58 
59 //____________________________
60 int StHbtThreeParticleAnalysis::Index(int binx, int biny, int binz)
61 {
62  if (binx==-1 || binx==mNumBinsX || biny==-1 || biny==mNumBinsY || binz==-1 || binz==mNumBinsZ)
63  return -1;
64  return binx + biny*mNumBinsX + binz*mNumBinsX*mNumBinsY;
65 }
66 
67 
68 
69 //This function sorts the event into a number of sectors if mIsSectoring==1;
70 
71 //____________________________
72 void StHbtThreeParticleAnalysis::SortHbtParticleCollection(StHbtParticleCut* partCut,
73  StHbtEvent* hbtEvent,
74  StHbtParticleCollection** SectoredPicoEvent)
75 {
76  int i,j,k;
77  switch (partCut->Type()) {
78  case hbtTrack: // cut is cutting on Tracks
79  {
80  StHbtTrackCut* pCut = (StHbtTrackCut*) partCut;
81  StHbtTrack* pParticle;
82  StHbtTrackIterator pIter;
83  StHbtTrackIterator startLoop = hbtEvent->TrackCollection()->begin();
84  StHbtTrackIterator endLoop = hbtEvent->TrackCollection()->end();
85  float pxave=0.0, pyave=0.0, pzave=0.0, pxhigh=0.0, pxlow=0.0, pyhigh=0.0, pylow=0.0, pzhigh=0.0, pzlow=0.0;
86  int num=0;
87 
88  for (pIter=startLoop;pIter!=endLoop;pIter++){
89  pParticle = *pIter;
90  bool tmpPassParticle = pCut->Pass(pParticle);
91  pCut->FillCutMonitor(pParticle, tmpPassParticle);
92  if (tmpPassParticle){
93  StHbtParticle* particle = new StHbtParticle(pParticle,pCut->Mass());
94  i = (int) floor((particle->FourMomentum().px()-mPXmin)/mDeltaP);
95  j = (int) floor((particle->FourMomentum().py()-mPYmin)/mDeltaP);
96  k = (int) floor((particle->FourMomentum().pz()-mPZmin)/mDeltaP);
97  if (particle->FourMomentum().px() > pxhigh) pxhigh = particle->FourMomentum().px();
98  if (particle->FourMomentum().px() < pxlow) pxlow = particle->FourMomentum().px();
99  if (particle->FourMomentum().py() > pyhigh) pyhigh = particle->FourMomentum().py();
100  if (particle->FourMomentum().py() < pylow) pylow = particle->FourMomentum().py();
101  if (particle->FourMomentum().pz() > pzhigh) pzhigh = particle->FourMomentum().pz();
102  if (particle->FourMomentum().pz() < pzlow) pzlow = particle->FourMomentum().pz();
103  if (i<mNumBinsX && j<mNumBinsY && k<mNumBinsZ && i>=0 && j>=0 && k>=0) {
104  SectoredPicoEvent[Index(i,j,k)]->push_back(particle);
105  pxave+=particle->FourMomentum().px();
106  pyave+=particle->FourMomentum().py();
107  pzave+=particle->FourMomentum().pz();
108  num++;
109  }
110  }
111  }
112  pxave /= num;
113  pyave /= num;
114  pzave /= num;
115 
116  cout << "AVE px=" << pxave << " AVE py="<<pyave<<" AVE pz="<<pzave<< endl;
117  cout << "HIGH px=" << pxhigh << " LOW px="<<pxlow<< " HIGH py=" << pyhigh << " LOW py="<<pylow<<" HIGH pz=" << pzhigh << " LOW pz="<<pzlow<<endl;
118 
119  break;
120  }
121 
122  case hbtV0: // cut is cutting on V0s
123  {
124  StHbtV0Cut* pCut = (StHbtV0Cut*) partCut;
125  StHbtV0* pParticle;
126  StHbtV0Iterator pIter;
127  StHbtV0Iterator startLoop = hbtEvent->V0Collection()->begin();
128  StHbtV0Iterator endLoop = hbtEvent->V0Collection()->end();
129  // this following "for" loop is identical to the one above, but because of scoping, I can's see how to avoid repetition...
130 
131  for (pIter=startLoop;pIter!=endLoop;pIter++){
132  pParticle = *pIter;
133  bool tmpPassV0 = pCut->Pass(pParticle);
134  pCut->FillCutMonitor(pParticle, tmpPassV0);
135  if (tmpPassV0){
136  StHbtParticle* particle = new StHbtParticle(pParticle,pCut->Mass());
137  i = (int) floor((particle->FourMomentum().px()-mPXmin)/mDeltaP);
138  j = (int) floor((particle->FourMomentum().py()-mPYmin)/mDeltaP);
139  k = (int) floor((particle->FourMomentum().pz()-mPZmin)/mDeltaP);
140  SectoredPicoEvent[Index(i,j,k)]->push_back(particle);
141  }
142  }
143  break;
144  }
145  default:
146  cout << "FillHbtParticleCollection function (in StHbtThreeParticleAnalysis.cxx) - undefined Particle Cut type!!! \n";
147  }
148 }
149 
150 // this little function used to apply ParticleCuts (TrackCuts or V0Cuts) and fill
151 // ParticleCollections of picoEvent. It is declared below, and defined in StHbtAnalysis.cxx
152 // It is called from StHbtThreeParticleAnalysis::ProcessEvent() if mIsSectoring==0.
153 void FillHbtParticleCollection(StHbtParticleCut* partCut,
154  StHbtEvent* hbtEvent,
155  StHbtParticleCollection* partCollection);
156 
157 
158 //____________________________
159 StHbtThreeParticleAnalysis::StHbtThreeParticleAnalysis(){
160  // mControlSwitch = 0;
161  mEventCut = 0;
162  mFirstParticleCut = 0;
163  mSecondParticleCut = 0;
164  mThirdParticleCut = 0;
165  mTripletCut = 0;
166  mCorrFctnCollection= 0;
167  mNeventsProcessed = 0;
168  mCalcCosPhi = false;
169  mPXmax = 5.0;
170  mPXmin = -5.0;
171  mPYmax = 5.0;
172  mPYmin = -5.0;
173  mPZmax = 5.0;
174  mPZmin = -5.0;
175  mDeltaP = 1.0;
176  mNumBinsX = 1;
177  mNumBinsY = 1;
178  mNumBinsZ = 1;
179  mNormFactor = 1.0;
180  mCorrFctnCollection = new StHbtCorrFctnCollection;
181  mMixingBuffer = new StHbtPicoEventCollection;
182  mSectoredMixingBuffer = new StHbtSectoredPicoEventCollection;
183  mIsSectoring=0;
184 }
185 
186 //____________________________
187 StHbtThreeParticleAnalysis::~StHbtThreeParticleAnalysis(){
188  //delete mControlSwitch ;
189  delete mEventCut ;
190  delete mFirstParticleCut ;
191  delete mSecondParticleCut ;
192  delete mThirdParticleCut ;
193  delete mTripletCut ;
194  // now delete every CorrFunction in the Collection, and then the Collection itself
195  StHbtCorrFctnIterator iter;
196  for (iter=mCorrFctnCollection->begin(); iter!=mCorrFctnCollection->end();iter++){
197  delete *iter;
198  }
199  delete mCorrFctnCollection;
200  // now delete every PicoEvent in the EventMixingBuffer and then the Buffer itself
201  StHbtPicoEventIterator piter;
202  for (piter=mMixingBuffer->begin();piter!=mMixingBuffer->end();piter++){
203  delete *piter;
204  }
205  delete mMixingBuffer;
206  // now delete every sectored PicoEvent in the SectoredEventMixingBuffer and then the Buffer itself
207  StHbtSectoredPicoEventIterator spiter;
208  for (spiter=mSectoredMixingBuffer->begin();spiter!=mSectoredMixingBuffer->end();spiter++){
209  delete *spiter;
210  }
211  delete mSectoredMixingBuffer;
212 
213  mNeventsProcessed = 0;
214 
215  if (mCalcCosPhi) delete mCosPhiE;
216 }
217 //______________________
218 StHbtCorrFctn* StHbtThreeParticleAnalysis::CorrFctn(int n){ // return pointer to n-th correlation function
219  if ( n<0 || n > (int)mCorrFctnCollection->size() )
220  return NULL;
221  StHbtCorrFctnIterator iter=mCorrFctnCollection->begin();
222  for (int i=0; i<n ;i++){
223  iter++;
224  }
225  return *iter;
226 }
227 //____________________________
228 StHbtString StHbtThreeParticleAnalysis::Report()
229 {
230  char Sect[100];
231  sprintf(Sect, "Sector bounds: %4.2f<px<%4.2f %4.2f<py<%4.2f %4.2f<pz<%4.2f, DeltaP=%4.2f", mPXmin, mPXmax, mPYmin, mPYmax, mPZmin, mPZmax, mDeltaP);
232  if (mCalcCosPhi)
233  cout << "CosPhiAnalysis - constructing Report..."<<endl;
234  else
235  cout << "StHbtThreeParticleAnalysis - constructing Report..."<<endl;
236  string temp = "-----------\nHbt Analysis Report:\n";
237  temp += "\nEvent Cuts:\n";
238  temp += mEventCut->Report();
239  temp += "\nParticle Cuts - First Particle:\n";
240  temp += mFirstParticleCut->Report();
241  temp += "\nParticle Cuts - Second Particle:\n";
242  temp += mSecondParticleCut->Report();
243  temp += "\nParticle Cuts - Third Particle:\n";
244  temp += mThirdParticleCut->Report();
245  temp += "\nTriplet Cuts:\n";
246  temp += mTripletCut->Report();
247  temp += "\nCorrelation Functions:\n";
248  StHbtCorrFctnIterator iter;
249  if ( mCorrFctnCollection->size()==0 ) {
250  cout << "StHbtThreeParticleAnalysis-Warning : no correlations functions in this analysis " << endl;
251  }
252  for (iter=mCorrFctnCollection->begin(); iter!=mCorrFctnCollection->end();iter++){
253  temp += (*iter)->Report();
254  temp += "\n";
255  }
256 
257  if (mIsSectoring) {
258  temp += "\nSectoring Information:\n";
259  temp += Sect;
260  }
261  temp += "-------------\n";
262  StHbtString returnThis=temp;
263  return returnThis;
264 }
265 //_________________________
266 void StHbtThreeParticleAnalysis::ProcessEvent(const StHbtEvent* hbtEvent) {
267 
268  // Add event to processed events
269  AddEventProcessed();
270 
271  if (!mIsSectoring) {
272 
273  if (!mCalcCosPhi) {
274  int NumTriplets=0, NumParticles, NumMixedParticles1, NumMixedParticles2;
275  double TotalRealTriplets, TotalMixedTriplets=0.0, Rate;
276  time_t TotalTime, StartTime;
277  // startup for EbyE
278  EventBegin(hbtEvent);
279  // event cut and event cut monitor
280  bool tmpPassEvent = mEventCut->Pass(hbtEvent);
281  mEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
282  if (tmpPassEvent) {
283  cout << "StHbtThreeParticleAnalysis::ProcessEvent() - Event has passed cut - build picoEvent from " <<
284  hbtEvent->TrackCollection()->size() << " tracks in TrackCollection" << endl;
285  // OK, analysis likes the event-- build a pico event from it, using tracks the analysis likes...
286  StHbtPicoEvent* picoEvent = new StHbtPicoEvent; // this is what we will make Triplets from and put in Mixing Buffer
287  FillHbtParticleCollection(mFirstParticleCut,(StHbtEvent*)hbtEvent,picoEvent->FirstParticleCollection());
288  if ( !(AnalyzeIdenticalParticles()) ) {
289  FillHbtParticleCollection(mSecondParticleCut,(StHbtEvent*)hbtEvent,picoEvent->SecondParticleCollection());
290  FillHbtParticleCollection(mThirdParticleCut,(StHbtEvent*)hbtEvent,picoEvent->ThirdParticleCollection());}
291  cout <<"StHbtThreeParticleAnalysis::ProcessEvent - #particles in First, Second, Third Collections: " <<
292  picoEvent->FirstParticleCollection()->size() << " " <<
293  picoEvent->SecondParticleCollection()->size() << " " <<
294  picoEvent->ThirdParticleCollection()->size() << endl;
295 
296  NumParticles = picoEvent->FirstParticleCollection()->size();
297  TotalRealTriplets = (NumParticles)*(NumParticles-1.0)*(NumParticles-2.0)/6.0;
298 
299  // OK, pico event is built
300  // make real Triplets...
301 
302  // Fabrice points out that we do not need to keep creating/deleting Triplets all the time
303  // We only ever need ONE Triplet, and we can just keep changing internal pointers
304  // this should help speed things up
305  StHbtTriplet* TheTriplet = new StHbtTriplet;
306 
307  StHbtParticleIterator PartIter1;
308  StHbtParticleIterator PartIter2;
309  StHbtParticleIterator PartIter3;
310  StHbtCorrFctnIterator CorrFctnIter;
311  StHbtParticleIterator StartOuterLoop = picoEvent->FirstParticleCollection()->begin(); // always
312  StHbtParticleIterator EndOuterLoop = picoEvent->FirstParticleCollection()->end(); // will be two less if identical
313  StHbtParticleIterator StartMiddleLoop;
314  StHbtParticleIterator EndMiddleLoop;
315  StHbtParticleIterator StartInnerLoop;
316  StHbtParticleIterator EndInnerLoop;
317 
318  if (NumParticles>2) {
319  if (AnalyzeIdenticalParticles()) { // only use First collection
320  EndOuterLoop--;
321  EndOuterLoop--; // outer loop goes to next-to-next-to-last particle in First collection
322  EndMiddleLoop = picoEvent->FirstParticleCollection()->end() ; // middle loop goes to next-to-last particle in First collection
323  EndMiddleLoop--;
324  EndInnerLoop = picoEvent->FirstParticleCollection()->end() ; // inner loop goes to last particle in First collection
325  }
326  else { // nonidentical - loop over First, Second and Third collections
327  StartMiddleLoop = picoEvent->SecondParticleCollection()->begin(); // middle loop starts at first particle in Second collection
328  EndMiddleLoop = picoEvent->SecondParticleCollection()->end() ; // middle loop goes to last particle in Second collection
329  StartInnerLoop = picoEvent->ThirdParticleCollection()->begin(); // inner loop starts at first particle in Third collection
330  EndInnerLoop = picoEvent->ThirdParticleCollection()->end() ; // inner loop goes to last particle in Third collection
331  }
332 
333  StartTime = time(NULL);
334  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
335  if (AnalyzeIdenticalParticles()){
336  StartMiddleLoop = PartIter1;
337  StartMiddleLoop++;
338  }
339  TheTriplet->SetTrack1(*PartIter1);
340  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
341  if (AnalyzeIdenticalParticles()){
342  StartInnerLoop = PartIter2;
343  StartInnerLoop++;
344  }
345  TheTriplet->SetTrack2(*PartIter2);
346  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
347  TheTriplet->SetTrack3(*PartIter3);
348  // The following lines have to be uncommented if you want TripletCutMonitors
349  // they are not in for speed reasons
350  // bool tmpPassTriplet = mTripletCut->Pass(TheTriplet);
351  // mTripletCut->FillCutMonitor(TheTriplet, tmpPassTriplet);
352  // if ( tmpPassTriplet ) {
353  if (mTripletCut->Pass(TheTriplet)){
354  NumTriplets++;
355  for (CorrFctnIter=mCorrFctnCollection->begin();
356  CorrFctnIter!=mCorrFctnCollection->end();CorrFctnIter++){
357  StHbtCorrFctn* CorrFctn = *CorrFctnIter;
358  CorrFctn->AddRealTriplet(TheTriplet);
359  }
360  } // if passed Triplet cut
361  } // loop over third particle
362  } // loop over second particle
363  } // loop over first particle
364  } // Make sure we have at least 3 particles in collection
365 
366  // ok, now make mixed Triplets, if the Mixing buffer is full
367 
368  TotalTime = time(NULL) - StartTime;
369  Rate = (double)NumTriplets/(double)TotalTime;
370 
371  cout << "StHbtThreeParticleAnalysis::ProcessEvent() - reals done, " << NumTriplets << " triplets entered out of " << TotalRealTriplets << "." << endl;
372  cout << "Total time: " << TotalTime << "seconds. This is " << Rate << " triplets per second." << endl;
373  if (MixingBufferFull()){
374  cout << "Mixing Buffer is full - lets rock and roll" << endl;
375  }
376  else {
377  cout << "Mixing Buffer not full -gotta wait " << MixingBuffer()->size() << endl;
378  }
379  NumTriplets=0;
380  if (MixingBufferFull()){
381  StartOuterLoop = picoEvent->FirstParticleCollection()->begin();
382  EndOuterLoop = picoEvent->FirstParticleCollection()->end();
383  StHbtPicoEvent* storedEvent1;
384  StHbtPicoEventIterator picoEventIter1;
385  StHbtPicoEventIterator BeginPicoEvent1;
386  StHbtPicoEventIterator EndPicoEvent1;
387  StHbtPicoEvent* storedEvent2;
388  StHbtPicoEventIterator picoEventIter2;
389  StHbtPicoEventIterator BeginPicoEvent2;
390  StHbtPicoEventIterator EndPicoEvent2;
391 
392 
393  BeginPicoEvent1 = MixingBuffer()->begin();
394  EndPicoEvent1 = MixingBuffer()->end();
395  EndPicoEvent1--;
396  EndPicoEvent2 = MixingBuffer()->end();
397 
398  StartTime = time(NULL);
399  for (picoEventIter1=BeginPicoEvent1;picoEventIter1!=EndPicoEvent1;picoEventIter1++){
400  cout << "Outer Event Done" << endl;
401  BeginPicoEvent2 = picoEventIter1;
402  BeginPicoEvent2++;
403  for (picoEventIter2=BeginPicoEvent2;picoEventIter2!=EndPicoEvent2;picoEventIter2++){
404  storedEvent1 = *picoEventIter1;
405  storedEvent2 = *picoEventIter2;
406  if (AnalyzeIdenticalParticles()){
407  StartMiddleLoop = storedEvent1->FirstParticleCollection()->begin();
408  EndMiddleLoop = storedEvent1->FirstParticleCollection()->end();
409  StartInnerLoop = storedEvent2->FirstParticleCollection()->begin();
410  EndInnerLoop = storedEvent2->FirstParticleCollection()->end();
411  }
412  else{
413  // This is WRONG...Will have to be fixed later.
414  StartMiddleLoop = storedEvent1->SecondParticleCollection()->begin();
415  EndMiddleLoop = storedEvent1->SecondParticleCollection()->end();
416  StartInnerLoop = storedEvent2->ThirdParticleCollection()->begin();
417  EndInnerLoop = storedEvent2->ThirdParticleCollection()->end();
418  }
419 
420  NumMixedParticles1 = storedEvent1->FirstParticleCollection()->size();
421  NumMixedParticles2 = storedEvent2->FirstParticleCollection()->size();
422  TotalMixedTriplets += NumParticles*NumMixedParticles1*NumMixedParticles2;
423 
424  if (NumMixedParticles1>0 && NumMixedParticles2>0) {
425  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
426  TheTriplet->SetTrack1(*PartIter1);
427  for (PartIter2=StartMiddleLoop;PartIter2!=EndMiddleLoop;PartIter2++){
428  TheTriplet->SetTrack2(*PartIter2);
429  for (PartIter3=StartInnerLoop;PartIter3!=EndInnerLoop;PartIter3++){
430  TheTriplet->SetTrack3(*PartIter3);
431  if (mTripletCut->Pass(TheTriplet)){
432  // testing... cout << " TheTriplet passed TripletCut... ";
433  NumTriplets++;
434  for (CorrFctnIter=mCorrFctnCollection->begin();
435  CorrFctnIter!=mCorrFctnCollection->end();CorrFctnIter++){
436  StHbtCorrFctn* CorrFctn = *CorrFctnIter;
437  CorrFctn->AddMixedTriplet(TheTriplet);
438  // testing...cout << " TheTriplet has been added to MixedTriplet method " << endl;
439  }
440  } // if passed Triplet cut
441  } // loop over third particle
442  } // loop over second particle
443  } // loop over first particle
444  } // make sure there are particles in the collections
445  } // loop over subset of pico-events stored in Mixing buffer
446  } // loop over set of pico-events stored in Mixing buffer
447 
448  TotalTime = time(NULL) - StartTime;
449  Rate = NumTriplets/(double)TotalTime;
450  // Now get rid of oldest stored pico-event in buffer.
451  // This means (1) delete the event from memory, (2) "pop" the pointer to it from the MixingBuffer
452  picoEventIter1 = MixingBuffer()->end();
453  picoEventIter1--; // bug fixed malisa 27jul99 - end() is one BEYOND the end! (besides crashing on linux, this was a memory leak)
454  delete *picoEventIter1;
455  MixingBuffer()->pop_back();
456  } // if mixing buffer is full
457  delete TheTriplet;
458  MixingBuffer()->push_front(picoEvent); // store the current pico-event in buffer
459  } // if currentEvent is accepted by currentAnalysis
460  EventEnd(hbtEvent); // cleanup for EbyE
461  cout << "StHbtThreeParticleAnalysis::ProcessEvent() - return to caller ... " << NumTriplets << " triplets out of " << TotalMixedTriplets << " in mixing buffer." << endl;
462  cout << "Total time: " << TotalTime << "seconds. This is " << Rate << " triplets per second." << endl;
463  } // normal correlation function calculation
464 
465  if (mCalcCosPhi) {
466  double C2Q12,C2Q23,C2Q31,C3,arg1,arg2,arg3,arg4,cosphi,cosphiError,termt;
467  int bin1,bin2,bin3,bin4;
468  // event cut and event cut monitor
469  bool tmpPassEvent = mEventCut->Pass(hbtEvent);
470  mEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
471  if (tmpPassEvent) {
472  cout << "StHbtThreeParticleAnalysis::ProcessCosPhi() - Event has passed cut - build picoEvent from " <<
473  hbtEvent->TrackCollection()->size() << " tracks in TrackCollection" << endl;
474  // OK, analysis likes the event-- build a pico event from it, using tracks the analysis likes...
475  StHbtPicoEvent* picoEvent = new StHbtPicoEvent; // this is what we will make Triplets from
476  FillHbtParticleCollection(mFirstParticleCut,(StHbtEvent*)hbtEvent,picoEvent->FirstParticleCollection());
477  if ( !(AnalyzeIdenticalParticles()) ) {
478  FillHbtParticleCollection(mSecondParticleCut,(StHbtEvent*)hbtEvent,picoEvent->SecondParticleCollection());
479  FillHbtParticleCollection(mThirdParticleCut,(StHbtEvent*)hbtEvent,picoEvent->ThirdParticleCollection());}
480  cout <<"StHbtThreeParticleAnalysis::ProcessCosPhi - #particles in First, Second, Third Collections: " <<
481  picoEvent->FirstParticleCollection()->size() << " " <<
482  picoEvent->SecondParticleCollection()->size() << " " <<
483  picoEvent->ThirdParticleCollection()->size() << endl;
484 
485  // OK, pico event is built
486  // make real Triplets...
487 
488  // Fabrice points out that we do not need to keep creating/deleting Triplets all the time
489  // We only ever need ONE Triplet, and we can just keep changing internal pointers
490  // this should help speed things up
491  StHbtTriplet* TheTriplet = new StHbtTriplet;
492 
493  if (picoEvent->FirstParticleCollection()->size()>2) {
494 
495  StHbtParticleIterator PartIter1;
496  StHbtParticleIterator PartIter2;
497  StHbtParticleIterator PartIter3;
498  StHbtParticleIterator StartOuterLoop = picoEvent->FirstParticleCollection()->begin(); // always
499  StHbtParticleIterator EndOuterLoop = picoEvent->FirstParticleCollection()->end(); // will be two less if identical
500  StHbtParticleIterator StartMiddleLoop;
501  StHbtParticleIterator EndMiddleLoop;
502  StHbtParticleIterator StartInnerLoop;
503  StHbtParticleIterator EndInnerLoop;
504  if (AnalyzeIdenticalParticles()) { // only use First collection
505  EndOuterLoop--;
506  EndOuterLoop--; // outer loop goes to next-to-next-to-last particle in First collection
507  EndMiddleLoop = picoEvent->FirstParticleCollection()->end() ; // middle loop goes to next-to-last particle in First collection
508  EndMiddleLoop--;
509  EndInnerLoop = picoEvent->FirstParticleCollection()->end() ; // inner loop goes to last particle in First collection
510  }
511  else { // nonidentical - loop over First, Second and Third collections
512  StartMiddleLoop = picoEvent->SecondParticleCollection()->begin(); // middle loop starts at first particle in Second collection
513  EndMiddleLoop = picoEvent->SecondParticleCollection()->end() ; // middle loop goes to last particle in Second collection
514  StartInnerLoop = picoEvent->ThirdParticleCollection()->begin(); // inner loop starts at first particle in Third collection
515  EndInnerLoop = picoEvent->ThirdParticleCollection()->end() ; // inner loop goes to last particle in Third collection
516  }
517  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
518  if (AnalyzeIdenticalParticles()){
519  StartMiddleLoop = PartIter1;
520  StartMiddleLoop++;
521  }
522  TheTriplet->SetTrack1(*PartIter1);
523  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
524  if (AnalyzeIdenticalParticles()){
525  StartInnerLoop = PartIter2;
526  StartInnerLoop++;
527  }
528  TheTriplet->SetTrack2(*PartIter2);
529  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
530  TheTriplet->SetTrack3(*PartIter3);
531  if (mTripletCut->Pass(TheTriplet)){
532  bin1 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv12());
533  bin2 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv23());
534  bin3 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv31());
535  bin4 = mQ3CF->GetXaxis()->FindBin(TheTriplet->qInv());
536 
537  C2Q12 = mQ2CF->GetBinContent(bin1);
538  C2Q23 = mQ2CF->GetBinContent(bin2);
539  C2Q31 = mQ2CF->GetBinContent(bin3);
540  C3 = mNormFactor*(mQ3CF->GetBinContent(bin4));
541 
542  termt = ::sqrt(fabs((C2Q12-1.0)*(C2Q23-1.0)*(C2Q31-1.0)));
543  if (termt) {
544  cosphi = ((C3-1.0) - (C2Q12-1.0) - (C2Q23-1.0) - (C2Q31-1.0))/(2.0*termt);
545  // cout << "CosPhi = " << cosphi << " Q12 = " << TheTriplet->qInv12() << " Q23 = " << TheTriplet->qInv23() << " Q31 = " << TheTriplet->qInv31() << " Q3 = " << TheTriplet->qInv() << " termt = " << termt << endl;
546  arg1 = mQ3CF->GetBinError(bin4);
547  arg2 = mQ2CF->GetBinError(bin1)*(1.0 + (C2Q23-1.0)*(C2Q31-1.0)*cosphi/termt);
548  arg3 = mQ2CF->GetBinError(bin2)*(1.0 + (C2Q31-1.0)*(C2Q12-1.0)*cosphi/termt);
549  arg4 = mQ2CF->GetBinError(bin3)*(1.0 + (C2Q12-1.0)*(C2Q23-1.0)*cosphi/termt);
550  cosphiError = (1.0/(2.0*termt))*::sqrt(arg1*arg1+arg2*arg2+arg3*arg3+arg4*arg4);
551  }
552  else {
553  cosphi = 0.0;
554  cosphiError = 0.0;
555  }
556 
557  // cout << "Q12:" << TheTriplet->qInv12() << " Q23:" << TheTriplet->qInv23() << " Q31:" << TheTriplet->qInv31() << " bin12:" << bin1 << " bin23:" << bin2 << " bin31:" << bin3 << " Q3:" << TheTriplet->qInv() << "
558  mCosPhi->AddBinContent(bin4, cosphi);
559  mCosPhiE->AddBinContent(bin4, cosphiError);
560  mCosPhiN->AddBinContent(bin4);
561  }// if passed Triplet cut
562  } // loop over third particle
563  } // loop over second particle
564  } // loop over first particle
565 
566  cout << "StHbtThreeParticleAnalysis::ProcessCosPhi() - done" << endl;
567 
568  } // if there are at least 3 particles in our collection
569 
570  delete picoEvent;
571  delete TheTriplet;
572  } // if currentEvent is accepted by currentAnalysis
573  cout << "StHbtThreeParticleAnalysis::ProcessCosPhi() - return to caller ... " << endl;
574  } // Calculate CosPhi
575 
576  } // Non-Sectored Analysis
577 
578 
579  if (mIsSectoring) {
580 
581  if (!mCalcCosPhi) {
582  int NumTriplets=0, NumParticles, i, j, k, i2, j2, k2, i3, j3, k3;
583  int size1=0, size2=0, size3=0;
584  double TotalRealTriplets, Rate;
585  time_t TotalTime, StartTime;
586 
587  // startup for EbyE
588  EventBegin(hbtEvent);
589  // event cut and event cut monitor
590  bool tmpPassEvent = mEventCut->Pass(hbtEvent);
591  mEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
592  if (tmpPassEvent) {
593  cout << "StHbtThreeParticleAnalysis::ProcessEvent() - Event has passed cut - build sectoredpicoEvent from " <<
594  hbtEvent->TrackCollection()->size() << " tracks in TrackCollection" << endl;
595  // OK, analysis likes the event-- build a sectored pico event from it, using tracks the analysis likes...
596 
597  // this is what we will make Triplets from and put in Mixing Buffer
598  StHbtSectoredPicoEvent* picoEvent = new StHbtSectoredPicoEvent(mNumBinsX, mNumBinsY, mNumBinsZ);
599 
600  SortHbtParticleCollection(mFirstParticleCut,(StHbtEvent*)hbtEvent,picoEvent->FirstSectoredCollection());
601  if ( !(AnalyzeIdenticalParticles()) ) {
602  SortHbtParticleCollection(mSecondParticleCut,(StHbtEvent*)hbtEvent,picoEvent->SecondSectoredCollection());
603  SortHbtParticleCollection(mThirdParticleCut,(StHbtEvent*)hbtEvent,picoEvent->ThirdSectoredCollection());
604  }
605 
606  for (k=0; k<mNumBinsZ; k++)
607  for (j=0; j<mNumBinsY; j++)
608  for (i=0; i<mNumBinsX; i++) {
609  size1 += picoEvent->FirstSectoredCollection()[Index(i,j,k)]->size();
610  // cout << "There are " << picoEvent->FirstSectoredCollection()[Index(i,j,k)]->size() << " particles in C["<<i<<"]["<<j<<"]["<<k<<"]"<<endl;
611  if ( !(AnalyzeIdenticalParticles()) ) {
612  size2 += picoEvent->SecondSectoredCollection()[Index(i,j,k)]->size();
613  size3 += picoEvent->ThirdSectoredCollection()[Index(i,j,k)]->size();
614  }
615  }
616 
617  cout <<"StHbtThreeParticleAnalysis::ProcessEvent (Sect) - #particles in First, Second, Third Collections: " <<
618  size1 << " " << size2 << " " << size3 << endl;
619 
620  NumParticles = size1;
621  TotalRealTriplets = (NumParticles)*(NumParticles-1.0)*(NumParticles-2.0)/6.0;
622 
623  StartTime = time(NULL);
624 
625  if (AnalyzeIdenticalParticles()) {
626  // Real Triplets--Identical Particles
627  for (k=0; k<mNumBinsZ; k++)
628  for (j=0; j<mNumBinsY; j++)
629  for (i=0; i<mNumBinsX; i++) {
630 
631  // 71 Terms! This was computed earlier -- Avoids any repetition
632 
633  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j-1,k+1));
634  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j-1,k+1), Index(i-1,j,k+1));
635  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j-1,k+1), Index(i,j-1,k+1));
636  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j-1,k+1), Index(i,j,k+1));
637  //-------------//
638  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j-1,k+1));
639  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j-1,k+1), Index(i-1,j,k+1));
640  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j-1,k+1), Index(i+1,j-1,k+1));
641  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j-1,k+1), Index(i,j,k+1));
642  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j-1,k+1), Index(i+1,j,k));
643  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j-1,k+1), Index(i+1,j,k+1));
644  //-------------//
645  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j-1,k+1));
646  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j-1,k+1), Index(i,j,k+1));
647  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j-1,k+1), Index(i+1,j,k));
648  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j-1,k+1), Index(i+1,j,k+1));
649  //-------------//
650  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j,k+1));
651  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j,k+1), Index(i-1,j+1,k));
652  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j,k+1), Index(i-1,j+1,k+1));
653  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j,k+1), Index(i,j,k+1));
654  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j,k+1), Index(i,j+1,k));
655  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j,k+1), Index(i,j+1,k+1));
656  //-------------//
657  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j,k+1));
658  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i-1,j+1,k));
659  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i-1,j+1,k+1));
660  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i,j+1,k));
661  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i,j+1,k+1));
662  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i+1,j,k));
663  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i+1,j,k+1));
664  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i+1,j+1,k));
665  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i+1,j+1,k+1));
666  //-------------//
667  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j,k+1));
668  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j,k+1), Index(i,j+1,k));
669  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j,k+1), Index(i,j+1,k+1));
670  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j,k+1), Index(i+1,j,k));
671  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j,k+1), Index(i+1,j+1,k));
672  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j,k+1), Index(i+1,j+1,k+1));
673  //-------------//
674  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j+1,k+1));
675  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j+1,k+1), Index(i-1,j+1,k));
676  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j+1,k+1), Index(i,j+1,k));
677  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j+1,k+1), Index(i,j+1,k+1));
678  //-------------//
679  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j+1,k+1));
680  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j+1,k+1), Index(i-1,j+1,k));
681  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j+1,k+1), Index(i,j+1,k));
682  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j+1,k+1), Index(i+1,j,k));
683  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j+1,k+1), Index(i+1,j+1,k));
684  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j+1,k+1), Index(i+1,j+1,k+1));
685  //-------------//
686  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j+1,k+1));
687  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j+1,k+1), Index(i,j+1,k));
688  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j+1,k+1), Index(i+1,j,k));
689  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j+1,k+1), Index(i+1,j+1,k));
690 
691  //-------------//
692  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j+1,k));
693  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i-1,j+1,k), Index(i,j+1,k));
694  //-------------//
695  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j+1,k));
696  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j+1,k), Index(i+1,j,k));
697  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i,j+1,k), Index(i+1,j+1,k));
698  //-------------//
699  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j+1,k));
700  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j+1,k), Index(i+1,j,k));
701 
702  //-------------//
703  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k), Index(i+1,j,k));
704 
705  //Mix sector with itself
706  NumTriplets+=CreateRealTriplets(picoEvent, Index(i,j,k));
707  }
708  cout << "done" << endl;
709  }
710  else {
711  // Real Triplets--Non-Identical Particles
712  for (i=0; i<mNumBinsX; i++)
713  for (j=0; j<mNumBinsY; j++)
714  for (k=0; k<mNumBinsZ; k++)
715  for (i2=i-1; i2<i+1; i2++)
716  for (j2=j-1; j2<j+1; j2++)
717  for (k2=k-1; k2<k+1; k2++)
718  if (i2<mNumBinsX && j2<mNumBinsY && k2<mNumBinsZ && i2!=-1 && j2!=-1 && k2!=-1)
719  for (i3=i-1; i3<i+1; i3++)
720  for (j3=j-1; j3<j+1; j3++)
721  for (k3=k-1; k3<k+1; k3++)
722  if (i3<mNumBinsX && j3<mNumBinsY && k3<mNumBinsZ && i3!=-1 && j3!=-1 && k3!=-1)
723  CreateRealTriplets(picoEvent, Index(i,j,k), Index(i2,j2,k2), Index(i3,j3,k3));
724  }
725 
726  TotalTime = time(NULL) - StartTime;
727  Rate = (double)NumTriplets/(double)TotalTime;
728 
729  cout << "StHbtThreeParticleAnalysis::ProcessEvent() (Sect) - reals done, " << NumTriplets << " triplets entered out of " << (int)TotalRealTriplets << "." << endl;
730  cout << "Total time: " << TotalTime << "seconds. This is " << Rate << " triplets per second." << endl;
731 
732  // ok, now make mixed Triplets, if the Mixing buffer is full
733 
734  if (SectoredMixingBufferFull()){
735  cout << "Mixing Buffer is full - lets rock and roll" << endl;
736  }
737  else {
738  cout << "Mixing Buffer not full -gotta wait " << SectoredMixingBuffer()->size() << endl;
739  }
740  NumTriplets=0;
741  if (SectoredMixingBufferFull()){
742  StHbtSectoredPicoEvent* storedEvent1;
743  StHbtSectoredPicoEventIterator picoEventIter1;
744  StHbtSectoredPicoEventIterator BeginPicoEvent1;
745  StHbtSectoredPicoEventIterator EndPicoEvent1;
746  StHbtSectoredPicoEvent* storedEvent2;
747  StHbtSectoredPicoEventIterator picoEventIter2;
748  StHbtSectoredPicoEventIterator BeginPicoEvent2;
749  StHbtSectoredPicoEventIterator EndPicoEvent2;
750 
751  BeginPicoEvent1 = SectoredMixingBuffer()->begin();
752  EndPicoEvent1 = SectoredMixingBuffer()->end();
753  EndPicoEvent1--;
754  EndPicoEvent2 = SectoredMixingBuffer()->end();
755 
756  StartTime = time(NULL);
757  for (picoEventIter1=BeginPicoEvent1;picoEventIter1!=EndPicoEvent1;picoEventIter1++){
758  BeginPicoEvent2 = picoEventIter1;
759  BeginPicoEvent2++;
760  for (picoEventIter2=BeginPicoEvent2;picoEventIter2!=EndPicoEvent2;picoEventIter2++){
761  storedEvent1 = *picoEventIter1;
762  storedEvent2 = *picoEventIter2;
763  if (AnalyzeIdenticalParticles()){
764  // Mixed Triplets--Identical Particles
765  for (i=0; i<mNumBinsX; i++)
766  for (j=0; j<mNumBinsY; j++)
767  for (k=0; k<mNumBinsZ; k++)
768  for (i2=i-1; i2<=i+1; i2++)
769  for (j2=j-1; j2<=j+1; j2++)
770  for (k2=k-1; k2<=k+1; k2++)
771  if (i2<mNumBinsX && j2<mNumBinsY && k2<mNumBinsZ && i2!=-1 && j2!=-1 && k2!=-1)
772  for (i3=i-1; i3<=i+1; i3++)
773  for (j3=j-1; j3<=j+1; j3++)
774  for (k3=k-1; k3<=k+1; k3++)
775  if (i3<mNumBinsX && j3<mNumBinsY && k3<mNumBinsZ && i3!=-1 && j3!=-1 && k3!=-1 && i3!=i2-2 && i3!=i2+2 && j3!=j2-2 && j3!=j2+2 && k3!=k2-2 && k3!=k2+2)
776  NumTriplets+=CreateMixedTriplets(picoEvent->FirstSectoredCollection()[Index(i,j,k)], storedEvent1->FirstSectoredCollection()[Index(i2,j2,k2)], storedEvent2->FirstSectoredCollection()[Index(i3,j3,k3)]);
777  }
778  else{
779  // This is WRONG...Will have to be fixed later.
780  // Mixed Triplets--Identical Particles
781  for (i=0; i<mNumBinsX; i++)
782  for (j=0; j<mNumBinsY; j++)
783  for (k=0; k<mNumBinsZ; k++)
784  for (i2=i-1; i2<i+1; i2++)
785  for (j2=j-1; j2<j+1; j2++)
786  for (k2=k-1; k2<k+1; k2++)
787  if (i2<mNumBinsX && j2<mNumBinsY && k2<mNumBinsZ && i2!=-1 && j2!=-1 && k2!=-1)
788  for (i3=i-1; i3<i+1; i3++)
789  for (j3=j-1; j3<j+1; j3++)
790  for (k3=k-1; k3<k+1; k3++)
791  if (i3<mNumBinsX && j3<mNumBinsY && k3<mNumBinsZ && i3!=-1 && j3!=-1 && k3!=-1)
792  CreateMixedTriplets(picoEvent->FirstSectoredCollection()[Index(i,j,k)], storedEvent1->SecondSectoredCollection()[Index(i2,j2,k2)], storedEvent2->ThirdSectoredCollection()[Index(i3,j3,k3)]);
793  }
794  }
795  }
796 
797  TotalTime = time(NULL) - StartTime;
798  Rate = NumTriplets/(double)TotalTime;
799 
800  // Now get rid of oldest stored pico-event in buffer.
801  // This means (1) delete the event from memory, (2) "pop" the pointer to it from the MixingBuffer
802  picoEventIter1 = SectoredMixingBuffer()->end();
803  picoEventIter1--; // bug fixed malisa 27jul99 - end() is one BEYOND the end! (besides crashing on linux, this was a memory leak)
804  delete *picoEventIter1;
805 
806  SectoredMixingBuffer()->pop_back();
807  } // if mixing buffer is full
808  SectoredMixingBuffer()->push_front(picoEvent); // store the current pico-event in buffer
809  } // if currentEvent is accepted by currentAnalysis
810  EventEnd(hbtEvent); // cleanup for EbyE
811  cout << "StHbtThreeParticleAnalysis::ProcessEvent() (Sect) - return to caller ... " << NumTriplets << " triplets out of about " << size1*size1*size1 << " in mixing buffer." << endl;
812  cout << "Total time: " << TotalTime << "seconds. This is " << Rate << " triplets per second." << endl;
813 
814  } // normal correlation function calculation
815 
816  if (mCalcCosPhi) {
817  int NumTriplets=0, NumParticles, i, j, k, i2, j2, k2, i3, j3, k3;
818  int size1=0, size2=0, size3=0;
819  double TotalRealTriplets, Rate;
820  time_t TotalTime, StartTime;
821 
822  // event cut and event cut monitor
823  bool tmpPassEvent = mEventCut->Pass(hbtEvent);
824  mEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
825  if (tmpPassEvent) {
826  cout << "StHbtThreeParticleAnalysis::ProcessEvent() - Event has passed cut - build sectoredpicoEvent from " <<
827  hbtEvent->TrackCollection()->size() << " tracks in TrackCollection" << endl;
828  // OK, analysis likes the event-- build a sectored pico event from it, using tracks the analysis likes...
829 
830  // this is what we will make Triplets from and put in Mixing Buffer
831  StHbtSectoredPicoEvent* picoEvent = new StHbtSectoredPicoEvent(mNumBinsX, mNumBinsY, mNumBinsZ);
832 
833  SortHbtParticleCollection(mFirstParticleCut,(StHbtEvent*)hbtEvent,picoEvent->FirstSectoredCollection());
834  if ( !(AnalyzeIdenticalParticles()) ) {
835  SortHbtParticleCollection(mSecondParticleCut,(StHbtEvent*)hbtEvent,picoEvent->SecondSectoredCollection());
836  SortHbtParticleCollection(mThirdParticleCut,(StHbtEvent*)hbtEvent,picoEvent->ThirdSectoredCollection());
837  }
838 
839  for (k=0; k<mNumBinsZ; k++)
840  for (j=0; j<mNumBinsY; j++)
841  for (i=0; i<mNumBinsX; i++) {
842  size1 += picoEvent->FirstSectoredCollection()[Index(i,j,k)]->size();
843  // cout << "There are " << picoEvent->FirstSectoredCollection()[Index(i,j,k)]->size() << " particles in C["<<i<<"]["<<j<<"]["<<k<<"]"<<endl;
844  if ( !(AnalyzeIdenticalParticles()) ) {
845  size2 += picoEvent->SecondSectoredCollection()[Index(i,j,k)]->size();
846  size3 += picoEvent->ThirdSectoredCollection()[Index(i,j,k)]->size();
847  }
848  }
849 
850  cout <<"StHbtThreeParticleAnalysis::ProcessEvent (Sect) - #particles in First, Second, Third Collections: " <<
851  size1 << " " << size2 << " " << size3 << endl;
852 
853  NumParticles = size1;
854  TotalRealTriplets = (NumParticles)*(NumParticles-1.0)*(NumParticles-2.0)/6.0;
855 
856  StartTime = time(NULL);
857 
858  if (AnalyzeIdenticalParticles()) {
859  // Real Triplets--Identical Particles
860  for (k=0; k<mNumBinsZ; k++)
861  for (j=0; j<mNumBinsY; j++)
862  for (i=0; i<mNumBinsX; i++) {
863 
864  // 71 Terms! This was computed earlier -- Avoids any repetition
865 
866  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j-1,k+1));
867  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j-1,k+1), Index(i-1,j,k+1));
868  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j-1,k+1), Index(i,j-1,k+1));
869  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j-1,k+1), Index(i,j,k+1));
870  //-------------//
871  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j-1,k+1));
872  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j-1,k+1), Index(i-1,j,k+1));
873  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j-1,k+1), Index(i+1,j-1,k+1));
874  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j-1,k+1), Index(i,j,k+1));
875  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j-1,k+1), Index(i+1,j,k));
876  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j-1,k+1), Index(i+1,j,k+1));
877  //-------------//
878  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j-1,k+1));
879  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j-1,k+1), Index(i,j,k+1));
880  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j-1,k+1), Index(i+1,j,k));
881  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j-1,k+1), Index(i+1,j,k+1));
882  //-------------//
883  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j,k+1));
884  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j,k+1), Index(i-1,j+1,k));
885  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j,k+1), Index(i-1,j+1,k+1));
886  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j,k+1), Index(i,j,k+1));
887  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j,k+1), Index(i,j+1,k));
888  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j,k+1), Index(i,j+1,k+1));
889  //-------------//
890  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j,k+1));
891  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i-1,j+1,k));
892  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i-1,j+1,k+1));
893  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i,j+1,k));
894  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i,j+1,k+1));
895  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i+1,j,k));
896  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i+1,j,k+1));
897  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i+1,j+1,k));
898  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j,k+1), Index(i+1,j+1,k+1));
899  //-------------//
900  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j,k+1));
901  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j,k+1), Index(i,j+1,k));
902  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j,k+1), Index(i,j+1,k+1));
903  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j,k+1), Index(i+1,j,k));
904  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j,k+1), Index(i+1,j+1,k));
905  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j,k+1), Index(i+1,j+1,k+1));
906  //-------------//
907  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j+1,k+1));
908  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j+1,k+1), Index(i-1,j+1,k));
909  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j+1,k+1), Index(i,j+1,k));
910  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j+1,k+1), Index(i,j+1,k+1));
911  //-------------//
912  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j+1,k+1));
913  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j+1,k+1), Index(i-1,j+1,k));
914  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j+1,k+1), Index(i,j+1,k));
915  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j+1,k+1), Index(i+1,j,k));
916  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j+1,k+1), Index(i+1,j+1,k));
917  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j+1,k+1), Index(i+1,j+1,k+1));
918  //-------------//
919  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j+1,k+1));
920  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j+1,k+1), Index(i,j+1,k));
921  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j+1,k+1), Index(i+1,j,k));
922  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j+1,k+1), Index(i+1,j+1,k));
923 
924  //-------------//
925  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j+1,k));
926  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i-1,j+1,k), Index(i,j+1,k));
927  //-------------//
928  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j+1,k));
929  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j+1,k), Index(i+1,j,k));
930  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i,j+1,k), Index(i+1,j+1,k));
931  //-------------//
932  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j+1,k));
933  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j+1,k), Index(i+1,j,k));
934 
935  //-------------//
936  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k), Index(i+1,j,k));
937 
938  //Mix sector with itself
939  NumTriplets+=CalculateCosPhi(picoEvent, Index(i,j,k));
940  }
941  cout << "done" << endl;
942  }
943  else {
944  // Real Triplets--Non-Identical Particles
945  for (i=0; i<mNumBinsX; i++)
946  for (j=0; j<mNumBinsY; j++)
947  for (k=0; k<mNumBinsZ; k++)
948  for (i2=i-1; i2<i+1; i2++)
949  for (j2=j-1; j2<j+1; j2++)
950  for (k2=k-1; k2<k+1; k2++)
951  if (i2<mNumBinsX && j2<mNumBinsY && k2<mNumBinsZ && i2!=-1 && j2!=-1 && k2!=-1)
952  for (i3=i-1; i3<i+1; i3++)
953  for (j3=j-1; j3<j+1; j3++)
954  for (k3=k-1; k3<k+1; k3++)
955  if (i3<mNumBinsX && j3<mNumBinsY && k3<mNumBinsZ && i3!=-1 && j3!=-1 && k3!=-1)
956  CalculateCosPhi(picoEvent, Index(i,j,k), Index(i2,j2,k2), Index(i3,j3,k3));
957  }
958 
959  TotalTime = time(NULL) - StartTime;
960  Rate = (double)NumTriplets/(double)TotalTime;
961 
962  cout << "StHbtThreeParticleAnalysis::ProcessCosPhi() (Sect) - done, " << NumTriplets << " triplets entered out of " << (int)TotalRealTriplets << "." << endl;
963  cout << "Total time: " << TotalTime << "seconds. This is " << Rate << " triplets per second." << endl;
964 
965  delete picoEvent;
966 
967  } // if currentEvent is accepted by currentAnalysis
968  cout << "StHbtThreeParticleAnalysis::ProcessCosPhi() (Sect) - return to caller ... " << endl;
969  } // Calculate CosPhi
970 
971  } // Sectoring Analysis
972 
973 } // ProcessEvent
974 
975 //_________________________
976 int StHbtThreeParticleAnalysis::CreateRealTriplets(StHbtSectoredPicoEvent *picoEvent, int Index1) {
977 
978  int NumTriplets=0;
979 
980  StHbtParticleCollection *partColl1 = picoEvent->FirstSectoredCollection()[Index1];
981 
982  if (partColl1->size()<3) return 0;
983 
984  StHbtTriplet* TheTriplet = new StHbtTriplet;
985 
986  StHbtParticleIterator PartIter1;
987  StHbtParticleIterator PartIter2;
988  StHbtParticleIterator PartIter3;
989  StHbtCorrFctnIterator CorrFctnIter;
990  StHbtParticleIterator StartOuterLoop = partColl1->begin(); // always
991  StHbtParticleIterator EndOuterLoop = partColl1->end(); // will be two less
992  StHbtParticleIterator StartMiddleLoop;
993  StHbtParticleIterator EndMiddleLoop;
994  StHbtParticleIterator StartInnerLoop;
995  StHbtParticleIterator EndInnerLoop;
996 
997  EndOuterLoop--;
998  EndOuterLoop--; // outer loop goes to next-to-next-to-last particle in collection
999  EndMiddleLoop = partColl1->end() ; // middle loop goes to next-to-last particle in collection
1000  EndMiddleLoop--;
1001  EndInnerLoop = partColl1->end() ; // inner loop goes to last particle in collection
1002 
1003  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
1004  StartMiddleLoop = PartIter1;
1005  StartMiddleLoop++;
1006  TheTriplet->SetTrack1(*PartIter1);
1007  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
1008  StartInnerLoop = PartIter2;
1009  StartInnerLoop++;
1010  TheTriplet->SetTrack2(*PartIter2);
1011  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
1012  TheTriplet->SetTrack3(*PartIter3);
1013  // The following lines have to be uncommented if you want TripletCutMonitors
1014  // they are not in for speed reasons
1015  // bool tmpPassTriplet = mTripletCut->Pass(TheTriplet);
1016  // mTripletCut->FillCutMonitor(TheTriplet, tmpPassTriplet);
1017  // if ( tmpPassTriplet ) {
1018  if (mTripletCut->Pass(TheTriplet)){
1019  NumTriplets++;
1020  for (CorrFctnIter=mCorrFctnCollection->begin();
1021  CorrFctnIter!=mCorrFctnCollection->end();CorrFctnIter++){
1022  StHbtCorrFctn* CorrFctn = *CorrFctnIter;
1023  CorrFctn->AddRealTriplet(TheTriplet);
1024  }
1025  } // if passed Triplet cut
1026  } // loop over third particle
1027  } // loop over second particle
1028  } // loop over first particle
1029 
1030  delete TheTriplet;
1031 
1032  return NumTriplets;
1033 
1034 }
1035 
1036 
1037 
1038 //_________________________
1039 int StHbtThreeParticleAnalysis::CreateRealTriplets(StHbtSectoredPicoEvent *picoEvent, int Index1, int Index2) {
1040 
1041  int NumTriplets=0;
1042 
1043  if (Index1==-1 || Index2==-1) return 0;
1044 
1045  StHbtParticleCollection *partColl1 = picoEvent->FirstSectoredCollection()[Index1];
1046  StHbtParticleCollection *partColl2 = picoEvent->FirstSectoredCollection()[Index2];
1047 
1048  StHbtTriplet* TheTriplet = new StHbtTriplet;
1049 
1050  // First, take two particles from the first collection and one particle from the second.
1051 
1052  if (partColl1->size()>=2 && partColl2->size()>=1) {
1053 
1054  StHbtParticleIterator PartIter1;
1055  StHbtParticleIterator PartIter2;
1056  StHbtParticleIterator PartIter3;
1057  StHbtCorrFctnIterator CorrFctnIter;
1058  StHbtParticleIterator StartOuterLoop = partColl1->begin(); // always
1059  StHbtParticleIterator EndOuterLoop = partColl1->end(); // will be one less
1060  StHbtParticleIterator StartMiddleLoop; // depends on start value of outer loop
1061  StHbtParticleIterator EndMiddleLoop = partColl1->end();
1062  StHbtParticleIterator StartInnerLoop = partColl2->begin();
1063  StHbtParticleIterator EndInnerLoop = partColl2->end();
1064 
1065  EndOuterLoop--; // outer loop goes to next-to-last particle in First collection
1066 
1067  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
1068  StartMiddleLoop = PartIter1;
1069  StartMiddleLoop++;
1070  TheTriplet->SetTrack1(*PartIter1);
1071  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
1072  TheTriplet->SetTrack2(*PartIter2);
1073  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
1074  TheTriplet->SetTrack3(*PartIter3);
1075  // The following lines have to be uncommented if you want TripletCutMonitors
1076  // they are not in for speed reasons
1077  // bool tmpPassTriplet = mTripletCut->Pass(TheTriplet);
1078  // mTripletCut->FillCutMonitor(TheTriplet, tmpPassTriplet);
1079  // if ( tmpPassTriplet ) {
1080  if (mTripletCut->Pass(TheTriplet)){
1081  NumTriplets++;
1082  for (CorrFctnIter=mCorrFctnCollection->begin();
1083  CorrFctnIter!=mCorrFctnCollection->end();CorrFctnIter++){
1084  StHbtCorrFctn* CorrFctn = *CorrFctnIter;
1085  CorrFctn->AddRealTriplet(TheTriplet);
1086  }
1087  } // if passed Triplet cut
1088  } // loop over third particle
1089  } // loop over second particle
1090  } // loop over first particle
1091 
1092  }
1093 
1094 
1095  // Then, take two particles from the second collection and one particle from the first.
1096 
1097  if (partColl1->size()>=1 && partColl2->size()>=2) {
1098 
1099  StHbtParticleIterator PartIter1;
1100  StHbtParticleIterator PartIter2;
1101  StHbtParticleIterator PartIter3;
1102  StHbtCorrFctnIterator CorrFctnIter;
1103  StHbtParticleIterator StartOuterLoop = partColl2->begin(); // always
1104  StHbtParticleIterator EndOuterLoop = partColl2->end(); // will be one less
1105  StHbtParticleIterator StartMiddleLoop; // depends on start value of outer loop
1106  StHbtParticleIterator EndMiddleLoop = partColl2->end();
1107  StHbtParticleIterator StartInnerLoop = partColl1->begin();
1108  StHbtParticleIterator EndInnerLoop = partColl1->end();
1109 
1110  EndOuterLoop--; // outer loop goes to next-to-last particle in First collection
1111 
1112  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
1113  StartMiddleLoop = PartIter1;
1114  StartMiddleLoop++;
1115  TheTriplet->SetTrack1(*PartIter1);
1116  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
1117  TheTriplet->SetTrack2(*PartIter2);
1118  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
1119  TheTriplet->SetTrack3(*PartIter3);
1120  // The following lines have to be uncommented if you want TripletCutMonitors
1121  // they are not in for speed reasons
1122  // bool tmpPassTriplet = mTripletCut->Pass(TheTriplet);
1123  // mTripletCut->FillCutMonitor(TheTriplet, tmpPassTriplet);
1124  // if ( tmpPassTriplet ) {
1125  if (mTripletCut->Pass(TheTriplet)){
1126  NumTriplets++;
1127  for (CorrFctnIter=mCorrFctnCollection->begin();
1128  CorrFctnIter!=mCorrFctnCollection->end();CorrFctnIter++){
1129  StHbtCorrFctn* CorrFctn = *CorrFctnIter;
1130  CorrFctn->AddRealTriplet(TheTriplet);
1131  }
1132  } // if passed Triplet cut
1133  } // loop over third particle
1134  } // loop over second particle
1135  } // loop over first particle
1136 
1137  }
1138 
1139  delete TheTriplet;
1140 
1141  return NumTriplets;
1142 
1143 }
1144 
1145 
1146 //_________________________
1147 int StHbtThreeParticleAnalysis::CreateRealTriplets(StHbtSectoredPicoEvent *picoEvent, int Index1, int Index2, int Index3) {
1148 
1149  int NumTriplets=0;
1150 
1151  if (Index1==-1 || Index2==-1 || Index3==-1) return 0;
1152 
1153  StHbtParticleCollection *partColl1 = picoEvent->FirstSectoredCollection()[Index1];
1154  StHbtParticleCollection *partColl2 = picoEvent->FirstSectoredCollection()[Index2];
1155  StHbtParticleCollection *partColl3 = picoEvent->FirstSectoredCollection()[Index3];
1156 
1157  if (!partColl1->size() || !partColl2->size() || !partColl2->size()) return 0;
1158 
1159  StHbtTriplet* TheTriplet = new StHbtTriplet;
1160 
1161  StHbtParticleIterator PartIter1;
1162  StHbtParticleIterator PartIter2;
1163  StHbtParticleIterator PartIter3;
1164  StHbtCorrFctnIterator CorrFctnIter;
1165  StHbtParticleIterator StartOuterLoop = partColl1->begin();
1166  StHbtParticleIterator EndOuterLoop = partColl1->end();
1167  StHbtParticleIterator StartMiddleLoop = partColl2->begin();
1168  StHbtParticleIterator EndMiddleLoop = partColl2->end();
1169  StHbtParticleIterator StartInnerLoop = partColl3->begin();
1170  StHbtParticleIterator EndInnerLoop = partColl3->end();
1171 
1172  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
1173  TheTriplet->SetTrack1(*PartIter1);
1174  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
1175  TheTriplet->SetTrack2(*PartIter2);
1176  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
1177  TheTriplet->SetTrack3(*PartIter3);
1178  // The following lines have to be uncommented if you want TripletCutMonitors
1179  // they are not in for speed reasons
1180  // bool tmpPassTriplet = mTripletCut->Pass(TheTriplet);
1181  // mTripletCut->FillCutMonitor(TheTriplet, tmpPassTriplet);
1182  // if ( tmpPassTriplet ) {
1183  if (mTripletCut->Pass(TheTriplet)){
1184  NumTriplets++;
1185  for (CorrFctnIter=mCorrFctnCollection->begin();
1186  CorrFctnIter!=mCorrFctnCollection->end();CorrFctnIter++){
1187  StHbtCorrFctn* CorrFctn = *CorrFctnIter;
1188  CorrFctn->AddRealTriplet(TheTriplet);
1189  }
1190  } // if passed Triplet cut
1191  } // loop over third particle
1192  } // loop over second particle
1193  } // loop over first particle
1194 
1195  delete TheTriplet;
1196 
1197  return NumTriplets;
1198 
1199 }
1200 
1201 
1202 //_________________________
1203 int StHbtThreeParticleAnalysis::CalculateCosPhi(StHbtSectoredPicoEvent *picoEvent, int Index1) {
1204 
1205  double C2Q12,C2Q23,C2Q31,C3,arg1,arg2,arg3,arg4,cosphi,cosphiError,termt;
1206  int bin1,bin2,bin3,bin4,NumTriplets=0;
1207 
1208  StHbtParticleCollection *partColl1 = picoEvent->FirstSectoredCollection()[Index1];
1209 
1210  if (partColl1->size()<3) return 0;
1211 
1212  StHbtTriplet* TheTriplet = new StHbtTriplet;
1213 
1214  StHbtParticleIterator PartIter1;
1215  StHbtParticleIterator PartIter2;
1216  StHbtParticleIterator PartIter3;
1217  StHbtParticleIterator StartOuterLoop = partColl1->begin(); // always
1218  StHbtParticleIterator EndOuterLoop = partColl1->end(); // will be two less
1219  StHbtParticleIterator StartMiddleLoop;
1220  StHbtParticleIterator EndMiddleLoop;
1221  StHbtParticleIterator StartInnerLoop;
1222  StHbtParticleIterator EndInnerLoop;
1223 
1224  EndOuterLoop--;
1225  EndOuterLoop--; // outer loop goes to next-to-next-to-last particle in collection
1226  EndMiddleLoop = partColl1->end() ; // middle loop goes to next-to-last particle in collection
1227  EndMiddleLoop--;
1228  EndInnerLoop = partColl1->end() ; // inner loop goes to last particle in collection
1229 
1230  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
1231  StartMiddleLoop = PartIter1;
1232  StartMiddleLoop++;
1233  TheTriplet->SetTrack1(*PartIter1);
1234  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
1235  StartInnerLoop = PartIter2;
1236  StartInnerLoop++;
1237  TheTriplet->SetTrack2(*PartIter2);
1238  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
1239  TheTriplet->SetTrack3(*PartIter3);
1240  // The following lines have to be uncommented if you want TripletCutMonitors
1241  // they are not in for speed reasons
1242  // bool tmpPassTriplet = mTripletCut->Pass(TheTriplet);
1243  // mTripletCut->FillCutMonitor(TheTriplet, tmpPassTriplet);
1244  // if ( tmpPassTriplet ) {
1245  if (mTripletCut->Pass(TheTriplet)){
1246  NumTriplets++;
1247  bin1 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv12());
1248  bin2 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv23());
1249  bin3 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv31());
1250  bin4 = mQ3CF->GetXaxis()->FindBin(TheTriplet->qInv());
1251 
1252  C2Q12 = mQ2CF->GetBinContent(bin1);
1253  C2Q23 = mQ2CF->GetBinContent(bin2);
1254  C2Q31 = mQ2CF->GetBinContent(bin3);
1255  C3 = mNormFactor*(mQ3CF->GetBinContent(bin4));
1256 
1257  termt = ::sqrt(fabs((C2Q12-1.0)*(C2Q23-1.0)*(C2Q31-1.0)));
1258  if (termt) {
1259  cosphi = ((C3-1.0) - (C2Q12-1.0) - (C2Q23-1.0) - (C2Q31-1.0))/(2.0*termt);
1260  arg1 = mQ3CF->GetBinError(bin4);
1261  arg2 = mQ2CF->GetBinError(bin1)*(1.0 + (C2Q23-1.0)*(C2Q31-1.0)*cosphi/termt);
1262  arg3 = mQ2CF->GetBinError(bin2)*(1.0 + (C2Q31-1.0)*(C2Q12-1.0)*cosphi/termt);
1263  arg4 = mQ2CF->GetBinError(bin3)*(1.0 + (C2Q12-1.0)*(C2Q23-1.0)*cosphi/termt);
1264  cosphiError = (1.0/(2.0*termt))*::sqrt(arg1*arg1+arg2*arg2+arg3*arg3+arg4*arg4);
1265  }
1266  else {
1267  cosphi = 0.0;
1268  cosphiError = 0.0;
1269  }
1270 
1271  mCosPhi->AddBinContent(bin4, cosphi);
1272  mCosPhiE->AddBinContent(bin4, cosphiError);
1273  mCosPhiN->AddBinContent(bin4);
1274 
1275  }// if passed Triplet cut
1276  } // loop over third particle
1277  } // loop over second particle
1278  } // loop over first particle
1279 
1280  delete TheTriplet;
1281 
1282  return NumTriplets;
1283 
1284 }
1285 
1286 
1287 
1288 //_________________________
1289 int StHbtThreeParticleAnalysis::CalculateCosPhi(StHbtSectoredPicoEvent *picoEvent, int Index1, int Index2) {
1290 
1291  double C2Q12,C2Q23,C2Q31,C3,arg1,arg2,arg3,arg4,cosphi,cosphiError,termt;
1292  int bin1,bin2,bin3,bin4,NumTriplets=0;
1293 
1294  if (Index1==-1 || Index2==-1) return 0;
1295 
1296  StHbtParticleCollection *partColl1 = picoEvent->FirstSectoredCollection()[Index1];
1297  StHbtParticleCollection *partColl2 = picoEvent->FirstSectoredCollection()[Index2];
1298 
1299  StHbtTriplet* TheTriplet = new StHbtTriplet;
1300 
1301  // First, take two particles from the first collection and one particle from the second.
1302 
1303  if (partColl1->size()>=2 && partColl2->size()>=1) {
1304 
1305  StHbtParticleIterator PartIter1;
1306  StHbtParticleIterator PartIter2;
1307  StHbtParticleIterator PartIter3;
1308  StHbtParticleIterator StartOuterLoop = partColl1->begin(); // always
1309  StHbtParticleIterator EndOuterLoop = partColl1->end(); // will be one less
1310  StHbtParticleIterator StartMiddleLoop; // depends on start value of outer loop
1311  StHbtParticleIterator EndMiddleLoop = partColl1->end();
1312  StHbtParticleIterator StartInnerLoop = partColl2->begin();
1313  StHbtParticleIterator EndInnerLoop = partColl2->end();
1314 
1315  EndOuterLoop--; // outer loop goes to next-to-last particle in First collection
1316 
1317  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
1318  StartMiddleLoop = PartIter1;
1319  StartMiddleLoop++;
1320  TheTriplet->SetTrack1(*PartIter1);
1321  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
1322  TheTriplet->SetTrack2(*PartIter2);
1323  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
1324  TheTriplet->SetTrack3(*PartIter3);
1325  // The following lines have to be uncommented if you want TripletCutMonitors
1326  // they are not in for speed reasons
1327  // bool tmpPassTriplet = mTripletCut->Pass(TheTriplet);
1328  // mTripletCut->FillCutMonitor(TheTriplet, tmpPassTriplet);
1329  // if ( tmpPassTriplet ) {
1330  if (mTripletCut->Pass(TheTriplet)){
1331  NumTriplets++;
1332  bin1 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv12());
1333  bin2 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv23());
1334  bin3 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv31());
1335  bin4 = mQ3CF->GetXaxis()->FindBin(TheTriplet->qInv());
1336 
1337  C2Q12 = mQ2CF->GetBinContent(bin1);
1338  C2Q23 = mQ2CF->GetBinContent(bin2);
1339  C2Q31 = mQ2CF->GetBinContent(bin3);
1340  C3 = mNormFactor*(mQ3CF->GetBinContent(bin4));
1341 
1342  termt = ::sqrt(fabs((C2Q12-1.0)*(C2Q23-1.0)*(C2Q31-1.0)));
1343  if (termt) {
1344  cosphi = ((C3-1.0) - (C2Q12-1.0) - (C2Q23-1.0) - (C2Q31-1.0))/(2.0*termt);
1345  arg1 = mQ3CF->GetBinError(bin4);
1346  arg2 = mQ2CF->GetBinError(bin1)*(1.0 + (C2Q23-1.0)*(C2Q31-1.0)*cosphi/termt);
1347  arg3 = mQ2CF->GetBinError(bin2)*(1.0 + (C2Q31-1.0)*(C2Q12-1.0)*cosphi/termt);
1348  arg4 = mQ2CF->GetBinError(bin3)*(1.0 + (C2Q12-1.0)*(C2Q23-1.0)*cosphi/termt);
1349  cosphiError = (1.0/(2.0*termt))*::sqrt(arg1*arg1+arg2*arg2+arg3*arg3+arg4*arg4);
1350  }
1351  else {
1352  cosphi = 0.0;
1353  cosphiError = 0.0;
1354  }
1355 
1356  mCosPhi->AddBinContent(bin4, cosphi);
1357  mCosPhiE->AddBinContent(bin4, cosphiError);
1358  mCosPhiN->AddBinContent(bin4);
1359 
1360  }// if passed Triplet cut
1361  } // loop over third particle
1362  } // loop over second particle
1363  } // loop over first particle
1364 
1365  }
1366 
1367 
1368  // Then, take two particles from the second collection and one particle from the first.
1369 
1370  if (partColl1->size()>=1 && partColl2->size()>=2) {
1371 
1372  StHbtParticleIterator PartIter1;
1373  StHbtParticleIterator PartIter2;
1374  StHbtParticleIterator PartIter3;
1375  StHbtParticleIterator StartOuterLoop = partColl2->begin(); // always
1376  StHbtParticleIterator EndOuterLoop = partColl2->end(); // will be one less
1377  StHbtParticleIterator StartMiddleLoop; // depends on start value of outer loop
1378  StHbtParticleIterator EndMiddleLoop = partColl2->end();
1379  StHbtParticleIterator StartInnerLoop = partColl1->begin();
1380  StHbtParticleIterator EndInnerLoop = partColl1->end();
1381 
1382  EndOuterLoop--; // outer loop goes to next-to-last particle in First collection
1383 
1384  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
1385  StartMiddleLoop = PartIter1;
1386  StartMiddleLoop++;
1387  TheTriplet->SetTrack1(*PartIter1);
1388  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
1389  TheTriplet->SetTrack2(*PartIter2);
1390  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
1391  TheTriplet->SetTrack3(*PartIter3);
1392  // The following lines have to be uncommented if you want TripletCutMonitors
1393  // they are not in for speed reasons
1394  // bool tmpPassTriplet = mTripletCut->Pass(TheTriplet);
1395  // mTripletCut->FillCutMonitor(TheTriplet, tmpPassTriplet);
1396  // if ( tmpPassTriplet ) {
1397  if (mTripletCut->Pass(TheTriplet)){
1398  NumTriplets++;
1399  bin1 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv12());
1400  bin2 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv23());
1401  bin3 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv31());
1402  bin4 = mQ3CF->GetXaxis()->FindBin(TheTriplet->qInv());
1403 
1404  C2Q12 = mQ2CF->GetBinContent(bin1);
1405  C2Q23 = mQ2CF->GetBinContent(bin2);
1406  C2Q31 = mQ2CF->GetBinContent(bin3);
1407  C3 = mNormFactor*(mQ3CF->GetBinContent(bin4));
1408 
1409  termt = ::sqrt(fabs((C2Q12-1.0)*(C2Q23-1.0)*(C2Q31-1.0)));
1410  if (termt) {
1411  cosphi = ((C3-1.0) - (C2Q12-1.0) - (C2Q23-1.0) - (C2Q31-1.0))/(2.0*termt);
1412  arg1 = mQ3CF->GetBinError(bin4);
1413  arg2 = mQ2CF->GetBinError(bin1)*(1.0 + (C2Q23-1.0)*(C2Q31-1.0)*cosphi/termt);
1414  arg3 = mQ2CF->GetBinError(bin2)*(1.0 + (C2Q31-1.0)*(C2Q12-1.0)*cosphi/termt);
1415  arg4 = mQ2CF->GetBinError(bin3)*(1.0 + (C2Q12-1.0)*(C2Q23-1.0)*cosphi/termt);
1416  cosphiError = (1.0/(2.0*termt))*::sqrt(arg1*arg1+arg2*arg2+arg3*arg3+arg4*arg4);
1417  }
1418  else {
1419  cosphi = 0.0;
1420  cosphiError = 0.0;
1421  }
1422 
1423  mCosPhi->AddBinContent(bin4, cosphi);
1424  mCosPhiE->AddBinContent(bin4, cosphiError);
1425  mCosPhiN->AddBinContent(bin4);
1426 
1427  } // if passed Triplet cut
1428  } // loop over third particle
1429  } // loop over second particle
1430  } // loop over first particle
1431 
1432  }
1433 
1434  delete TheTriplet;
1435 
1436  return NumTriplets;
1437 
1438 }
1439 
1440 
1441 //_________________________
1442 int StHbtThreeParticleAnalysis::CalculateCosPhi(StHbtSectoredPicoEvent *picoEvent, int Index1, int Index2, int Index3) {
1443 
1444  double C2Q12,C2Q23,C2Q31,C3,arg1,arg2,arg3,arg4,cosphi,cosphiError,termt;
1445  int bin1,bin2,bin3,bin4,NumTriplets=0;
1446 
1447  if (Index1==-1 || Index2==-1 || Index3==-1) return 0;
1448 
1449  StHbtParticleCollection *partColl1 = picoEvent->FirstSectoredCollection()[Index1];
1450  StHbtParticleCollection *partColl2 = picoEvent->FirstSectoredCollection()[Index2];
1451  StHbtParticleCollection *partColl3 = picoEvent->FirstSectoredCollection()[Index3];
1452 
1453  if (!partColl1->size() || !partColl2->size() || !partColl2->size()) return 0;
1454 
1455  StHbtTriplet* TheTriplet = new StHbtTriplet;
1456 
1457  StHbtParticleIterator PartIter1;
1458  StHbtParticleIterator PartIter2;
1459  StHbtParticleIterator PartIter3;
1460  StHbtParticleIterator StartOuterLoop = partColl1->begin();
1461  StHbtParticleIterator EndOuterLoop = partColl1->end();
1462  StHbtParticleIterator StartMiddleLoop = partColl2->begin();
1463  StHbtParticleIterator EndMiddleLoop = partColl2->end();
1464  StHbtParticleIterator StartInnerLoop = partColl3->begin();
1465  StHbtParticleIterator EndInnerLoop = partColl3->end();
1466 
1467  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
1468  TheTriplet->SetTrack1(*PartIter1);
1469  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
1470  TheTriplet->SetTrack2(*PartIter2);
1471  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
1472  TheTriplet->SetTrack3(*PartIter3);
1473  // The following lines have to be uncommented if you want TripletCutMonitors
1474  // they are not in for speed reasons
1475  // bool tmpPassTriplet = mTripletCut->Pass(TheTriplet);
1476  // mTripletCut->FillCutMonitor(TheTriplet, tmpPassTriplet);
1477  // if ( tmpPassTriplet ) {
1478  if (mTripletCut->Pass(TheTriplet)){
1479  NumTriplets++;
1480  bin1 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv12());
1481  bin2 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv23());
1482  bin3 = mQ2CF->GetXaxis()->FindBin(TheTriplet->qInv31());
1483  bin4 = mQ3CF->GetXaxis()->FindBin(TheTriplet->qInv());
1484 
1485  C2Q12 = mQ2CF->GetBinContent(bin1);
1486  C2Q23 = mQ2CF->GetBinContent(bin2);
1487  C2Q31 = mQ2CF->GetBinContent(bin3);
1488  C3 = mNormFactor*(mQ3CF->GetBinContent(bin4));
1489 
1490  termt = ::sqrt(fabs((C2Q12-1.0)*(C2Q23-1.0)*(C2Q31-1.0)));
1491  if (termt) {
1492  cosphi = ((C3-1.0) - (C2Q12-1.0) - (C2Q23-1.0) - (C2Q31-1.0))/(2.0*termt);
1493  arg1 = mQ3CF->GetBinError(bin4);
1494  arg2 = mQ2CF->GetBinError(bin1)*(1.0 + (C2Q23-1.0)*(C2Q31-1.0)*cosphi/termt);
1495  arg3 = mQ2CF->GetBinError(bin2)*(1.0 + (C2Q31-1.0)*(C2Q12-1.0)*cosphi/termt);
1496  arg4 = mQ2CF->GetBinError(bin3)*(1.0 + (C2Q12-1.0)*(C2Q23-1.0)*cosphi/termt);
1497  cosphiError = (1.0/(2.0*termt))*::sqrt(arg1*arg1+arg2*arg2+arg3*arg3+arg4*arg4);
1498  }
1499  else {
1500  cosphi = 0.0;
1501  cosphiError = 0.0;
1502  }
1503 
1504  mCosPhi->AddBinContent(bin4, cosphi);
1505  mCosPhiE->AddBinContent(bin4, cosphiError);
1506  mCosPhiN->AddBinContent(bin4);
1507 
1508  } // if passed Triplet cut
1509  } // loop over third particle
1510  } // loop over second particle
1511  } // loop over first particle
1512 
1513  delete TheTriplet;
1514 
1515  return NumTriplets;
1516 
1517 }
1518 
1519 //_________________________
1520 int StHbtThreeParticleAnalysis::CreateMixedTriplets(StHbtParticleCollection *partColl1, StHbtParticleCollection *partColl2, StHbtParticleCollection *partColl3) {
1521 
1522  int NumTriplets=0;
1523 
1524  if (!partColl1->size() || !partColl1->size() || !partColl1->size()) return 0;
1525 
1526  StHbtTriplet* TheTriplet = new StHbtTriplet;
1527 
1528  StHbtParticleIterator PartIter1;
1529  StHbtParticleIterator PartIter2;
1530  StHbtParticleIterator PartIter3;
1531  StHbtCorrFctnIterator CorrFctnIter;
1532  StHbtParticleIterator StartOuterLoop = partColl1->begin();
1533  StHbtParticleIterator EndOuterLoop = partColl1->end();
1534  StHbtParticleIterator StartMiddleLoop = partColl2->begin();
1535  StHbtParticleIterator EndMiddleLoop = partColl2->end();
1536  StHbtParticleIterator StartInnerLoop = partColl3->begin();
1537  StHbtParticleIterator EndInnerLoop = partColl3->end();
1538 
1539  for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
1540  TheTriplet->SetTrack1(*PartIter1);
1541  for (PartIter2 = StartMiddleLoop; PartIter2!=EndMiddleLoop;PartIter2++){
1542  TheTriplet->SetTrack2(*PartIter2);
1543  for (PartIter3 = StartInnerLoop; PartIter3!=EndInnerLoop;PartIter3++){
1544  TheTriplet->SetTrack3(*PartIter3);
1545  // The following lines have to be uncommented if you want TripletCutMonitors
1546  // they are not in for speed reasons
1547  // bool tmpPassTriplet = mTripletCut->Pass(TheTriplet);
1548  // mTripletCut->FillCutMonitor(TheTriplet, tmpPassTriplet);
1549  // if ( tmpPassTriplet ) {
1550  if (mTripletCut->Pass(TheTriplet)){
1551  NumTriplets++;
1552  for (CorrFctnIter=mCorrFctnCollection->begin();
1553  CorrFctnIter!=mCorrFctnCollection->end();CorrFctnIter++){
1554  StHbtCorrFctn* CorrFctn = *CorrFctnIter;
1555  CorrFctn->AddMixedTriplet(TheTriplet);
1556  }
1557  } // if passed Triplet cut
1558  } // loop over third particle
1559  } // loop over second particle
1560  } // loop over first particle
1561 
1562  delete TheTriplet;
1563 
1564  return NumTriplets;
1565 
1566 }
1567 
1568 //_________________________
1569 void StHbtThreeParticleAnalysis::EventBegin(const StHbtEvent* ev){
1570  mFirstParticleCut->EventBegin(ev);
1571  mSecondParticleCut->EventBegin(ev);
1572  mThirdParticleCut->EventBegin(ev);
1573  for (StHbtCorrFctnIterator iter=mCorrFctnCollection->begin(); iter!=mCorrFctnCollection->end();iter++){
1574  (*iter)->EventBegin(ev);
1575  }
1576 }
1577 //_________________________
1578 void StHbtThreeParticleAnalysis::EventEnd(const StHbtEvent* ev){
1579  mFirstParticleCut->EventBegin(ev);
1580  mSecondParticleCut->EventBegin(ev);
1581  mThirdParticleCut->EventBegin(ev);
1582  for (StHbtCorrFctnIterator iter=mCorrFctnCollection->begin(); iter!=mCorrFctnCollection->end();iter++){
1583  (*iter)->EventEnd(ev);
1584  }
1585 }
1586 //_________________________
1587 void StHbtThreeParticleAnalysis::Finish(){
1588  if (mCalcCosPhi) {
1589 
1590  mCosPhiE->Divide(mCosPhiN);
1591  mCosPhi->Divide(mCosPhiN);
1592  for (int Iter1=1; Iter1<=mCosPhiN->GetNbinsX();Iter1++)
1593  mCosPhi->SetBinError(Iter1, mCosPhiE->GetBinContent(Iter1));
1594 
1595  cout << "Writing to file: CosPhi" << endl;
1596  TFile cosfile(mSaveFile, "NEW");
1597  mCosPhi->Write("cosphi");
1598  mCosPhiN->Write("cosphiN");
1599  cosfile.Close();
1600  cout << "Done writing: CosPhi" << endl;
1601  }
1602  StHbtCorrFctnIterator iter;
1603  for (iter=mCorrFctnCollection->begin(); iter!=mCorrFctnCollection->end();iter++){
1604  (*iter)->Finish();
1605  }
1606 }
1607 
1608 //_________________________
1609 void StHbtThreeParticleAnalysis::AddEventProcessed() {
1610  mNeventsProcessed++;
1611 }
1612 
StHbtThreeParticleAnalysis::ProcessEvent
virtual void ProcessEvent(const StHbtEvent *)
returns reports of all cuts applied and correlation functions being done
Definition: StHbtThreeParticleAnalysis.cxx:266
Generated by   doxygen 1.8.5