StRoot: StHbtMaker/Infrastructure/StHbtSectoredAnalysis.cxx Source File

00001 /***************************************************************************
00002  *
00003  * $Id: StHbtSectoredAnalysis.cxx,v 1.3 2000/08/11 16:35:41 rcwells Exp $
00004  *
00005  * Author: Robert Willson, Ohio State, willson@bnl.gov
00006  ***************************************************************************
00007  *  
00008  * Description: part of STAR HBT Framework: StHbtMaker package.
00009  *      This is a derived class for two-particle analysis objects.  
00010  *      In this analysis, momentum space is sectored into a number of
00011  *      cubes, which are then mixed with nearest neighbors.  Several
00012  *      initializations are needed to set up the sectoring, see 
00013  *      StHbtSectoredExample.C in the /doc directory for use.
00014  *
00015  ***************************************************************************
00016  *
00017  * $Log: StHbtSectoredAnalysis.cxx,v $
00018  * Revision 1.3  2000/08/11 16:35:41  rcwells
00019  * Added number of events processed to each HBT analysis
00020  *
00021  * Revision 1.2  2000/07/16 21:38:23  laue
00022  * StHbtCoulomb.cxx StHbtSectoredAnalysis.cxx : updated for standalone version
00023  * StHbtV0.cc StHbtV0.hh : some cast to prevent compiling warnings
00024  * StHbtParticle.cc StHbtParticle.hh : pointers mTrack,mV0 initialized to 0
00025  * StHbtIOBinary.cc : some printouts in #ifdef STHBTDEBUG
00026  * StHbtEvent.cc : B-Field set to 0.25Tesla, we have to think about a better
00027  *                 solution
00028  *
00029  * Revision 1.1  2000/04/12 01:46:18  willson
00030  * Initial Installation
00031  *
00032  * 
00033  ***************************************************************************/
00034 
00035 #include "StHbtMaker/Infrastructure/StHbtSectoredAnalysis.h"
00036 #include "StHbtMaker/Infrastructure/StHbtParticleCollection.hh"
00037 #include "StHbtMaker/Base/StHbtTrackCut.h"
00038 #include "StHbtMaker/Base/StHbtV0Cut.h"
00039 
00040 #ifdef __ROOT__ 
00041 ClassImp(StHbtSectoredAnalysis)
00042 #endif
00043 
00044 //This function gives the correct linear transformation from a 3-D coordinate to a 1-D array index. 
00045 
00046 //____________________________
00047 int StHbtSectoredAnalysis::Index(int binx, int biny, int binz)
00048 {
00049   return binx + biny*mNumBinsX + binz*mNumBinsX*mNumBinsY;
00050 }
00051 
00052 
00053 
00054 //This function sorts the event into a number of sectors.  The overflow is the last array element.
00055 
00056 //____________________________
00057 void StHbtSectoredAnalysis::SortHbtParticleCollection(StHbtParticleCut*         partCut,
00058                                                       StHbtEvent*               hbtEvent,
00059                                                       StHbtParticleCollection** SectoredPicoEvent)
00060 {
00061   int i,j,k;
00062   switch (partCut->Type()) {
00063   case hbtTrack:       // cut is cutting on Tracks
00064     {
00065       StHbtTrackCut* pCut = (StHbtTrackCut*) partCut;
00066       StHbtTrack* pParticle;
00067       StHbtTrackIterator pIter;
00068       StHbtTrackIterator startLoop = hbtEvent->TrackCollection()->begin();
00069       StHbtTrackIterator endLoop   = hbtEvent->TrackCollection()->end();
00070 
00071       for (pIter=startLoop;pIter!=endLoop;pIter++){
00072         pParticle = *pIter;
00073         bool tmpPassParticle = pCut->Pass(pParticle);
00074         pCut->FillCutMonitor(pParticle, tmpPassParticle);
00075         if (tmpPassParticle){
00076           StHbtParticle* particle = new StHbtParticle(pParticle,pCut->Mass());
00077           if (particle->FourMomentum().px()<mPXmin || particle->FourMomentum().px()>mPXmax ||
00078               particle->FourMomentum().py()<mPYmin || particle->FourMomentum().py()>mPYmax ||
00079               particle->FourMomentum().pz()<mPZmin || particle->FourMomentum().pz()>mPZmax)
00080             SectoredPicoEvent[mNumBinsX*mNumBinsY*mNumBinsZ]->push_back(particle);
00081           else {
00082             i = (int) floor((particle->FourMomentum().px()-mPXmin)/mDeltaP);
00083             j = (int) floor((particle->FourMomentum().py()-mPYmin)/mDeltaP);
00084             k = (int) floor((particle->FourMomentum().pz()-mPZmin)/mDeltaP);
00085             SectoredPicoEvent[Index(i,j,k)]->push_back(particle);
00086           }
00087         }
00088       }
00089       break;      
00090     }
00091     
00092   case hbtV0:          // cut is cutting on V0s
00093     {
00094       StHbtV0Cut* pCut = (StHbtV0Cut*) partCut;
00095       StHbtV0* pParticle;
00096       StHbtV0Iterator pIter;
00097       StHbtV0Iterator startLoop = hbtEvent->V0Collection()->begin();
00098       StHbtV0Iterator endLoop   = hbtEvent->V0Collection()->end();
00099       // this following "for" loop is identical to the one above, but because of scoping, I can's see how to avoid repetition...
00100       
00101       for (pIter=startLoop;pIter!=endLoop;pIter++){
00102         pParticle = *pIter;
00103         bool tmpPassV0 = pCut->Pass(pParticle);
00104         pCut->FillCutMonitor(pParticle, tmpPassV0);
00105         if (tmpPassV0){
00106           StHbtParticle* particle = new StHbtParticle(pParticle,pCut->Mass());
00107           if (particle->FourMomentum().px()<mPXmin || particle->FourMomentum().px()>mPXmax ||
00108               particle->FourMomentum().py()<mPYmin || particle->FourMomentum().py()>mPYmax ||
00109               particle->FourMomentum().pz()<mPZmin || particle->FourMomentum().pz()>mPZmax)
00110             SectoredPicoEvent[mNumBinsX*mNumBinsY*mNumBinsZ]->push_back(particle);
00111           else {
00112             i = (int) floor((particle->FourMomentum().px()-mPXmin)/mDeltaP);
00113             j = (int) floor((particle->FourMomentum().py()-mPYmin)/mDeltaP);
00114             k = (int) floor((particle->FourMomentum().pz()-mPZmin)/mDeltaP);
00115             SectoredPicoEvent[Index(i,j,k)]->push_back(particle);
00116           }
00117         }
00118       }
00119       break;
00120     }
00121   default:
00122     cout << "FillHbtParticleCollection function (in StHbtSectoredAnalysis.cxx) - undefined Particle Cut type!!! \n";
00123   } 
00124 }
00125 
00126 
00127 //____________________________
00128 StHbtSectoredAnalysis::StHbtSectoredAnalysis(){
00129   //  mControlSwitch     = 0;
00130   mEventCut          = 0;
00131   mFirstParticleCut  = 0;
00132   mSecondParticleCut = 0;
00133   mPairCut           = 0;
00134   mCorrFctnCollection= 0;
00135   mNeventsProcessed = 0;
00136   // Defaults for sectoring
00137   mPXmax = 5.0;
00138   mPXmin = -5.0;
00139   mPYmax = 5.0;
00140   mPYmin = -5.0;
00141   mPZmax = 5.0;
00142   mPZmin = -5.0;
00143   mDeltaP = 10.0;
00144   mNumBinsX = 1;
00145   mNumBinsY = 1;
00146   mNumBinsZ = 1;
00147   mCorrFctnCollection = new StHbtCorrFctnCollection;
00148   mSectoredMixingBuffer = new StHbtSectoredPicoEventCollection;
00149 }
00150 
00151 //____________________________
00152 StHbtSectoredAnalysis::StHbtSectoredAnalysis(const StHbtSectoredAnalysis& a) : StHbtBaseAnalysis() {
00153   //  mControlSwitch     = 0;
00154   mEventCut          = 0;
00155   mFirstParticleCut  = 0;
00156   mSecondParticleCut = 0;
00157   mPairCut           = 0;
00158   mCorrFctnCollection= 0;
00159   mNeventsProcessed = 0;
00160   mPXmax = 5.0;
00161   mPXmin = -5.0;
00162   mPYmax = 5.0;
00163   mPYmin = -5.0;
00164   mPZmax = 5.0;
00165   mPZmin = -5.0;
00166   mDeltaP = 10.0;
00167   mNumBinsX = 1;
00168   mNumBinsY = 1;
00169   mNumBinsZ = 1;
00170   mCorrFctnCollection = new StHbtCorrFctnCollection;
00171   mSectoredMixingBuffer = new StHbtSectoredPicoEventCollection;
00172   
00173   // find the right event cut
00174   mEventCut = a.mEventCut->Clone();
00175   // find the right first particle cut
00176   mFirstParticleCut = a.mFirstParticleCut->Clone();
00177   // find the right second particle cut
00178   if (a.mFirstParticleCut==a.mSecondParticleCut) 
00179     SetSecondParticleCut(mFirstParticleCut); // identical particle hbt
00180   else
00181   mSecondParticleCut = a.mSecondParticleCut->Clone();
00182 
00183   mPairCut = a.mPairCut->Clone();
00184   
00185   if ( mEventCut ) {
00186       SetEventCut(mEventCut); // this will set the myAnalysis pointer inside the cut
00187       cout << " StHbtAnalysis::StHbtAnalysis(const StHbtAnalysis& a) - event cut set " << endl;
00188   }
00189   if ( mFirstParticleCut ) {
00190       SetFirstParticleCut(mFirstParticleCut); // this will set the myAnalysis pointer inside the cut
00191       cout << " StHbtAnalysis::StHbtAnalysis(const StHbtAnalysis& a) - first particle cut set " << endl;
00192   }
00193   if ( mSecondParticleCut ) {
00194       SetSecondParticleCut(mSecondParticleCut); // this will set the myAnalysis pointer inside the cut
00195       cout << " StHbtAnalysis::StHbtAnalysis(const StHbtAnalysis& a) - second particle cut set " << endl;
00196   }  if ( mPairCut ) {
00197       SetPairCut(mPairCut); // this will set the myAnalysis pointer inside the cut
00198       cout << " StHbtAnalysis::StHbtAnalysis(const StHbtAnalysis& a) - pair cut set " << endl;
00199   }
00200 
00201   StHbtCorrFctnIterator iter;
00202   for (iter=a.mCorrFctnCollection->begin(); iter!=a.mCorrFctnCollection->end();iter++){
00203     cout << " StHbtAnalysis::StHbtAnalysis(const StHbtAnalysis& a) - looking for correlation functions " << endl;
00204     StHbtCorrFctn* fctn = (*iter)->Clone();
00205     if (fctn) AddCorrFctn(fctn);
00206     else cout << " StHbtAnalysis::StHbtAnalysis(const StHbtAnalysis& a) - correlation function not found " << endl;
00207   }
00208 
00209   mNumEventsToMix = a.mNumEventsToMix;
00210 
00211   cout << " StHbtSectoredAnalysis::StHbtSectoredAnalysis(const StHbtSectoredAnalysis& a) - analysis copied " << endl;
00212   
00213 }
00214 
00215 //____________________________
00216 StHbtSectoredAnalysis::~StHbtSectoredAnalysis(){
00217   delete mEventCut          ;
00218   delete mFirstParticleCut  ;
00219   delete mSecondParticleCut ;
00220   delete mPairCut           ;
00221   // now delete every CorrFunction in the Collection, and then the Collection itself
00222   StHbtCorrFctnIterator iter;
00223   for (iter=mCorrFctnCollection->begin(); iter!=mCorrFctnCollection->end();iter++){
00224     delete *iter;
00225   }
00226   delete mCorrFctnCollection;
00227   // now delete every PicoEvent in the EventMixingBuffer and then the Buffer itself
00228   StHbtSectoredPicoEventIterator piter;
00229   for (piter=mSectoredMixingBuffer->begin();piter!=mSectoredMixingBuffer->end();piter++){
00230     delete *piter;
00231   }
00232   delete mSectoredMixingBuffer;
00233   
00234 }
00235 //______________________
00236 StHbtCorrFctn* StHbtSectoredAnalysis::CorrFctn(int n){  // return pointer to n-th correlation function
00237   if ( n<0 || n > (int)mCorrFctnCollection->size() )
00238     return NULL;
00239   StHbtCorrFctnIterator iter=mCorrFctnCollection->begin();
00240   for (int i=0; i<n ;i++){
00241     iter++;
00242   }
00243   return *iter;
00244 }
00245 //____________________________
00246 StHbtString StHbtSectoredAnalysis::Report()
00247 {
00248   char Sect[100];
00249   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);
00250   cout << "StHbtSectoredAnalysis - constructing Report..."<<endl;
00251   string temp = "-----------\nHbt Analysis Report:\n";
00252   temp += "\nEvent Cuts:\n";
00253   temp += mEventCut->Report();
00254   temp += "\nParticle Cuts - First Particle:\n";
00255   temp += mFirstParticleCut->Report();
00256   temp += "\nParticle Cuts - Second Particle:\n";
00257   temp += mSecondParticleCut->Report();
00258   temp += "\nPair Cuts:\n";
00259   temp += mPairCut->Report();
00260   temp += "\nCorrelation Functions:\n";
00261   StHbtCorrFctnIterator iter;
00262   if ( mCorrFctnCollection->size()==0 ) {
00263     cout << "StHbtSectoredAnalysis-Warning : no correlations functions in this analysis " << endl;
00264   }
00265   for (iter=mCorrFctnCollection->begin(); iter!=mCorrFctnCollection->end();iter++){
00266     temp += (*iter)->Report();
00267     temp += "\n";
00268   }
00269   temp += "\nSectoring Information:\n";
00270   temp += Sect;
00271   temp += "\n-------------\n";
00272   StHbtString returnThis=temp;
00273   return returnThis;
00274 }
00275 //_________________________
00276 void StHbtSectoredAnalysis::ProcessEvent(const StHbtEvent* hbtEvent) {
00277 
00278   // Add event to processed events
00279   AddEventProcessed();
00280 
00281   int i, j, k, i2, j2, k2;   // used in loops over sectors
00282   int size1=0, size2=0;      // To get the total size of all sectors
00283   bool DidOverflow;          // To only count the overflow bin once per sector
00284    
00285   // startup for EbyE 
00286   EventBegin(hbtEvent);  
00287   // event cut and event cut monitor
00288   bool tmpPassEvent = mEventCut->Pass(hbtEvent);
00289   mEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
00290   if (tmpPassEvent) {
00291     cout << "StHbtSectoredAnalysis::ProcessEvent() - Event has passed cut - build sectoredpicoEvent from " <<
00292       hbtEvent->TrackCollection()->size() << " tracks in TrackCollection" << endl;
00293     // OK, analysis likes the event-- build a sectored pico event from it, using tracks the analysis likes...
00294     
00295     // this is what we will make Pairs from and put in Mixing Buffer
00296     StHbtSectoredPicoEvent* picoEvent = new StHbtSectoredPicoEvent(mNumBinsX, mNumBinsY, mNumBinsZ);     
00297     
00298     SortHbtParticleCollection(mFirstParticleCut,(StHbtEvent*)hbtEvent,picoEvent->FirstSectoredCollection());
00299     if ( !(AnalyzeIdenticalParticles()) ) 
00300       SortHbtParticleCollection(mSecondParticleCut,(StHbtEvent*)hbtEvent,picoEvent->SecondSectoredCollection());
00301     
00302     //Obtain combined size of all sectors
00303     
00304     for (i=0; i<mNumBinsX; i++) 
00305       for (j=0; j<mNumBinsY; j++) 
00306         for (k=0; k<mNumBinsZ; k++) {
00307           size1 += picoEvent->FirstSectoredCollection()[Index(i,j,k)]->size();
00308           if ( !(AnalyzeIdenticalParticles()) ) 
00309             size2 += picoEvent->SecondSectoredCollection()[Index(i,j,k)]->size();
00310         }
00311     
00312     cout <<"StHbtSectoredAnalysis::ProcessEvent - #particles in First, Second Collections: " <<
00313       size1 << " " << size2 << endl;
00314     
00315     if (AnalyzeIdenticalParticles()) {
00316       // Real Pairs--Identical Particles
00317       for (i=0; i<mNumBinsX; i++) 
00318         for (j=0; j<mNumBinsY; j++)
00319           for (k=0; k<mNumBinsZ; k++) {
00320             DidOverflow=false;
00321             
00322             // A sector sweep of a single sectoredpicoevent - 9 terms with i+1, 3 term with i, j+1, 
00323             // and 1 term with i, j, k+1.  Plus another term for mixing a sector with itself.
00324 
00325             //The 9 terms
00326               for (j2=j-1; j2<=j+1; j2++)
00327                 for (k2=k-1; k2<=k+1; k2++) {
00328                   if (j2<mNumBinsY && k2<mNumBinsZ && j2!=-1 && k2!=-1 && i+1<mNumBinsX) {
00329                     CreateRealPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], picoEvent->FirstSectoredCollection()[Index(i+1,j2,k2)]);
00330                   }
00331                   else 
00332                     if (!DidOverflow) {
00333                       CreateRealPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], picoEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ]);
00334                       DidOverflow=true;
00335                     }
00336                 }
00337 
00338               //The 3 terms
00339               for (k2=k-1; k2<=k+1; k2++) {
00340                 if (k2<mNumBinsZ && k2!=-1 && j+1<mNumBinsY) {
00341                   CreateRealPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], picoEvent->FirstSectoredCollection()[Index(i,j+1,k2)]);
00342                 }
00343                 else 
00344                   if (!DidOverflow) {
00345                     CreateRealPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], picoEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ]);
00346                     DidOverflow=true;
00347                   }
00348               }
00349 
00350               //The 1 term
00351               if (k+1<mNumBinsZ) {
00352                 CreateRealPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], picoEvent->FirstSectoredCollection()[Index(i,j,k+1)]);
00353               }
00354               else 
00355                 if (!DidOverflow) {
00356                   CreateRealPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], picoEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ]);
00357                   DidOverflow=true;
00358                 }
00359                     
00360               //Mix sector with itself
00361               CreateRealPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)]);
00362           }
00363       //Mix overflow bin with itself
00364       CreateRealPairs(picoEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ]);
00365           
00366     }  
00367     
00368     else {
00369       // Real Pairs--Non-Identical Particles
00370       for (i=0; i<mNumBinsX; i++) 
00371         for (j=0; j<mNumBinsY; j++)
00372           for (k=0; k<mNumBinsZ; k++) {
00373             DidOverflow=false;
00374             for (i2=i-1; i2<=i+1; i2++)
00375               for (j2=j-1; j2<=j+1; j2++)
00376                 for (k2=k-1; k2<=k+1; k2++)
00377                   if (i2<mNumBinsX && j2<mNumBinsY && k2<mNumBinsZ && i2!=-1 && j2!=-1 && k2!=-1)
00378                     CreateRealPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], picoEvent->SecondSectoredCollection()[Index(i2,j2,k2)]);
00379                   else 
00380                     if (!DidOverflow) {
00381                       CreateRealPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], picoEvent->SecondSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ]);
00382                       DidOverflow=true;
00383                     }
00384           }
00385     }
00386     //Need to mix with the overflow bin of FirstSectoredCollection:
00387     for (i=0; i<mNumBinsX; i++) 
00388       for (j=0; j<mNumBinsY; j++)
00389         for (k=0; k<mNumBinsZ; k++)
00390           if (i==0 || j==0 || k==0 || i==mNumBinsX-1 || j==mNumBinsY-1 || k==mNumBinsZ-1)
00391             CreateRealPairs(picoEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ], picoEvent->SecondSectoredCollection()[Index(i,j,k)]);
00392 
00393     //Need to mix both overflow bins
00394     CreateRealPairs(picoEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ], picoEvent->SecondSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ]);
00395     
00396     
00397     cout << "StHbtSectoredAnalysis::ProcessEvent() - reals done, " << endl; 
00398     
00399     // ok, now make mixed Pairs, if the Mixing buffer is full
00400     
00401     if (SectoredMixingBufferFull()){
00402       cout << "Mixing Buffer is full - lets rock and roll" << endl;
00403     }
00404     else {
00405       cout << "Mixing Buffer not full -gotta wait " << SectoredMixingBuffer()->size() << endl;
00406     }
00407     if (SectoredMixingBufferFull()){
00408       StHbtSectoredPicoEvent* storedEvent;
00409       StHbtSectoredPicoEventIterator picoEventIter;
00410       for (picoEventIter=SectoredMixingBuffer()->begin();picoEventIter!=SectoredMixingBuffer()->end();picoEventIter++){
00411         storedEvent = *picoEventIter;
00412         if (AnalyzeIdenticalParticles()){
00413           // Mixed Pairs--Identical Particles
00414           for (i=0; i<mNumBinsX; i++) 
00415             for (j=0; j<mNumBinsY; j++)
00416               for (k=0; k<mNumBinsZ; k++) {
00417                 DidOverflow=false;
00418                 for (i2=i-1; i2<=i+1; i2++)
00419                   for (j2=j-1; j2<=j+1; j2++)
00420                     for (k2=k-1; k2<=k+1; k2++)
00421                       if (i2<mNumBinsX && j2<mNumBinsY && k2<mNumBinsZ && i2!=-1 && j2!=-1 && k2!=-1) {
00422                         CreateMixedPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], storedEvent->FirstSectoredCollection()[Index(i2,j2,k2)]);
00423                       }
00424                       else 
00425                         if (!DidOverflow) {
00426                           CreateMixedPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], storedEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ]);
00427                           DidOverflow=true;
00428                         }
00429               }
00430           
00431           //Need to mix with the overflow bin of FirstSectoredCollection:
00432           for (i=0; i<mNumBinsX; i++) 
00433             for (j=0; j<mNumBinsY; j++)
00434               for (k=0; k<mNumBinsZ; k++)
00435                 if (i==0 || j==0 || k==0 || i==mNumBinsX-1 || j==mNumBinsY-1 || k==mNumBinsZ-1){
00436                   CreateMixedPairs(picoEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ], storedEvent->FirstSectoredCollection()[Index(i,j,k)]);
00437                 }
00438 
00439           //Need to mix the overflow bin of picoEvent with the overflow bin of storedEvent
00440           CreateMixedPairs(picoEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ], storedEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ]);
00441           
00442         }
00443           
00444         else{
00445           // Mixed Pairs--Non-Identical Particles
00446           for (i=0; i<mNumBinsX; i++) 
00447             for (j=0; j<mNumBinsY; j++)
00448               for (k=0; k<mNumBinsZ; k++) {
00449                 DidOverflow=false;
00450                 for (i2=i-1; i2<=i+1; i2++)
00451                   for (j2=j-1; j2<=j+1; j2++)
00452                     for (k2=k-1; k2<=k+1; k2++)
00453                       if (i2<mNumBinsX && j2<mNumBinsY && k2<mNumBinsZ && i2!=-1 && j2!=-1 && k2!=-1)
00454                         CreateMixedPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], storedEvent->SecondSectoredCollection()[Index(i2,j2,k2)]);     
00455                       else 
00456                         if (!DidOverflow) {
00457                           CreateMixedPairs(picoEvent->FirstSectoredCollection()[Index(i,j,k)], storedEvent->SecondSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ]);
00458                           DidOverflow=true;
00459                         }
00460               }
00461           //Need to mix with the overflow bin of FirstSectoredCollection:
00462           for (i=0; i<mNumBinsX; i++) 
00463             for (j=0; j<mNumBinsY; j++)
00464               for (k=0; k<mNumBinsZ; k++)
00465                 if (i==0 || j==0 || k==0 || i==mNumBinsX-1 || j==mNumBinsY-1 || k==mNumBinsZ-1)
00466                   CreateMixedPairs(picoEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ], storedEvent->SecondSectoredCollection()[Index(i,j,k)]);
00467           
00468           //Need to mix the overflow bin of picoEvent with the overflow bin of storedEvent
00469           CreateMixedPairs(picoEvent->FirstSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ], storedEvent->SecondSectoredCollection()[mNumBinsX*mNumBinsY*mNumBinsZ]);
00470          
00471         }
00472       }
00473       
00474       // Now get rid of oldest stored pico-event in buffer.
00475       // This means (1) delete the event from memory, (2) "pop" the pointer to it from the MixingBuffer
00476       picoEventIter = SectoredMixingBuffer()->end();
00477       picoEventIter--;   // bug fixed malisa 27jul99 - end() is one BEYOND the end! (besides crashing on linux, this was a memory leak)
00478       delete *picoEventIter;
00479       
00480       SectoredMixingBuffer()->pop_back();
00481     }  // if mixing buffer is full
00482     SectoredMixingBuffer()->push_front(picoEvent);  // store the current pico-event in buffer
00483   }    // if currentEvent is accepted by currentAnalysis
00484   EventEnd(hbtEvent);  // cleanup for EbyE 
00485   cout << "StHbtSectoredAnalysis::ProcessEvent() - return to caller ...  " << endl; 
00486     
00487 }
00488 
00489       
00490 //_________________________
00491 void StHbtSectoredAnalysis::CreateRealPairs(StHbtParticleCollection *partColl1) {
00492 
00493   if (partColl1->size()<2) return;
00494 
00495   StHbtPair* ThePair = new StHbtPair;
00496   
00497   StHbtParticleIterator PartIter1;
00498   StHbtParticleIterator PartIter2;
00499   StHbtCorrFctnIterator CorrFctnIter;
00500   StHbtParticleIterator StartOuterLoop = partColl1->begin();  // always
00501   StHbtParticleIterator EndOuterLoop   = partColl1->end();    // will be one less
00502   StHbtParticleIterator StartInnerLoop;
00503   StHbtParticleIterator EndInnerLoop;
00504         
00505   EndOuterLoop--;                     // outer loop goes to next-to-last particle in collection
00506   EndInnerLoop = partColl1->end() ; 
00507 
00508   for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
00509     StartInnerLoop = PartIter1;
00510     StartInnerLoop++;
00511     ThePair->SetTrack1(*PartIter1);
00512     for (PartIter2 = StartInnerLoop; PartIter2!=EndInnerLoop;PartIter2++){
00513       ThePair->SetTrack2(*PartIter2);
00514       // The following lines have to be uncommented if you want PairCutMonitors
00515       // they are not in for speed reasons
00516       // bool tmpPassPair = mPairCut->Pass(ThePair);
00517       // mPairCut->FillCutMonitor(ThePair, tmpPassPair);
00518       // if ( tmpPassPair ) {
00519       if (mPairCut->Pass(ThePair)){
00520         for (CorrFctnIter=mCorrFctnCollection->begin();
00521              CorrFctnIter!=mCorrFctnCollection->end();CorrFctnIter++){
00522           StHbtCorrFctn* CorrFctn = *CorrFctnIter;
00523           CorrFctn->AddRealPair(ThePair);
00524         }
00525       }  // if passed Pair cut
00526     }    // loop over second particle
00527   }      // loop over first particle   
00528   
00529   delete ThePair;
00530  
00531 }
00532 
00533 //_________________________
00534 void StHbtSectoredAnalysis::CreateRealPairs(StHbtParticleCollection *partColl1, StHbtParticleCollection *partColl2) {
00535 
00536   if (partColl1->size()<1 || partColl2->size()<1) return;
00537 
00538   StHbtPair* ThePair = new StHbtPair;
00539         
00540   // We are taking one particle from each of the two collections
00541 
00542   StHbtParticleIterator PartIter1;
00543   StHbtParticleIterator PartIter2;
00544   StHbtCorrFctnIterator CorrFctnIter;
00545   StHbtParticleIterator StartOuterLoop = partColl1->begin(); 
00546   StHbtParticleIterator EndOuterLoop   = partColl1->end();    
00547   StHbtParticleIterator StartInnerLoop = partColl2->begin();  
00548   StHbtParticleIterator EndInnerLoop   = partColl2->end();
00549   
00550   for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
00551     ThePair->SetTrack1(*PartIter1);
00552     for (PartIter2 = StartInnerLoop; PartIter2!=EndInnerLoop;PartIter2++){
00553       ThePair->SetTrack2(*PartIter2);
00554       // The following lines have to be uncommented if you want PairCutMonitors
00555       // they are not in for speed reasons
00556       // bool tmpPassPair = mPairCut->Pass(ThePair);
00557       // mPairCut->FillCutMonitor(ThePair, tmpPassPair);
00558       // if ( tmpPassPair ) {
00559       if (mPairCut->Pass(ThePair)){
00560         for (CorrFctnIter=mCorrFctnCollection->begin();
00561              CorrFctnIter!=mCorrFctnCollection->end();CorrFctnIter++){
00562           StHbtCorrFctn* CorrFctn = *CorrFctnIter;
00563           CorrFctn->AddRealPair(ThePair);
00564         }
00565       }  // if passed Pair cut
00566     }    // loop over second particle
00567   }      // loop over first particle
00568   
00569   delete ThePair;
00570 
00571 }
00572 
00573 //_________________________
00574 void StHbtSectoredAnalysis::CreateMixedPairs(StHbtParticleCollection *partColl1, StHbtParticleCollection *partColl2) {
00575 
00576   if (partColl1->size()<1 || partColl2->size()<1) return;
00577 
00578   StHbtPair* ThePair = new StHbtPair;
00579         
00580   // We are taking one particle from each of the two collections
00581 
00582   StHbtParticleIterator PartIter1;
00583   StHbtParticleIterator PartIter2;
00584   StHbtCorrFctnIterator CorrFctnIter;
00585   StHbtParticleIterator StartOuterLoop = partColl1->begin(); 
00586   StHbtParticleIterator EndOuterLoop   = partColl1->end();    
00587   StHbtParticleIterator StartInnerLoop = partColl2->begin();  
00588   StHbtParticleIterator EndInnerLoop   = partColl2->end();
00589   
00590   for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
00591     ThePair->SetTrack1(*PartIter1);
00592     for (PartIter2 = StartInnerLoop; PartIter2!=EndInnerLoop;PartIter2++){
00593       ThePair->SetTrack2(*PartIter2);
00594       // The following lines have to be uncommented if you want PairCutMonitors
00595       // they are not in for speed reasons
00596       // bool tmpPassPair = mPairCut->Pass(ThePair);
00597       // mPairCut->FillCutMonitor(ThePair, tmpPassPair);
00598       // if ( tmpPassPair ) {
00599       if (mPairCut->Pass(ThePair)){
00600         for (CorrFctnIter=mCorrFctnCollection->begin();
00601              CorrFctnIter!=mCorrFctnCollection->end();CorrFctnIter++){
00602           StHbtCorrFctn* CorrFctn = *CorrFctnIter;
00603           CorrFctn->AddMixedPair(ThePair);
00604         }
00605       }  // if passed Pair cut
00606     }    // loop over second particle
00607   }      // loop over first particle
00608   
00609   delete ThePair;
00610 
00611 }
00612  
00613 //_________________________
00614 void StHbtSectoredAnalysis::EventBegin(const StHbtEvent* ev){
00615   mFirstParticleCut->EventBegin(ev);
00616   mSecondParticleCut->EventBegin(ev);
00617   mPairCut->EventBegin(ev);
00618   for (StHbtCorrFctnIterator iter=mCorrFctnCollection->begin(); iter!=mCorrFctnCollection->end();iter++){
00619     (*iter)->EventBegin(ev);
00620   }
00621 }
00622 //_________________________
00623 void StHbtSectoredAnalysis::EventEnd(const StHbtEvent* ev){
00624   mFirstParticleCut->EventBegin(ev);
00625   mSecondParticleCut->EventBegin(ev);
00626   mPairCut->EventEnd(ev);
00627   for (StHbtCorrFctnIterator iter=mCorrFctnCollection->begin(); iter!=mCorrFctnCollection->end();iter++){
00628     (*iter)->EventEnd(ev);
00629   }
00630 }
00631 //_________________________
00632 void StHbtSectoredAnalysis::Finish(){
00633   StHbtCorrFctnIterator iter;
00634   for (iter=mCorrFctnCollection->begin(); iter!=mCorrFctnCollection->end();iter++){
00635     (*iter)->Finish();
00636   }
00637 }
00638 //_________________________
00639 void StHbtSectoredAnalysis::AddEventProcessed() {
00640   mNeventsProcessed++;
00641 }
00642 
00643