Initial check-in of the model classes
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / Analysis / AliFemtoLikeSignAnalysis.cxx
1 /***************************************************************************
2  *
3  * $Id$
4  *
5  * Author: Frank Laue, Ohio State, Laue@mps.ohio-state.edu
6  ***************************************************************************
7  *
8  * Description: part of STAR HBT Framework: AliFemtoMaker package
9  *      This is the Class for Analysis objects.  Each of the simultaneous
10  *      Analyses running should have one of these instantiated.  They link
11  *      into the Manager in an Analysis Collection.
12  *
13  ***************************************************************************/
14
15 #include "Analysis/AliFemtoLikeSignAnalysis.h"
16 #include "Infrastructure/AliFemtoParticleCollection.h"
17 #include "Base/AliFemtoTrackCut.h"
18 #include "Base/AliFemtoV0Cut.h"
19 #include "Infrastructure/AliFemtoPicoEventCollectionVector.h"
20 #include "Infrastructure/AliFemtoPicoEventCollectionVectorHideAway.h"
21
22 #ifdef __ROOT__ 
23 ClassImp(AliFemtoLikeSignAnalysis)
24 #endif
25
26 // this little function used to apply ParticleCuts (TrackCuts or V0Cuts) and fill ParticleCollections of picoEvent
27 //  it is called from AliFemtoAnalysis::ProcessEvent()
28
29
30 extern void FillHbtParticleCollection(AliFemtoParticleCut*         partCut,
31                                      AliFemtoEvent*               hbtEvent,
32                                      AliFemtoParticleCollection*  partCollection);
33
34  
35 //____________________________
36 AliFemtoLikeSignAnalysis::AliFemtoLikeSignAnalysis(unsigned int bins, double min, double max) : 
37   AliFemtoAnalysis(),
38   fVertexBins(0),
39   fOverFlow(0),  
40   fUnderFlow(0)  
41 {
42   fVertexBins = bins;
43   fVertexZ[0] = min;
44   fVertexZ[1] = max;
45   fUnderFlow = 0; 
46   fOverFlow = 0; 
47   if (fMixingBuffer) delete fMixingBuffer;
48   fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexBins,fVertexZ[0],fVertexZ[1]);
49     /* no-op */
50 }
51 //____________________________
52 AliFemtoLikeSignAnalysis::AliFemtoLikeSignAnalysis(const AliFemtoLikeSignAnalysis& a) : 
53   AliFemtoAnalysis(a) ,
54   fVertexBins(0),
55   fOverFlow(0),  
56   fUnderFlow(0)  
57 {
58   fVertexBins = a.fVertexBins; 
59   fVertexZ[0] = a.fVertexZ[0]; 
60   fVertexZ[1] = a.fVertexZ[1];
61   fUnderFlow = 0; 
62   fOverFlow = 0; 
63   if (fMixingBuffer) delete fMixingBuffer;
64   fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexBins,fVertexZ[0],fVertexZ[1]);
65  }
66 //____________________________ 
67 AliFemtoLikeSignAnalysis::~AliFemtoLikeSignAnalysis(){
68   delete fPicoEventCollectionVectorHideAway; fPicoEventCollectionVectorHideAway=0;
69 }
70 //____________________________
71 AliFemtoString AliFemtoLikeSignAnalysis::Report()
72 {  
73   char Ctemp[200];
74   cout << "AliFemtoLikeSignAnalysis - constructing Report..."<<endl;
75   AliFemtoString temp = "-----------\nHbt Analysis Report:\n";
76   sprintf(Ctemp,"Events are mixed in %d bins in the range %E cm to %E cm.\n",fVertexBins,fVertexZ[0],fVertexZ[1]);
77   temp += Ctemp;
78   sprintf(Ctemp,"Events underflowing: %d\n",fUnderFlow);
79   temp += Ctemp;
80   sprintf(Ctemp,"Events overflowing: %d\n",fOverFlow);
81   temp += Ctemp;
82   sprintf(Ctemp,"Now adding AliFemtoAnalysis(base) Report\n");
83   temp += Ctemp; 
84   temp += "Adding AliFemtoAnalysis(base) Report now:\n";
85   temp += AliFemtoAnalysis::Report();
86   temp += "-------------\n";
87   AliFemtoString returnThis=temp;
88   return returnThis;
89 }
90 //_________________________
91 void AliFemtoLikeSignAnalysis::ProcessEvent(const AliFemtoEvent* hbtEvent) {
92   // get right mixing buffer
93   double vertexZ = hbtEvent->PrimVertPos().z();
94   fMixingBuffer = fPicoEventCollectionVectorHideAway->PicoEventCollection(vertexZ); 
95   if (!fMixingBuffer) {
96     if ( vertexZ < fVertexZ[0] ) fUnderFlow++;
97     if ( vertexZ > fVertexZ[1] ) fOverFlow++;
98     return;
99   }
100
101   // startup for EbyE 
102   EventBegin(hbtEvent);  
103   // event cut and event cut monitor
104   bool tmpPassEvent = fEventCut->Pass(hbtEvent);
105   fEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
106   if (tmpPassEvent) {
107       fNeventsProcessed++;
108       cout << "AliFemtoLikeSignAnalysis::ProcessEvent() - " << hbtEvent->TrackCollection()->size();
109       cout << " #track=" << hbtEvent->TrackCollection()->size();
110       // OK, analysis likes the event-- build a pico event from it, using tracks the analysis likes...
111       AliFemtoPicoEvent* picoEvent = new AliFemtoPicoEvent;       // this is what we will make pairs from and put in Mixing Buffer
112       FillHbtParticleCollection(fFirstParticleCut,(AliFemtoEvent*)hbtEvent,picoEvent->FirstParticleCollection());
113       if ( !(AnalyzeIdenticalParticles()) )
114         FillHbtParticleCollection(fSecondParticleCut,(AliFemtoEvent*)hbtEvent,picoEvent->SecondParticleCollection());
115       cout <<"   #particles in First, Second Collections: " <<
116         picoEvent->FirstParticleCollection()->size() << " " <<
117         picoEvent->SecondParticleCollection()->size() << endl;
118       
119       if (picoEvent->SecondParticleCollection()->size()*picoEvent->FirstParticleCollection()->size()==0) {
120         delete picoEvent;
121         cout << "AliFemtoLikeSignAnalysis - picoEvent deleted due to empty collection " <<endl; 
122         return;
123       }
124       // OK, pico event is built
125       // make real pairs...
126       
127       // Fabrice points out that we do not need to keep creating/deleting pairs all the time
128       // We only ever need ONE pair, and we can just keep changing internal pointers
129       // this should help speed things up
130       AliFemtoPair* ThePair = new AliFemtoPair;
131       
132       AliFemtoParticleIterator PartIter1;
133       AliFemtoParticleIterator PartIter2;
134       AliFemtoCorrFctnIterator CorrFctnIter;
135       AliFemtoParticleIterator StartOuterLoop = picoEvent->FirstParticleCollection()->begin();  // always
136       AliFemtoParticleIterator EndOuterLoop   = picoEvent->FirstParticleCollection()->end();    // will be one less if identical
137       AliFemtoParticleIterator StartInnerLoop;
138       AliFemtoParticleIterator EndInnerLoop;
139       if (AnalyzeIdenticalParticles()) {             // only use First collection
140         EndOuterLoop--;                                               // outer loop goes to next-to-last particle in First collection
141         EndInnerLoop = picoEvent->FirstParticleCollection()->end() ;  // inner loop goes to last particle in First collection
142       }
143       else {                                                          // nonidentical - loop over First and Second collections
144         StartInnerLoop = picoEvent->SecondParticleCollection()->begin(); // inner loop starts at first particle in Second collection
145         EndInnerLoop   = picoEvent->SecondParticleCollection()->end() ;  // inner loop goes to last particle in Second collection
146       }
147       // real pairs
148       for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
149         if (AnalyzeIdenticalParticles()){
150           StartInnerLoop = PartIter1;
151           StartInnerLoop++;
152         }
153         ThePair->SetTrack1(*PartIter1);
154         for (PartIter2 = StartInnerLoop; PartIter2!=EndInnerLoop;PartIter2++){
155           ThePair->SetTrack2(*PartIter2);
156           // The following lines have to be uncommented if you want pairCutMonitors
157           // they are not in for speed reasons
158           // bool tmpPassPair = mPairCut->Pass(ThePair);
159           // mPairCut->FillCutMonitor(ThePair, tmpPassPair);
160           // if ( tmpPassPair ) {
161           if (fPairCut->Pass(ThePair)){
162             for (CorrFctnIter=fCorrFctnCollection->begin();
163                  CorrFctnIter!=fCorrFctnCollection->end();CorrFctnIter++){
164               AliFemtoLikeSignCorrFctn* CorrFctn = dynamic_cast<AliFemtoLikeSignCorrFctn*>(*CorrFctnIter);
165               if (CorrFctn) CorrFctn->AddRealPair(ThePair);
166             }
167           }  // if passed pair cut
168         }    // loop over second particle
169       }      // loop over first particle
170 #ifdef STHBTDEBUG
171       cout << "AliFemtoLikeSignAnalysis::ProcessEvent() - reals done" << endl;
172 #endif
173
174       AliFemtoParticleIterator nextIter;
175       AliFemtoParticleIterator prevIter;
176
177       // like sign first partilce collection pairs
178       prevIter = EndOuterLoop;
179       prevIter--;
180       for (PartIter1=StartOuterLoop;PartIter1!=prevIter;PartIter1++){
181         ThePair->SetTrack1(*PartIter1);
182         nextIter = PartIter1;
183         nextIter++;
184         for (PartIter2 = nextIter; PartIter2!=EndOuterLoop;PartIter2++){
185           ThePair->SetTrack2(*PartIter2);
186           // The following lines have to be uncommented if you want pairCutMonitors
187           // they are not in for speed reasons
188           // bool tmpPassPair = mPairCut->Pass(ThePair);
189           // mPairCut->FillCutMonitor(ThePair, tmpPassPair);
190           // if ( tmpPassPair ) {
191           if (fPairCut->Pass(ThePair)){
192             for (CorrFctnIter=fCorrFctnCollection->begin();
193                  CorrFctnIter!=fCorrFctnCollection->end();CorrFctnIter++){
194               AliFemtoLikeSignCorrFctn* CorrFctn = dynamic_cast<AliFemtoLikeSignCorrFctn*>(*CorrFctnIter);
195               if (CorrFctn) CorrFctn->AddLikeSignPositivePair(ThePair);
196             }
197           }  // if passed pair cut
198         }    // loop over second particle
199       }      // loop over first particle
200 #ifdef STHBTDEBUG
201       cout << "AliFemtoLikeSignAnalysis::ProcessEvent() - like sign first collection done" << endl;
202 #endif
203       // like sign second partilce collection pairs
204       prevIter = EndInnerLoop;
205       prevIter--;
206       for (PartIter1=StartInnerLoop;PartIter1!=prevIter;PartIter1++){
207         ThePair->SetTrack1(*PartIter1);
208         nextIter = PartIter1;
209         nextIter++;
210         for (PartIter2 = nextIter; PartIter2!=EndInnerLoop;PartIter2++){
211           ThePair->SetTrack2(*PartIter2);
212           // The following lines have to be uncommented if you want pairCutMonitors
213           // they are not in for speed reasons
214           // bool tmpPassPair = mPairCut->Pass(ThePair);
215           // mPairCut->FillCutMonitor(ThePair, tmpPassPair);
216           // if ( tmpPassPair ) {
217           if (fPairCut->Pass(ThePair)){
218             for (CorrFctnIter=fCorrFctnCollection->begin();
219                  CorrFctnIter!=fCorrFctnCollection->end();CorrFctnIter++){
220               AliFemtoLikeSignCorrFctn* CorrFctn = dynamic_cast<AliFemtoLikeSignCorrFctn*>(*CorrFctnIter);
221               if (CorrFctn) CorrFctn->AddLikeSignNegativePair(ThePair);
222             }
223           }  // if passed pair cut
224         }    // loop over second particle
225       }      // loop over first particle
226 #ifdef STHBTDEBUG
227       cout << "AliFemtoLikeSignAnalysis::ProcessEvent() - like sign second collection done" << endl;
228 #endif
229       
230       if (MixingBufferFull()){
231 #ifdef STHBTDEBUG
232         cout << "Mixing Buffer is full - lets rock and roll" << endl;
233 #endif
234       }
235       else {
236         cout << "Mixing Buffer not full -gotta wait " << MixingBuffer()->size() << endl;
237       }
238       if (MixingBufferFull()){
239         StartOuterLoop = picoEvent->FirstParticleCollection()->begin();
240         EndOuterLoop   = picoEvent->FirstParticleCollection()->end();
241         AliFemtoPicoEvent* storedEvent;
242         AliFemtoPicoEventIterator picoEventIter;
243         for (picoEventIter=MixingBuffer()->begin();picoEventIter!=MixingBuffer()->end();picoEventIter++){
244           storedEvent = *picoEventIter;
245           if (AnalyzeIdenticalParticles()){
246             StartInnerLoop = storedEvent->FirstParticleCollection()->begin();
247             EndInnerLoop = storedEvent->FirstParticleCollection()->end();
248           }
249           else{
250             StartInnerLoop = storedEvent->SecondParticleCollection()->begin();
251             EndInnerLoop = storedEvent->SecondParticleCollection()->end();
252           }
253           for (PartIter1=StartOuterLoop;PartIter1!=EndOuterLoop;PartIter1++){
254             ThePair->SetTrack1(*PartIter1);
255             for (PartIter2=StartInnerLoop;PartIter2!=EndInnerLoop;PartIter2++){
256               ThePair->SetTrack2(*PartIter2);
257               // testing...           cout << "ThePair defined... going to pair cut... ";
258               if (fPairCut->Pass(ThePair)){
259                 // testing...           cout << " ThePair passed PairCut... ";
260                 for (CorrFctnIter=fCorrFctnCollection->begin();
261                      CorrFctnIter!=fCorrFctnCollection->end();CorrFctnIter++){
262                   AliFemtoLikeSignCorrFctn* CorrFctn = dynamic_cast<AliFemtoLikeSignCorrFctn*>(*CorrFctnIter);
263                   if (CorrFctn) { 
264                     CorrFctn->AddMixedPair(ThePair);
265                     //cout << " ThePair has been added to MixedPair method " << endl;
266                   }
267                 }
268               }  // if passed pair cut
269             }    // loop over second particle
270           }      // loop over first particle
271         }        // loop over pico-events stored in Mixing buffer
272         // Now get rid of oldest stored pico-event in buffer.
273         // This means (1) delete the event from memory, (2) "pop" the pointer to it from the MixingBuffer
274         delete MixingBuffer()->back();
275         MixingBuffer()->pop_back();
276       }  // if mixing buffer is full
277       delete ThePair;
278       MixingBuffer()->push_front(picoEvent);  // store the current pico-event in buffer
279     }   // if currentEvent is accepted by currentAnalysis
280     EventEnd(hbtEvent);  // cleanup for EbyE 
281     //    cout << "AliFemtoLikeSignAnalysis::ProcessEvent() - return to caller ... " << endl;
282 }
283
284
285