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