Add Reaction Plane aware analysis
[u/mrichter/AliRoot.git] / PWG2 / FEMTOSCOPY / AliFemto / AliFemtoAnalysisReactionPlane.cxx
1 ////////////////////////////////////////////////////////////////////////////
2 //                                                                        //
3 // AliFemtoAnalysisReactionPlane - Femtoscopic analysis which mixes event //
4 // with respect to the z position of the primary vertex and event total   //
5 // multiplicity and uses only events in certain reaction plane angle bin  //
6 //                                                                        //
7 ////////////////////////////////////////////////////////////////////////////
8
9 #include <TMath.h>
10 #include "AliFemtoAnalysisReactionPlane.h"
11 #include "AliFemtoParticleCollection.h"
12 #include "AliFemtoTrackCut.h"
13 #include "AliFemtoV0Cut.h"
14 #include "AliFemtoKinkCut.h"
15 #include "AliFemtoPicoEventCollectionVector.h"
16 #include "AliFemtoPicoEventCollectionVectorHideAway.h"
17
18 #ifdef __ROOT__ 
19 ClassImp(AliFemtoAnalysisReactionPlane)
20 #endif
21
22 extern void FillHbtParticleCollection(AliFemtoParticleCut*         partCut,
23                                       AliFemtoEvent*               hbtEvent,
24                                       AliFemtoParticleCollection*  partCollection);
25
26
27 //____________________________
28 AliFemtoAnalysisReactionPlane::AliFemtoAnalysisReactionPlane(unsigned int binsVertex, double minVertex, double maxVertex,
29                                                        unsigned int binsMult, double minMult, double maxMult, unsigned short binsRP) 
30   : 
31   fVertexZBins(binsVertex), 
32   fOverFlowVertexZ(0), 
33   fUnderFlowVertexZ(0),
34   fMultBins(binsMult) ,
35   fOverFlowMult(0),    
36   fUnderFlowMult(0),
37   fRPBins(binsRP),
38   fCurrentRP(0)
39 {
40   //  mControlSwitch     = 0;
41   fEventCut          = 0;
42   fFirstParticleCut  = 0;
43   fSecondParticleCut = 0;
44   fPairCut           = 0;
45   fCorrFctnCollection= 0;
46   fCorrFctnCollection = new AliFemtoCorrFctnCollection;
47   fVertexZ[0] = minVertex;
48   fVertexZ[1] = maxVertex;
49   fUnderFlowVertexZ = 0; 
50   fOverFlowVertexZ = 0; 
51   fMult[0] = minMult;
52   fMult[1] = maxMult;
53   fUnderFlowMult = 0; 
54   fOverFlowMult = 0; 
55   if (fMixingBuffer) delete fMixingBuffer;
56   fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
57                                                                                      fMultBins,fMult[0],fMult[1],
58                                                                                      fRPBins,0.0,TMath::Pi());
59 }
60 //____________________________
61
62 AliFemtoAnalysisReactionPlane::AliFemtoAnalysisReactionPlane(const AliFemtoAnalysisReactionPlane& a) : 
63   AliFemtoSimpleAnalysis(),
64   fVertexZBins(a.fVertexZBins), 
65   fOverFlowVertexZ(0), 
66   fUnderFlowVertexZ(0),
67   fMultBins(a.fMultBins) ,
68   fOverFlowMult(0),    
69   fUnderFlowMult(0),
70   fRPBins(a.fRPBins),
71   fCurrentRP(0)
72 {
73   //AliFemtoAnalysisReactionPlane();
74   fEventCut          = 0;
75   fFirstParticleCut  = 0;
76   fSecondParticleCut = 0;
77   fPairCut           = 0;
78   fCorrFctnCollection= 0;
79   fCorrFctnCollection = new AliFemtoCorrFctnCollection;
80   fVertexZ[0] = a.fVertexZ[0]; 
81   fVertexZ[1] = a.fVertexZ[1];
82   fUnderFlowVertexZ = 0; 
83   fOverFlowVertexZ = 0; 
84   fMult[0] = a.fMult[0]; 
85   fMult[1] = a.fMult[1];
86   fUnderFlowMult = 0; 
87   fOverFlowMult = 0; 
88   if (fMixingBuffer) delete fMixingBuffer;
89   fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
90                                                                                      fMultBins,fMult[0],fMult[1],
91                                                                                      fRPBins,0.0,TMath::Pi());
92
93   // find the right event cut
94   fEventCut = a.fEventCut->Clone();
95   // find the right first particle cut
96   fFirstParticleCut = a.fFirstParticleCut->Clone();
97   // find the right second particle cut
98   if (a.fFirstParticleCut==a.fSecondParticleCut) 
99     SetSecondParticleCut(fFirstParticleCut); // identical particle hbt
100   else
101   fSecondParticleCut = a.fSecondParticleCut->Clone();
102
103   fPairCut = a.fPairCut->Clone();
104   
105   if ( fEventCut ) {
106       SetEventCut(fEventCut); // this will set the myAnalysis pointer inside the cut
107       cout << " AliFemtoAnalysisReactionPlane::AliFemtoAnalysisReactionPlane(const AliFemtoAnalysisReactionPlane& a) - event cut set " << endl;
108   }
109   if ( fFirstParticleCut ) {
110       SetFirstParticleCut(fFirstParticleCut); // this will set the myAnalysis pointer inside the cut
111       cout << " AliFemtoAnalysisReactionPlane::AliFemtoAnalysisReactionPlane(const AliFemtoAnalysisReactionPlane& a) - first particle cut set " << endl;
112   }
113   if ( fSecondParticleCut ) {
114       SetSecondParticleCut(fSecondParticleCut); // this will set the myAnalysis pointer inside the cut
115       cout << " AliFemtoAnalysisReactionPlane::AliFemtoAnalysisReactionPlane(const AliFemtoAnalysisReactionPlane& a) - second particle cut set " << endl;
116   }  if ( fPairCut ) {
117       SetPairCut(fPairCut); // this will set the myAnalysis pointer inside the cut
118       cout << " AliFemtoAnalysisReactionPlane::AliFemtoAnalysisReactionPlane(const AliFemtoAnalysisReactionPlane& a) - pair cut set " << endl;
119   }
120
121   AliFemtoCorrFctnIterator iter;
122   for (iter=a.fCorrFctnCollection->begin(); iter!=a.fCorrFctnCollection->end();iter++){
123     cout << " AliFemtoAnalysisReactionPlane::AliFemtoAnalysisReactionPlane(const AliFemtoAnalysisReactionPlane& a) - looking for correlation functions " << endl;
124     AliFemtoCorrFctn* fctn = (*iter)->Clone();
125     if (fctn) AddCorrFctn(fctn);
126     else cout << " AliFemtoAnalysisReactionPlane::AliFemtoAnalysisReactionPlane(const AliFemtoAnalysisReactionPlane& a) - correlation function not found " << endl;
127   }
128
129   fNumEventsToMix = a.fNumEventsToMix;
130
131   cout << " AliFemtoAnalysisReactionPlane::AliFemtoAnalysisReactionPlane(const AliFemtoAnalysisReactionPlane& a) - analysis copied " << endl;
132
133 }
134 //____________________________
135 AliFemtoAnalysisReactionPlane::~AliFemtoAnalysisReactionPlane(){
136   // now delete every PicoEvent in the EventMixingBuffer and then the Buffer itself
137   delete fPicoEventCollectionVectorHideAway;
138 }
139
140 //____________________________
141 AliFemtoString AliFemtoAnalysisReactionPlane::Report()
142 {
143   // Prepare a report of the execution
144   cout << "AliFemtoAnalysisReactionPlane - constructing Report..."<<endl;
145   char ctemp[200];
146   AliFemtoString temp = "-----------\nHbt AliFemtoAnalysisReactionPlane Report:\n";
147   sprintf(ctemp,"Events are mixed in %d VertexZ bins in the range %E cm to %E cm.\n",fVertexZBins,fVertexZ[0],fVertexZ[1]);
148   temp += ctemp;
149   sprintf(ctemp,"Events underflowing: %d\n",fUnderFlowVertexZ);
150   temp += ctemp;
151   sprintf(ctemp,"Events overflowing: %d\n",fOverFlowVertexZ);
152   temp += ctemp;
153   sprintf(ctemp,"Events are mixed in %d Mult bins in the range %E cm to %E cm.\n",fMultBins,fMult[0],fMult[1]);
154   temp += ctemp;
155   sprintf(ctemp,"Events underflowing: %d\n",fUnderFlowMult);
156   temp += ctemp;
157   sprintf(ctemp,"Events overflowing: %d\n",fOverFlowMult);
158   temp += ctemp;
159   sprintf(ctemp,"Now adding AliFemtoSimpleAnalysis(base) Report\n");
160   temp += ctemp;
161   temp += AliFemtoSimpleAnalysis::Report();
162   AliFemtoString returnThis=temp;
163   return returnThis;
164 }
165 //_________________________
166 void AliFemtoAnalysisReactionPlane::ProcessEvent(const AliFemtoEvent* hbtEvent) {
167   // Perform event processing
168   // in bins of z vertex and multiplicity
169
170   // cout << " AliFemtoAnalysisReactionPlane::ProcessEvent(const AliFemtoEvent* hbtEvent) " << endl;
171   // get right mixing buffer
172   double vertexZ = hbtEvent->PrimVertPos().z();
173   double mult = hbtEvent->UncorrectedNumberOfPrimaries();
174   fCurrentRP = hbtEvent->ReactionPlaneAngle();
175   
176   fMixingBuffer = fPicoEventCollectionVectorHideAway->PicoEventCollection(vertexZ,mult,fCurrentRP); 
177   if (!fMixingBuffer) {
178     if ( vertexZ < fVertexZ[0] ) fUnderFlowVertexZ++;
179     if ( vertexZ > fVertexZ[1] ) fOverFlowVertexZ++;
180     if ( mult < fMult[0] ) fUnderFlowMult++;
181     if ( mult > fMult[1] ) fOverFlowMult++;
182     return;
183   }
184   // call ProcessEvent() from AliFemtoSimpleAnalysis-base
185   AliFemtoSimpleAnalysis::ProcessEvent(hbtEvent);
186 }
187
188 double AliFemtoAnalysisReactionPlane::GetCurrentReactionPlane()
189 {
190   return fCurrentRP;
191 }