1 ////////////////////////////////////////////////////////////////////////////
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 //
7 ////////////////////////////////////////////////////////////////////////////
9 #include "AliFemtoAnalysisAzimuthal.h"
13 #include "AliFemtoParticleCollection.h"
14 #include "AliFemtoTrackCut.h"
15 #include "AliFemtoV0Cut.h"
16 #include "AliFemtoPairCut.h"
18 #include "AliFemtoKinkCut.h"
19 #include "AliFemtoPicoEventCollectionVector.h"
20 #include "AliFemtoPicoEventCollectionVectorHideAway.h"
23 ClassImp(AliFemtoAnalysisAzimuthal)
26 extern void FillHbtParticleCollection(AliFemtoParticleCut* partCut,
27 AliFemtoEvent* hbtEvent,
28 AliFemtoParticleCollection* partCollection);
31 //____________________________
32 AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(unsigned int binsVertex, double minVertex, double maxVertex,
33 unsigned int binsMult, double minMult, double maxMult, unsigned short binsRP)
39 fVertexZBins(binsVertex),
48 // mControlSwitch = 0;
49 fCorrFctnCollection= 0;
50 fCorrFctnCollection = new AliFemtoCorrFctnCollection;
51 fVertexZ[0] = minVertex;
52 fVertexZ[1] = maxVertex;
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());
61 //____________________________
63 AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) :
64 AliFemtoSimpleAnalysis(),
69 fVertexZBins(a.fVertexZBins),
72 fMultBins(a.fMultBins) ,
79 fCorrFctnCollection= 0;
80 fCorrFctnCollection = new AliFemtoCorrFctnCollection;
81 fVertexZ[0] = a.fVertexZ[0];
82 fVertexZ[1] = a.fVertexZ[1];
83 fMult[0] = a.fMult[0];
84 fMult[1] = a.fMult[1];
85 if (fMixingBuffer) delete fMixingBuffer;
86 fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
87 fMultBins,fMult[0],fMult[1],
88 fRPBins,0.0,TMath::Pi());
89 // find the right event cut
90 fEventCut = a.fEventCut->Clone();
91 // find the right femto particle cut
92 fFemtoParticleCut = a.fFemtoParticleCut->Clone();
93 // find the right flow particle cut
94 fFlowParticleCut = a.fFlowParticleCut->Clone();
95 // find the right pair cut
96 fPairCut = a.fPairCut->Clone();
99 SetEventCut(fEventCut); // this will set the myAnalysis pointer inside the cut
100 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - event cut set " << endl;
102 if ( fFemtoParticleCut ) {
103 SetFirstParticleCut(fFemtoParticleCut); // this will set the myAnalysis pointer inside the cut
104 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - femto particle cut set " << endl;
106 if ( fFlowParticleCut ) {
107 SetSecondParticleCut(fFlowParticleCut); // this will set the myAnalysis pointer inside the cut
108 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - flow particle cut set " << endl;
111 SetPairCut(fPairCut); // this will set the myAnalysis pointer inside the cut
112 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - pair cut set " << endl;
115 AliFemtoCorrFctnIterator iter;
116 for (iter=a.fCorrFctnCollection->begin(); iter!=a.fCorrFctnCollection->end();iter++){
117 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - looking for correlation functions " << endl;
118 AliFemtoCorrFctn* fctn = (*iter)->Clone();
119 if (fctn) AddCorrFctn(fctn);
120 else cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - correlation function not found " << endl;
122 fNumEventsToMix = a.fNumEventsToMix;
123 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - analysis copied " << endl;
126 AliFemtoAnalysisAzimuthal& AliFemtoAnalysisAzimuthal::operator=(const AliFemtoAnalysisAzimuthal& a)
128 // Assignment operator
132 fCorrFctnCollection= 0;
133 fCorrFctnCollection = new AliFemtoCorrFctnCollection;
134 fVertexZ[0] = a.fVertexZ[0];
135 fVertexZ[1] = a.fVertexZ[1];
136 fMult[0] = a.fMult[0];
137 fMult[1] = a.fMult[1];
138 if (fMixingBuffer) delete fMixingBuffer;
139 fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
140 fMultBins,fMult[0],fMult[1],
141 fRPBins,0.0,TMath::Pi());
142 // find the right event cut
143 fEventCut = a.fEventCut->Clone();
144 // find the right femto particle cut
145 fFemtoParticleCut = a.fFemtoParticleCut->Clone();
146 // find the right flow particle cut
147 fFlowParticleCut = a.fFlowParticleCut->Clone();
148 // find the right pair cut
149 fPairCut = a.fPairCut->Clone();
152 SetEventCut(fEventCut); // this will set the myAnalysis pointer inside the cut
153 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - event cut set " << endl;
155 if ( fFemtoParticleCut ) {
156 SetFirstParticleCut(fFemtoParticleCut); // this will set the myAnalysis pointer inside the cut
157 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - femto particle cut set " << endl;
159 if ( fFlowParticleCut ) {
160 SetSecondParticleCut(fFlowParticleCut); // this will set the myAnalysis pointer inside the cut
161 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - flow particle cut set " << endl;
164 SetPairCut(fPairCut); // this will set the myAnalysis pointer inside the cut
165 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - pair cut set " << endl;
168 AliFemtoCorrFctnIterator iter;
169 for (iter=a.fCorrFctnCollection->begin(); iter!=a.fCorrFctnCollection->end();iter++){
170 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - looking for correlation functions " << endl;
171 AliFemtoCorrFctn* fctn = (*iter)->Clone();
172 if (fctn) AddCorrFctn(fctn);
173 else cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - correlation function not found " << endl;
175 fNumEventsToMix = a.fNumEventsToMix;
176 cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - analysis copied " << endl;
182 //____________________________
183 AliFemtoAnalysisAzimuthal::~AliFemtoAnalysisAzimuthal(){
184 // now delete every PicoEvent in the EventMixingBuffer and then the Buffer itself
185 delete fPicoEventCollectionVectorHideAway;
188 //_________________________
189 void AliFemtoAnalysisAzimuthal::ProcessEvent(const AliFemtoEvent* hbtEvent) {
190 // Perform event processing in bins of z vertex, multiplicity and Reaction Plane angle
191 // cout << " AliFemtoAnalysisAzimuthal::ProcessEvent(const AliFemtoEvent* hbtEvent) " << endl;
193 // get right mixing buffer
194 double vertexZ = hbtEvent->PrimVertPos().z();
195 double mult = hbtEvent->UncorrectedNumberOfPrimaries();
196 //double rpangle = GetQVector() // Not possible at the moment to get rpangle for mixing wrt RP
198 fMixingBuffer = fPicoEventCollectionVectorHideAway->PicoEventCollection(vertexZ,mult,fPsi);
199 if (!fMixingBuffer) {
200 if ( vertexZ < fVertexZ[0] ) fUnderFlowVertexZ++;
201 if ( vertexZ > fVertexZ[1] ) fOverFlowVertexZ++;
202 if ( mult < fMult[0] ) fUnderFlowMult++;
203 if ( mult > fMult[1] ) fOverFlowMult++;
207 //****from AliFemtoSimpleAnalysis****
209 fPicoEvent=0; // we will get a new pico event, if not prevent corr. fctn to access old pico event
213 fFemtoParticleCut->EventBegin(hbtEvent);
214 fFlowParticleCut->EventBegin(hbtEvent);
215 fPairCut->EventBegin(hbtEvent);
217 for (AliFemtoCorrFctnIterator iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
218 (*iter)->EventBegin(hbtEvent);
221 // event cut and event cut monitor
222 bool tmpPassEvent = fEventCut->Pass(hbtEvent);
224 fEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
226 fPicoEvent = new AliFemtoPicoEvent; // this is what we will make pairs from and put in Mixing Buffer, no memory leak. we will delete picoevents when they come out of the mixing buffer
227 FillHbtParticleCollection(fFemtoParticleCut,(AliFemtoEvent*)hbtEvent,fPicoEvent->FirstParticleCollection());
228 FillHbtParticleCollection(fFlowParticleCut,(AliFemtoEvent*)hbtEvent,fPicoEvent->SecondParticleCollection());
230 cout << "#particles in Collection: " << fPicoEvent->FirstParticleCollection()->size() << endl;
232 //switch which allows only using events with ParticleCollections containing a minimum number of entries
233 if (fPicoEvent->FirstParticleCollection()->size() >= fMinSizePartCollection ) {
234 fEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
237 //------ Make real pairs (assume identical) ------//
238 MakePairs("real", fPicoEvent->FirstParticleCollection() );
239 cout << "AliFemtoAnalysisAzimuthal::ProcessEvent() - reals done ";
241 //---- Make pairs for mixed events, looping over events in mixingBuffer ----//
242 AliFemtoPicoEvent* storedEvent;
243 AliFemtoPicoEventIterator fPicoEventIter;
244 for (fPicoEventIter=MixingBuffer()->begin();fPicoEventIter!=MixingBuffer()->end();fPicoEventIter++){
245 storedEvent = *fPicoEventIter;
246 MakePairs("mixed",fPicoEvent->FirstParticleCollection(),
247 storedEvent->FirstParticleCollection() );
249 cout << " - mixed done " << endl;
251 //--------- If mixing buffer is full, delete oldest event ---------//
252 if ( MixingBufferFull() ) {
253 delete MixingBuffer()->back();
254 MixingBuffer()->pop_back();
257 //-------- Add current event (fPicoEvent) to mixing buffer --------//
258 MixingBuffer()->push_front(fPicoEvent);
260 } // if ParticleCollections are big enough (mal jun2002)
262 fEventCut->FillCutMonitor(hbtEvent, !tmpPassEvent);
266 // if currentEvent is accepted by currentAnalysis cleanup for EbyE
267 fFemtoParticleCut->EventEnd(hbtEvent);
268 fFlowParticleCut->EventEnd(hbtEvent);
269 fPairCut->EventEnd(hbtEvent);
270 for (AliFemtoCorrFctnIterator iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
271 (*iter)->EventEnd(hbtEvent);
275 //_______________________________________________________________________________
276 void AliFemtoAnalysisAzimuthal::MakePairs(const char* typeIn, AliFemtoParticleCollection *partCollection1,
277 AliFemtoParticleCollection *partCollection2){
278 // Build pairs, check pair cuts, and call CFs' AddRealPair() or AddMixedPair() methods.
279 // If no second particle collection is specfied, make pairs within first particle collection.
281 string type = typeIn;
283 // int swpart = ((long int) partCollection1) % 2;
284 int swpart = fNeventsProcessed % 2;
286 AliFemtoPair* tPair = new AliFemtoPair;
287 AliFemtoCorrFctnIterator tCorrFctnIter;
288 AliFemtoParticleIterator tPartIter1, tPartIter2;
290 AliFemtoParticleIterator tStartOuterLoop = partCollection1->begin(); // always
291 AliFemtoParticleIterator tEndOuterLoop = partCollection1->end(); // will be one less if identical
292 AliFemtoParticleIterator tStartInnerLoop;
293 AliFemtoParticleIterator tEndInnerLoop;
294 if (partCollection2) { // Two collections:
295 tStartInnerLoop = partCollection2->begin(); // Full inner & outer loops
296 tEndInnerLoop = partCollection2->end();
298 else { // One collection:
299 tEndOuterLoop--; // Outer loop goes to next-to-last particle
300 tEndInnerLoop = partCollection1->end() ; // Inner loop goes to last particle
302 for (tPartIter1=tStartOuterLoop;tPartIter1!=tEndOuterLoop;tPartIter1++) {
303 if (!partCollection2){
304 tStartInnerLoop = tPartIter1;
307 tPair->SetTrack1(*tPartIter1);
308 for (tPartIter2 = tStartInnerLoop; tPartIter2!=tEndInnerLoop;tPartIter2++) {
309 tPair->SetTrack2(*tPartIter2);
311 //For getting the pair angle wrt EP
314 float Psi=0, mQx=0, mQy=0, PairAngleEP=0;
315 Q = GetQVector(partCollection1);
317 mQx = Q.X() - cos(2*(tPair->Track1()->FourMomentum().Phi()))*(tPair->Track1()->Track()->Pt()) - cos(2*(tPair->Track2()->FourMomentum().Phi()))*(tPair->Track2()->Track()->Pt());
318 mQy = Q.Y() - sin(2*(tPair->Track1()->FourMomentum().Phi()))*(tPair->Track1()->Track()->Pt()) - sin(2*(tPair->Track2()->FourMomentum().Phi()))*(tPair->Track2()->Track()->Pt());
320 QVector.Set(mQx,mQy);
322 Psi=QVector.Phi()/2.;
323 if (Psi > TMath::Pi()) Psi -= TMath::Pi();
325 PairAngleEP = (tPair->EmissionAngle() - Psi);
326 tPair->SetPairAngleEP(PairAngleEP);
329 if (type == "mixed"){
330 float Psi1=0, Psi2=0, mQx1=0, mQx2=0,mQy1=0, mQy2=0, px1=0, px2=0, py1=0, py2=0, PairAngleEP=0;
331 TVector2 Q1, Q2, QVector1, QVector2, P;
333 Q1 = GetQVector(partCollection1);
334 Q2 = GetQVector(partCollection2);
336 mQx1 = Q1.X() - cos(2*(tPair->Track1()->FourMomentum().Phi()))*(tPair->Track1()->Track()->Pt());
337 mQx2 = Q2.X() - cos(2*(tPair->Track2()->FourMomentum().Phi()))*(tPair->Track2()->Track()->Pt());
338 mQy1 = Q1.Y() - sin(2*(tPair->Track1()->FourMomentum().Phi()))*(tPair->Track1()->Track()->Pt());
339 mQy2 = Q2.Y() - sin(2*(tPair->Track2()->FourMomentum().Phi()))*(tPair->Track2()->Track()->Pt());
341 QVector1.Set(mQx1,mQy1);
342 QVector2.Set(mQx2,mQy2);
344 Psi1=QVector1.Phi()/2.;
345 if (Psi1 > TMath::Pi()) Psi1 -= TMath::Pi();
347 Psi2=QVector2.Phi()/2.;
348 if (Psi2 > TMath::Pi()) Psi2 -= TMath::Pi();
350 px1 = (tPair->Track1()->Track()->Pt())*cos(tPair->Track1()->FourMomentum().Phi() - Psi1);
351 px2 = (tPair->Track2()->Track()->Pt())*cos(tPair->Track2()->FourMomentum().Phi() - Psi2);
352 py1 = (tPair->Track1()->Track()->Pt())*sin(tPair->Track1()->FourMomentum().Phi() - Psi1);
353 py2 = (tPair->Track2()->Track()->Pt())*sin(tPair->Track2()->FourMomentum().Phi() - Psi2);
355 P.Set(px1+px2, py1+py2);
356 PairAngleEP = P.Phi();
358 tPair->SetPairAngleEP(PairAngleEP);
361 // The following lines have to be uncommented if you want pairCutMonitors
362 // they are not in for speed reasons
363 // bool tmpPassPair = fPairCut->Pass(tPair);
364 // fPairCut->FillCutMonitor(tPair, tmpPassPair);
365 // if ( tmpPassPair )
367 //---- If pair passes cut, loop over CF's and add pair to real/mixed ----//
369 if (!partCollection2) {
371 tPair->SetTrack1(*tPartIter2);
372 tPair->SetTrack2(*tPartIter1);
376 tPair->SetTrack1(*tPartIter1);
377 tPair->SetTrack2(*tPartIter2);
383 if (fPairCut->Pass(tPair)){
384 for (tCorrFctnIter=fCorrFctnCollection->begin();
385 tCorrFctnIter!=fCorrFctnCollection->end();tCorrFctnIter++){
386 AliFemtoCorrFctn* tCorrFctn = *tCorrFctnIter;
388 tCorrFctn->AddRealPair(tPair);
389 else if(type == "mixed")
390 tCorrFctn->AddMixedPair(tPair);
392 cout << "Problem with pair type, type = " << type.c_str() << endl;
395 } // loop over second particle
396 } // loop over first particle
401 //_____________________________________________
402 TVector2 AliFemtoAnalysisAzimuthal::GetQVector(AliFemtoParticleCollection* particlecollection){
407 AliFemtoParticle* flowparticle;
408 AliFemtoParticleIterator pIter;
409 AliFemtoParticleIterator startLoop = particlecollection->begin();
410 AliFemtoParticleIterator endLoop = particlecollection->end();
411 for (pIter=startLoop;pIter!=endLoop;pIter++){
412 flowparticle = *pIter;
413 mQx += (cos(2*flowparticle->FourMomentum().Phi()))*(flowparticle->Track()->Pt());
414 mQy += (sin(2*flowparticle->FourMomentum().Phi()))*(flowparticle->Track()->Pt());
421 //__________________________________________________
422 double AliFemtoAnalysisAzimuthal::GetCurrentReactionPlane()
427 //_________________________
428 TList* AliFemtoAnalysisAzimuthal::GetOutputList()
430 // Collect the list of output objects to be written
432 TList *tOutputList = new TList();
433 TList *p1Cut = fFemtoParticleCut->GetOutputList();
435 TListIter nextp1(p1Cut);
436 while (TObject *obj = nextp1.Next()) {
437 tOutputList->Add(obj);
440 TList *pairCut = fPairCut->GetOutputList();
442 TIter nextpair(pairCut);
443 while (TObject *obj = nextpair()) {
444 tOutputList->Add(obj);
447 TList *eventCut = fEventCut->GetOutputList();
449 TIter nextevent(eventCut);
450 while (TObject *obj = nextevent()) {
451 tOutputList->Add(obj);
454 AliFemtoCorrFctnIterator iter;
455 for (iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
456 TList *tListCf = (*iter)->GetOutputList();
458 TIter nextListCf(tListCf);
459 while (TObject *obj = nextListCf()) {
460 tOutputList->Add(obj);