[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemto / AliFemtoAnalysisAzimuthal.cxx

////////////////////////////////////////////////////////////////////////////
//                                                                        //
// AliFemtoAnalysisReactionPlane - Femtoscopic analysis which mixes event //
// with respect to the z position of the primary vertex and event total   //
// multiplicity and uses only events in certain reaction plane angle bin  //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include "AliFemtoAnalysisAzimuthal.h"
#include <TMath.h>
#include <string>
#include <cstdio>
#include "AliFemtoParticleCollection.h"
#include "AliFemtoTrackCut.h"
#include "AliFemtoV0Cut.h"
#include "AliFemtoPairCut.h"
#include "TVector2.h"
#include "AliFemtoKinkCut.h"
#include "AliFemtoPicoEventCollectionVector.h"
#include "AliFemtoPicoEventCollectionVectorHideAway.h"

#ifdef __ROOT__ 
ClassImp(AliFemtoAnalysisAzimuthal)
#endif

extern void FillHbtParticleCollection(AliFemtoParticleCut*         partCut,
				      AliFemtoEvent*               hbtEvent,
				      AliFemtoParticleCollection*  partCollection);


//____________________________
AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(unsigned int binsVertex, double minVertex, double maxVertex,
						       unsigned int binsMult, double minMult, double maxMult, unsigned short binsRP) 
  : 
  fFemtoParticleCut(0),
  fFlowParticleCut(0),
  fVertexZBins(binsVertex), 
  fOverFlowVertexZ(0), 
  fUnderFlowVertexZ(0),
  fMultBins(binsMult) ,
  fOverFlowMult(0),    
  fUnderFlowMult(0),
  fRPBins(binsRP),
  fPsi(0)
{
  //  mControlSwitch     = 0;
  fCorrFctnCollection= 0;
  fCorrFctnCollection = new AliFemtoCorrFctnCollection;
  fVertexZ[0] = minVertex;
  fVertexZ[1] = maxVertex;
  fMult[0] = minMult;
  fMult[1] = maxMult;
  if (fMixingBuffer) delete fMixingBuffer;
  fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
										     fMultBins,fMult[0],fMult[1],
										     fRPBins,0.0,TMath::Pi());
}

//____________________________

AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) : 
  AliFemtoSimpleAnalysis(),
  fFemtoParticleCut(0),
  fFlowParticleCut(0),
  fVertexZBins(a.fVertexZBins), 
  fOverFlowVertexZ(0), 
  fUnderFlowVertexZ(0),
  fMultBins(a.fMultBins) ,
  fOverFlowMult(0),    
  fUnderFlowMult(0),
  fRPBins(a.fRPBins),
  fPsi(0)

{
  fCorrFctnCollection= 0;
  fCorrFctnCollection = new AliFemtoCorrFctnCollection;
  fVertexZ[0] = a.fVertexZ[0]; 
  fVertexZ[1] = a.fVertexZ[1];
  fMult[0] = a.fMult[0]; 
  fMult[1] = a.fMult[1];
  if (fMixingBuffer) delete fMixingBuffer;
  fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
										     fMultBins,fMult[0],fMult[1],
										     fRPBins,0.0,TMath::Pi());
  // find the right event cut
  fEventCut = a.fEventCut->Clone();
  // find the right femto particle cut
  fFemtoParticleCut = a.fFemtoParticleCut->Clone();
  // find the right flow particle cut
  fFlowParticleCut = a.fFlowParticleCut->Clone();
  // find the right pair cut
  fPairCut = a.fPairCut->Clone();
  
  if ( fEventCut ) {
      SetEventCut(fEventCut); // this will set the myAnalysis pointer inside the cut
      cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - event cut set " << endl;
  }
  if ( fFemtoParticleCut ) {
      SetFirstParticleCut(fFemtoParticleCut); // this will set the myAnalysis pointer inside the cut
      cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - femto particle cut set " << endl;
  }
  if ( fFlowParticleCut ) {
      SetSecondParticleCut(fFlowParticleCut); // this will set the myAnalysis pointer inside the cut
      cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - flow particle cut set " << endl;
  }
  if ( fPairCut ) {
      SetPairCut(fPairCut); // this will set the myAnalysis pointer inside the cut
      cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - pair cut set " << endl;
  }

  AliFemtoCorrFctnIterator iter;
  for (iter=a.fCorrFctnCollection->begin(); iter!=a.fCorrFctnCollection->end();iter++){
    cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - looking for correlation functions " << endl;
    AliFemtoCorrFctn* fctn = (*iter)->Clone();
    if (fctn) AddCorrFctn(fctn);
    else cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - correlation function not found " << endl;
  }
  fNumEventsToMix = a.fNumEventsToMix;
  cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - analysis copied " << endl;
}

AliFemtoAnalysisAzimuthal& AliFemtoAnalysisAzimuthal::operator=(const AliFemtoAnalysisAzimuthal& a)
{
  // Assignment operator
  if (this == &a)
    return *this;

  fCorrFctnCollection= 0;
  fCorrFctnCollection = new AliFemtoCorrFctnCollection;
  fVertexZ[0] = a.fVertexZ[0]; 
  fVertexZ[1] = a.fVertexZ[1];
  fMult[0] = a.fMult[0]; 
  fMult[1] = a.fMult[1];
  if (fMixingBuffer) delete fMixingBuffer;
  fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
										     fMultBins,fMult[0],fMult[1],
										     fRPBins,0.0,TMath::Pi());
  // find the right event cut
  fEventCut = a.fEventCut->Clone();
  // find the right femto particle cut
  fFemtoParticleCut = a.fFemtoParticleCut->Clone();
  // find the right flow particle cut
  fFlowParticleCut = a.fFlowParticleCut->Clone();
  // find the right pair cut
  fPairCut = a.fPairCut->Clone();
  
  if ( fEventCut ) {
      SetEventCut(fEventCut); // this will set the myAnalysis pointer inside the cut
      cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - event cut set " << endl;
  }
  if ( fFemtoParticleCut ) {
      SetFirstParticleCut(fFemtoParticleCut); // this will set the myAnalysis pointer inside the cut
      cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - femto particle cut set " << endl;
  }
  if ( fFlowParticleCut ) {
      SetSecondParticleCut(fFlowParticleCut); // this will set the myAnalysis pointer inside the cut
      cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - flow particle cut set " << endl;
  }
  if ( fPairCut ) {
      SetPairCut(fPairCut); // this will set the myAnalysis pointer inside the cut
      cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - pair cut set " << endl;
  }

  AliFemtoCorrFctnIterator iter;
  for (iter=a.fCorrFctnCollection->begin(); iter!=a.fCorrFctnCollection->end();iter++){
    cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - looking for correlation functions " << endl;
    AliFemtoCorrFctn* fctn = (*iter)->Clone();
    if (fctn) AddCorrFctn(fctn);
    else cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - correlation function not found " << endl;
  }
  fNumEventsToMix = a.fNumEventsToMix;
  cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - analysis copied " << endl;

  return *this;
  
}

//____________________________
AliFemtoAnalysisAzimuthal::~AliFemtoAnalysisAzimuthal(){
  // now delete every PicoEvent in the EventMixingBuffer and then the Buffer itself
  delete fPicoEventCollectionVectorHideAway;
}

//_________________________
void AliFemtoAnalysisAzimuthal::ProcessEvent(const AliFemtoEvent* hbtEvent) {
  // Perform event processing in bins of z vertex, multiplicity and Reaction Plane angle
  // cout << " AliFemtoAnalysisAzimuthal::ProcessEvent(const AliFemtoEvent* hbtEvent) " << endl;

  //****from AliFemtoSimpleAnalysis****

  fPicoEvent=0; // we will get a new pico event, if not prevent corr. fctn to access old pico event
  fNeventsProcessed++;

  // startup for EbyE 
  fFemtoParticleCut->EventBegin(hbtEvent);
  fFlowParticleCut->EventBegin(hbtEvent);
  fPairCut->EventBegin(hbtEvent);

  for (AliFemtoCorrFctnIterator iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
    (*iter)->EventBegin(hbtEvent);
  }

  // event cut and event cut monitor
  bool tmpPassEvent = fEventCut->Pass(hbtEvent);
  if (!tmpPassEvent) 
    fEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
  if (tmpPassEvent) {
    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
    FillHbtParticleCollection(fFemtoParticleCut,(AliFemtoEvent*)hbtEvent,fPicoEvent->FirstParticleCollection());
    FillHbtParticleCollection(fFlowParticleCut,(AliFemtoEvent*)hbtEvent,fPicoEvent->SecondParticleCollection());
    
      // get right mixing buffer
  double vertexZ = hbtEvent->PrimVertPos().z();
  double mult = hbtEvent->UncorrectedNumberOfPrimaries();
  TVector2 tQ = GetQVector(fPicoEvent->SecondParticleCollection());
  double tPsi=tQ.Phi()/2.;
  if (tPsi > TMath::Pi()) tPsi -= TMath::Pi();
      
   fMixingBuffer = fPicoEventCollectionVectorHideAway->PicoEventCollection(vertexZ,mult,fPsi); 
  if (!fMixingBuffer) {
    if ( vertexZ < fVertexZ[0] ) fUnderFlowVertexZ++;
    if ( vertexZ > fVertexZ[1] ) fOverFlowVertexZ++;
    if ( mult < fMult[0] ) fUnderFlowMult++;
    if ( mult > fMult[1] ) fOverFlowMult++;
    return;
  }   
    
    //cout << "#particles in Collection: " << fPicoEvent->FirstParticleCollection()->size() << endl;
    
    //switch which allows only using events with ParticleCollections containing a minimum number of entries
    if (fPicoEvent->FirstParticleCollection()->size() >= fMinSizePartCollection ) {
      fEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
    }

      //------ Make real pairs (assume identical) ------//
      MakePairs("real", fPicoEvent->FirstParticleCollection() );
      //cout << "AliFemtoAnalysisAzimuthal::ProcessEvent() - reals done ";

      //---- Make pairs for mixed events, looping over events in mixingBuffer ----//
      AliFemtoPicoEvent* storedEvent;
      AliFemtoPicoEventIterator fPicoEventIter;
      for (fPicoEventIter=MixingBuffer()->begin();fPicoEventIter!=MixingBuffer()->end();fPicoEventIter++){
        storedEvent = *fPicoEventIter;
        MakePairs("mixed",fPicoEvent->FirstParticleCollection(),
                            storedEvent->FirstParticleCollection() );
      }
      //cout << " - mixed done   " << endl;

      //--------- If mixing buffer is full, delete oldest event ---------//
      if ( MixingBufferFull() ) {
        delete MixingBuffer()->back();
        MixingBuffer()->pop_back();
      }

      //-------- Add current event (fPicoEvent) to mixing buffer --------//
      MixingBuffer()->push_front(fPicoEvent);

    }  // if ParticleCollections are big enough (mal jun2002)
    else{
      fEventCut->FillCutMonitor(hbtEvent, !tmpPassEvent);
      delete fPicoEvent;
    }

  // if currentEvent is accepted by currentAnalysis cleanup for EbyE 
  fFemtoParticleCut->EventEnd(hbtEvent);
  fFlowParticleCut->EventEnd(hbtEvent);
  fPairCut->EventEnd(hbtEvent);
  for (AliFemtoCorrFctnIterator iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
    (*iter)->EventEnd(hbtEvent);
  } 
}

//_______________________________________________________________________________
void AliFemtoAnalysisAzimuthal::MakePairs(const char* typeIn, AliFemtoParticleCollection *partCollection1,
				       AliFemtoParticleCollection *partCollection2){
// Build pairs, check pair cuts, and call CFs' AddRealPair() or AddMixedPair() methods. 
// If no second particle collection is specfied, make pairs within first particle collection.

  string type = typeIn;

  //  int swpart = ((long int) partCollection1) % 2;
  int swpart = fNeventsProcessed % 2;

  AliFemtoPair* tPair = new AliFemtoPair;
  AliFemtoCorrFctnIterator tCorrFctnIter;
  AliFemtoParticleIterator tPartIter1, tPartIter2;

  AliFemtoParticleIterator tStartOuterLoop = partCollection1->begin();  // always
  AliFemtoParticleIterator tEndOuterLoop   = partCollection1->end();    // will be one less if identical
  AliFemtoParticleIterator tStartInnerLoop;
  AliFemtoParticleIterator tEndInnerLoop;
  if (partCollection2) {                         //   Two collections:
    tStartInnerLoop = partCollection2->begin();  //   Full inner & outer loops
    tEndInnerLoop   = partCollection2->end();    
  }
  else {                                         //   One collection:
    tEndOuterLoop--;                             //   Outer loop goes to next-to-last particle
    tEndInnerLoop = partCollection1->end() ;     //   Inner loop goes to last particle
  }
  for (tPartIter1=tStartOuterLoop;tPartIter1!=tEndOuterLoop;tPartIter1++) {
    if (!partCollection2){
      tStartInnerLoop = tPartIter1;
      tStartInnerLoop++;
    }
    tPair->SetTrack1(*tPartIter1);
    for (tPartIter2 = tStartInnerLoop; tPartIter2!=tEndInnerLoop;tPartIter2++) {
      tPair->SetTrack2(*tPartIter2);

//For getting the pair angle wrt EP
  if (type == "real"){
	TVector2 tQ, tQVector;
	float tPsi=0, mQx=0, mQy=0, PairAngleEP=0;
	tQ = GetQVector(partCollection1);

	mQx = tQ.X() - cos(2*(tPair->Track1()->FourMomentum().Phi()))*(tPair->Track1()->Track()->Pt()) - cos(2*(tPair->Track2()->FourMomentum().Phi()))*(tPair->Track2()->Track()->Pt());
	mQy = tQ.Y() - sin(2*(tPair->Track1()->FourMomentum().Phi()))*(tPair->Track1()->Track()->Pt()) - sin(2*(tPair->Track2()->FourMomentum().Phi()))*(tPair->Track2()->Track()->Pt());
  
	tQVector.Set(mQx,mQy);

	tPsi=tQVector.Phi()/2.;
	if (tPsi > TMath::Pi()) tPsi -= TMath::Pi();

	PairAngleEP = (tPair->EmissionAngle() - tPsi);
	tPair->SetPairAngleEP(PairAngleEP);
  }

  if (type == "mixed"){
	float tPsi1=0, tPsi2=0, mQx1=0, mQx2=0,mQy1=0, mQy2=0, px1=0, px2=0, py1=0, py2=0, PairAngleEP=0;
	TVector2 tQ1, tQ2, tQVector1, tQVector2, tP;

	tQ1 = GetQVector(partCollection1);
	tQ2 = GetQVector(partCollection2);

	mQx1 = tQ1.X() - cos(2*(tPair->Track1()->FourMomentum().Phi()))*(tPair->Track1()->Track()->Pt());
	mQx2 = tQ2.X() - cos(2*(tPair->Track2()->FourMomentum().Phi()))*(tPair->Track2()->Track()->Pt());
	mQy1 = tQ1.Y() - sin(2*(tPair->Track1()->FourMomentum().Phi()))*(tPair->Track1()->Track()->Pt());
	mQy2 = tQ2.Y() - sin(2*(tPair->Track2()->FourMomentum().Phi()))*(tPair->Track2()->Track()->Pt());
  
	tQVector1.Set(mQx1,mQy1);
	tQVector2.Set(mQx2,mQy2);

	tPsi1=tQVector1.Phi()/2.;
	if (tPsi1 > TMath::Pi()) tPsi1 -= TMath::Pi();

	tPsi2=tQVector2.Phi()/2.;
	if (tPsi2 > TMath::Pi()) tPsi2 -= TMath::Pi();

	px1 = (tPair->Track1()->Track()->Pt())*cos(tPair->Track1()->FourMomentum().Phi() - tPsi1);
	px2 = (tPair->Track2()->Track()->Pt())*cos(tPair->Track2()->FourMomentum().Phi() - tPsi2);
	py1 = (tPair->Track1()->Track()->Pt())*sin(tPair->Track1()->FourMomentum().Phi() - tPsi1);
	py2 = (tPair->Track2()->Track()->Pt())*sin(tPair->Track2()->FourMomentum().Phi() - tPsi2);

	tP.Set(px1+px2, py1+py2);
	PairAngleEP = tP.Phi();

	tPair->SetPairAngleEP(PairAngleEP);
  }

      // The following lines have to be uncommented if you want pairCutMonitors
      // they are not in for speed reasons
      // bool tmpPassPair = fPairCut->Pass(tPair);
      // fPairCut->FillCutMonitor(tPair, tmpPassPair);
      // if ( tmpPassPair )

      //---- If pair passes cut, loop over CF's and add pair to real/mixed ----//

      if (!partCollection2) {
	if (swpart) {
 	  tPair->SetTrack1(*tPartIter2);
 	  tPair->SetTrack2(*tPartIter1);
 	  swpart = 0;
 	}
 	else {
 	  tPair->SetTrack1(*tPartIter1);
 	  tPair->SetTrack2(*tPartIter2);
 	  swpart = 1;
	}
      }


      if (fPairCut->Pass(tPair)){
        for (tCorrFctnIter=fCorrFctnCollection->begin();
             tCorrFctnIter!=fCorrFctnCollection->end();tCorrFctnIter++){
          AliFemtoCorrFctn* tCorrFctn = *tCorrFctnIter;
          if(type == "real")
            tCorrFctn->AddRealPair(tPair);
	  else if(type == "mixed")
            tCorrFctn->AddMixedPair(tPair);
          else
            cout << "Problem with pair type, type = " << type.c_str() << endl;
        }
      }
    }    // loop over second particle
  }      // loop over first particle

  delete tPair;
}

//_____________________________________________
TVector2 AliFemtoAnalysisAzimuthal::GetQVector(AliFemtoParticleCollection* particlecollection){
  
  TVector2 mQ;
  float mQx=0, mQy=0;

  if (!particlecollection) {
    mQ.Set(0.0, 0.0);
    return mQ;
  }

  AliFemtoParticle* flowparticle;
  AliFemtoParticleIterator pIter;
  AliFemtoParticleIterator startLoop = particlecollection->begin();
  AliFemtoParticleIterator endLoop   = particlecollection->end();
  for (pIter=startLoop;pIter!=endLoop;pIter++){
    flowparticle = *pIter;
    mQx += (cos(2*flowparticle->FourMomentum().Phi()))*(flowparticle->Track()->Pt());
    mQy += (sin(2*flowparticle->FourMomentum().Phi()))*(flowparticle->Track()->Pt());
  }
  
  mQ.Set(mQx,mQy);
  return mQ;
}

//__________________________________________________
double AliFemtoAnalysisAzimuthal::GetCurrentReactionPlane()
{
  return fPsi;
}

//_________________________
TList* AliFemtoAnalysisAzimuthal::GetOutputList()
{
  // Collect the list of output objects to be written

  TList *tOutputList = new TList();
  TList *p1Cut = fFemtoParticleCut->GetOutputList();

  TListIter nextp1(p1Cut);
  while (TObject *obj = nextp1.Next()) {
    tOutputList->Add(obj);
  }

  TList *pairCut = fPairCut->GetOutputList();

  TIter nextpair(pairCut);
  while (TObject *obj = nextpair()) {
    tOutputList->Add(obj);
  }

  TList *eventCut = fEventCut->GetOutputList();

  TIter nextevent(eventCut);
  while (TObject *obj = nextevent()) {
    tOutputList->Add(obj);
  }

  AliFemtoCorrFctnIterator iter;
  for (iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
    TList *tListCf = (*iter)->GetOutputList();
    
    TIter nextListCf(tListCf);
    while (TObject *obj = nextListCf()) {
      tOutputList->Add(obj);
    }
  }

  return tOutputList;
  
}
Commit	Line	Data
76ce4b5b	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 "AliFemtoAnalysisAzimuthal.h"
	10	#include <TMath.h>
	11	#include <string>
	12	#include <cstdio>
	13	#include "AliFemtoParticleCollection.h"
	14	#include "AliFemtoTrackCut.h"
	15	#include "AliFemtoV0Cut.h"
	16	#include "AliFemtoPairCut.h"
	17	#include "TVector2.h"
	18	#include "AliFemtoKinkCut.h"
	19	#include "AliFemtoPicoEventCollectionVector.h"
	20	#include "AliFemtoPicoEventCollectionVectorHideAway.h"
	21
	22	#ifdef __ROOT__
	23	ClassImp(AliFemtoAnalysisAzimuthal)
	24	#endif
	25
	26	extern void FillHbtParticleCollection(AliFemtoParticleCut* partCut,
	27	AliFemtoEvent* hbtEvent,
	28	AliFemtoParticleCollection* partCollection);
	29
	30
	31	//____________________________
	32	AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(unsigned int binsVertex, double minVertex, double maxVertex,
	33	unsigned int binsMult, double minMult, double maxMult, unsigned short binsRP)
	34	:
	35	fFemtoParticleCut(0),
	36	fFlowParticleCut(0),
	37	fVertexZBins(binsVertex),
	38	fOverFlowVertexZ(0),
	39	fUnderFlowVertexZ(0),
	40	fMultBins(binsMult) ,
	41	fOverFlowMult(0),
	42	fUnderFlowMult(0),
	43	fRPBins(binsRP),
	44	fPsi(0)
	45	{
	46	// mControlSwitch = 0;
	47	fCorrFctnCollection= 0;
	48	fCorrFctnCollection = new AliFemtoCorrFctnCollection;
	49	fVertexZ[0] = minVertex;
	50	fVertexZ[1] = maxVertex;
	51	fMult[0] = minMult;
	52	fMult[1] = maxMult;
	53	if (fMixingBuffer) delete fMixingBuffer;
	54	fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
	55	fMultBins,fMult[0],fMult[1],
	56	fRPBins,0.0,TMath::Pi());
	57	}
	58
	59	//____________________________
	60
	61	AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) :
	62	AliFemtoSimpleAnalysis(),
	63	fFemtoParticleCut(0),
	64	fFlowParticleCut(0),
65	fVertexZBins(a.fVertexZBins),
66	fOverFlowVertexZ(0),
67	fUnderFlowVertexZ(0),
68	fMultBins(a.fMultBins) ,
69	fOverFlowMult(0),
70	fUnderFlowMult(0),
71	fRPBins(a.fRPBins),
72	fPsi(0)
73
74	{
75	fCorrFctnCollection= 0;
76	fCorrFctnCollection = new AliFemtoCorrFctnCollection;
77	fVertexZ[0] = a.fVertexZ[0];
78	fVertexZ[1] = a.fVertexZ[1];
79	fMult[0] = a.fMult[0];
80	fMult[1] = a.fMult[1];
81	if (fMixingBuffer) delete fMixingBuffer;
82	fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
83	fMultBins,fMult[0],fMult[1],
84	fRPBins,0.0,TMath::Pi());
85	// find the right event cut
86	fEventCut = a.fEventCut->Clone();
87	// find the right femto particle cut
88	fFemtoParticleCut = a.fFemtoParticleCut->Clone();
89	// find the right flow particle cut
90	fFlowParticleCut = a.fFlowParticleCut->Clone();
91	// find the right pair cut
92	fPairCut = a.fPairCut->Clone();
93
94	if ( fEventCut ) {
95	SetEventCut(fEventCut); // this will set the myAnalysis pointer inside the cut
96	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - event cut set " << endl;
97	}
98	if ( fFemtoParticleCut ) {
99	SetFirstParticleCut(fFemtoParticleCut); // this will set the myAnalysis pointer inside the cut
100	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - femto particle cut set " << endl;
101	}
102	if ( fFlowParticleCut ) {
103	SetSecondParticleCut(fFlowParticleCut); // this will set the myAnalysis pointer inside the cut
104	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - flow particle cut set " << endl;
105	}
106	if ( fPairCut ) {
107	SetPairCut(fPairCut); // this will set the myAnalysis pointer inside the cut
108	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - pair cut set " << endl;
109	}
110
111	AliFemtoCorrFctnIterator iter;
112	for (iter=a.fCorrFctnCollection->begin(); iter!=a.fCorrFctnCollection->end();iter++){
113	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - looking for correlation functions " << endl;
114	AliFemtoCorrFctn* fctn = (*iter)->Clone();
115	if (fctn) AddCorrFctn(fctn);
116	else cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - correlation function not found " << endl;
117	}
118	fNumEventsToMix = a.fNumEventsToMix;
119	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - analysis copied " << endl;
120	}
121
122	AliFemtoAnalysisAzimuthal& AliFemtoAnalysisAzimuthal::operator=(const AliFemtoAnalysisAzimuthal& a)
123	{
124	// Assignment operator
125	if (this == &a)
126	return *this;
127
128	fCorrFctnCollection= 0;
129	fCorrFctnCollection = new AliFemtoCorrFctnCollection;
130	fVertexZ[0] = a.fVertexZ[0];
131	fVertexZ[1] = a.fVertexZ[1];
132	fMult[0] = a.fMult[0];
133	fMult[1] = a.fMult[1];
134	if (fMixingBuffer) delete fMixingBuffer;
135	fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
136	fMultBins,fMult[0],fMult[1],
137	fRPBins,0.0,TMath::Pi());
138	// find the right event cut
139	fEventCut = a.fEventCut->Clone();
140	// find the right femto particle cut
141	fFemtoParticleCut = a.fFemtoParticleCut->Clone();
142	// find the right flow particle cut
143	fFlowParticleCut = a.fFlowParticleCut->Clone();
144	// find the right pair cut
145	fPairCut = a.fPairCut->Clone();
146
147	if ( fEventCut ) {
148	SetEventCut(fEventCut); // this will set the myAnalysis pointer inside the cut
149	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - event cut set " << endl;
150	}
151	if ( fFemtoParticleCut ) {
152	SetFirstParticleCut(fFemtoParticleCut); // this will set the myAnalysis pointer inside the cut
153	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - femto particle cut set " << endl;
154	}
155	if ( fFlowParticleCut ) {
156	SetSecondParticleCut(fFlowParticleCut); // this will set the myAnalysis pointer inside the cut
157	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - flow particle cut set " << endl;
158	}
159	if ( fPairCut ) {
160	SetPairCut(fPairCut); // this will set the myAnalysis pointer inside the cut
161	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - pair cut set " << endl;
162	}
163
164	AliFemtoCorrFctnIterator iter;
165	for (iter=a.fCorrFctnCollection->begin(); iter!=a.fCorrFctnCollection->end();iter++){
166	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - looking for correlation functions " << endl;
167	AliFemtoCorrFctn* fctn = (*iter)->Clone();
168	if (fctn) AddCorrFctn(fctn);
169	else cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - correlation function not found " << endl;
170	}
171	fNumEventsToMix = a.fNumEventsToMix;
172	cout << " AliFemtoAnalysisAzimuthal::AliFemtoAnalysisAzimuthal(const AliFemtoAnalysisAzimuthal& a) - analysis copied " << endl;
173
174	return *this;
175
176	}
177
178	//____________________________
179	AliFemtoAnalysisAzimuthal::~AliFemtoAnalysisAzimuthal(){
180	// now delete every PicoEvent in the EventMixingBuffer and then the Buffer itself
181	delete fPicoEventCollectionVectorHideAway;
182	}
183
184	//_________________________
185	void AliFemtoAnalysisAzimuthal::ProcessEvent(const AliFemtoEvent* hbtEvent) {
186	// Perform event processing in bins of z vertex, multiplicity and Reaction Plane angle
187	// cout << " AliFemtoAnalysisAzimuthal::ProcessEvent(const AliFemtoEvent* hbtEvent) " << endl;
188
189	//**from AliFemtoSimpleAnalysis**
190
191	fPicoEvent=0; // we will get a new pico event, if not prevent corr. fctn to access old pico event
192	fNeventsProcessed++;
193
194	// startup for EbyE
195	fFemtoParticleCut->EventBegin(hbtEvent);
196	fFlowParticleCut->EventBegin(hbtEvent);
197	fPairCut->EventBegin(hbtEvent);
198
199	for (AliFemtoCorrFctnIterator iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
200	(*iter)->EventBegin(hbtEvent);
201	}
202
203	// event cut and event cut monitor
204	bool tmpPassEvent = fEventCut->Pass(hbtEvent);
205	if (!tmpPassEvent)
206	fEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
207	if (tmpPassEvent) {
208	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
209	FillHbtParticleCollection(fFemtoParticleCut,(AliFemtoEvent*)hbtEvent,fPicoEvent->FirstParticleCollection());
210	FillHbtParticleCollection(fFlowParticleCut,(AliFemtoEvent*)hbtEvent,fPicoEvent->SecondParticleCollection());
211
212	// get right mixing buffer
213	double vertexZ = hbtEvent->PrimVertPos().z();
214	double mult = hbtEvent->UncorrectedNumberOfPrimaries();
215	TVector2 tQ = GetQVector(fPicoEvent->SecondParticleCollection());
216	double tPsi=tQ.Phi()/2.;
217	if (tPsi > TMath::Pi()) tPsi -= TMath::Pi();
218
219	fMixingBuffer = fPicoEventCollectionVectorHideAway->PicoEventCollection(vertexZ,mult,fPsi);
220	if (!fMixingBuffer) {
221	if ( vertexZ < fVertexZ[0] ) fUnderFlowVertexZ++;
222	if ( vertexZ > fVertexZ[1] ) fOverFlowVertexZ++;
223	if ( mult < fMult[0] ) fUnderFlowMult++;
224	if ( mult > fMult[1] ) fOverFlowMult++;
225	return;
226	}
227
228	//cout << "#particles in Collection: " << fPicoEvent->FirstParticleCollection()->size() << endl;
229
230	//switch which allows only using events with ParticleCollections containing a minimum number of entries
231	if (fPicoEvent->FirstParticleCollection()->size() >= fMinSizePartCollection ) {
232	fEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
233	}
234
235	//------ Make real pairs (assume identical) ------//
236	MakePairs("real", fPicoEvent->FirstParticleCollection() );
237	//cout << "AliFemtoAnalysisAzimuthal::ProcessEvent() - reals done ";
238
239	//---- Make pairs for mixed events, looping over events in mixingBuffer ----//
240	AliFemtoPicoEvent* storedEvent;
241	AliFemtoPicoEventIterator fPicoEventIter;
242	for (fPicoEventIter=MixingBuffer()->begin();fPicoEventIter!=MixingBuffer()->end();fPicoEventIter++){
243	storedEvent = *fPicoEventIter;
244	MakePairs("mixed",fPicoEvent->FirstParticleCollection(),
245	storedEvent->FirstParticleCollection() );
246	}
247	//cout << " - mixed done " << endl;
248
249	//--------- If mixing buffer is full, delete oldest event ---------//
250	if ( MixingBufferFull() ) {
251	delete MixingBuffer()->back();
252	MixingBuffer()->pop_back();
253	}
254
255	//-------- Add current event (fPicoEvent) to mixing buffer --------//
256	MixingBuffer()->push_front(fPicoEvent);
257
258	} // if ParticleCollections are big enough (mal jun2002)
259	else{
260	fEventCut->FillCutMonitor(hbtEvent, !tmpPassEvent);
261	delete fPicoEvent;
262	}
263
264	// if currentEvent is accepted by currentAnalysis cleanup for EbyE
265	fFemtoParticleCut->EventEnd(hbtEvent);
266	fFlowParticleCut->EventEnd(hbtEvent);
267	fPairCut->EventEnd(hbtEvent);
268	for (AliFemtoCorrFctnIterator iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
269	(*iter)->EventEnd(hbtEvent);
270	}
271	}
272
273	//_______________________________________________________________________________
274	void AliFemtoAnalysisAzimuthal::MakePairs(const char* typeIn, AliFemtoParticleCollection *partCollection1,
275	AliFemtoParticleCollection *partCollection2){
276	// Build pairs, check pair cuts, and call CFs' AddRealPair() or AddMixedPair() methods.
277	// If no second particle collection is specfied, make pairs within first particle collection.
278
279	string type = typeIn;
280
281	// int swpart = ((long int) partCollection1) % 2;
282	int swpart = fNeventsProcessed % 2;
283
284	AliFemtoPair* tPair = new AliFemtoPair;
285	AliFemtoCorrFctnIterator tCorrFctnIter;
286	AliFemtoParticleIterator tPartIter1, tPartIter2;
287
288	AliFemtoParticleIterator tStartOuterLoop = partCollection1->begin(); // always
289	AliFemtoParticleIterator tEndOuterLoop = partCollection1->end(); // will be one less if identical
290	AliFemtoParticleIterator tStartInnerLoop;
291	AliFemtoParticleIterator tEndInnerLoop;
292	if (partCollection2) { // Two collections:
293	tStartInnerLoop = partCollection2->begin(); // Full inner & outer loops
294	tEndInnerLoop = partCollection2->end();
295	}
296	else { // One collection:
297	tEndOuterLoop--; // Outer loop goes to next-to-last particle
298	tEndInnerLoop = partCollection1->end() ; // Inner loop goes to last particle
299	}
300	for (tPartIter1=tStartOuterLoop;tPartIter1!=tEndOuterLoop;tPartIter1++) {
301	if (!partCollection2){
302	tStartInnerLoop = tPartIter1;
303	tStartInnerLoop++;
304	}
305	tPair->SetTrack1(*tPartIter1);
306	for (tPartIter2 = tStartInnerLoop; tPartIter2!=tEndInnerLoop;tPartIter2++) {
307	tPair->SetTrack2(*tPartIter2);
308
309	//For getting the pair angle wrt EP
310	if (type == "real"){
311	TVector2 tQ, tQVector;
312	float tPsi=0, mQx=0, mQy=0, PairAngleEP=0;
313	tQ = GetQVector(partCollection1);
314
315	mQx = tQ.X() - cos(2(tPair->Track1()->FourMomentum().Phi()))(tPair->Track1()->Track()->Pt()) - cos(2(tPair->Track2()->FourMomentum().Phi()))(tPair->Track2()->Track()->Pt());
316	mQy = tQ.Y() - sin(2(tPair->Track1()->FourMomentum().Phi()))(tPair->Track1()->Track()->Pt()) - sin(2(tPair->Track2()->FourMomentum().Phi()))(tPair->Track2()->Track()->Pt());
317
318	tQVector.Set(mQx,mQy);
319
320	tPsi=tQVector.Phi()/2.;
321	if (tPsi > TMath::Pi()) tPsi -= TMath::Pi();
322
323	PairAngleEP = (tPair->EmissionAngle() - tPsi);
324	tPair->SetPairAngleEP(PairAngleEP);
325	}
326
327	if (type == "mixed"){
328	float tPsi1=0, tPsi2=0, mQx1=0, mQx2=0,mQy1=0, mQy2=0, px1=0, px2=0, py1=0, py2=0, PairAngleEP=0;
329	TVector2 tQ1, tQ2, tQVector1, tQVector2, tP;
330
331	tQ1 = GetQVector(partCollection1);
332	tQ2 = GetQVector(partCollection2);
333
334	mQx1 = tQ1.X() - cos(2(tPair->Track1()->FourMomentum().Phi()))(tPair->Track1()->Track()->Pt());
335	mQx2 = tQ2.X() - cos(2(tPair->Track2()->FourMomentum().Phi()))(tPair->Track2()->Track()->Pt());
336	mQy1 = tQ1.Y() - sin(2(tPair->Track1()->FourMomentum().Phi()))(tPair->Track1()->Track()->Pt());
337	mQy2 = tQ2.Y() - sin(2(tPair->Track2()->FourMomentum().Phi()))(tPair->Track2()->Track()->Pt());
338
339	tQVector1.Set(mQx1,mQy1);
340	tQVector2.Set(mQx2,mQy2);
341
342	tPsi1=tQVector1.Phi()/2.;
343	if (tPsi1 > TMath::Pi()) tPsi1 -= TMath::Pi();
344
345	tPsi2=tQVector2.Phi()/2.;
346	if (tPsi2 > TMath::Pi()) tPsi2 -= TMath::Pi();
347
348	px1 = (tPair->Track1()->Track()->Pt())*cos(tPair->Track1()->FourMomentum().Phi() - tPsi1);
349	px2 = (tPair->Track2()->Track()->Pt())*cos(tPair->Track2()->FourMomentum().Phi() - tPsi2);
350	py1 = (tPair->Track1()->Track()->Pt())*sin(tPair->Track1()->FourMomentum().Phi() - tPsi1);
351	py2 = (tPair->Track2()->Track()->Pt())*sin(tPair->Track2()->FourMomentum().Phi() - tPsi2);
352
353	tP.Set(px1+px2, py1+py2);
354	PairAngleEP = tP.Phi();
355
356	tPair->SetPairAngleEP(PairAngleEP);
357	}
358
359	// The following lines have to be uncommented if you want pairCutMonitors
360	// they are not in for speed reasons
361	// bool tmpPassPair = fPairCut->Pass(tPair);
362	// fPairCut->FillCutMonitor(tPair, tmpPassPair);
363	// if ( tmpPassPair )
364
365	//---- If pair passes cut, loop over CF's and add pair to real/mixed ----//
366
367	if (!partCollection2) {
368	if (swpart) {
369	tPair->SetTrack1(*tPartIter2);
370	tPair->SetTrack2(*tPartIter1);
371	swpart = 0;
372	}
373	else {
374	tPair->SetTrack1(*tPartIter1);
375	tPair->SetTrack2(*tPartIter2);
376	swpart = 1;
377	}
378	}
379
380
381	if (fPairCut->Pass(tPair)){
382	for (tCorrFctnIter=fCorrFctnCollection->begin();
383	tCorrFctnIter!=fCorrFctnCollection->end();tCorrFctnIter++){
384	AliFemtoCorrFctn* tCorrFctn = *tCorrFctnIter;
385	if(type == "real")
386	tCorrFctn->AddRealPair(tPair);
387	else if(type == "mixed")
388	tCorrFctn->AddMixedPair(tPair);
389	else
390	cout << "Problem with pair type, type = " << type.c_str() << endl;
391	}
392	}
393	} // loop over second particle
394	} // loop over first particle
395
396	delete tPair;
397	}
398
399	//_____________________________________________
400	TVector2 AliFemtoAnalysisAzimuthal::GetQVector(AliFemtoParticleCollection* particlecollection){
401
402	TVector2 mQ;
403	float mQx=0, mQy=0;
404
405	if (!particlecollection) {
406	mQ.Set(0.0, 0.0);
407	return mQ;
408	}
409
410	AliFemtoParticle* flowparticle;
411	AliFemtoParticleIterator pIter;
412	AliFemtoParticleIterator startLoop = particlecollection->begin();
413	AliFemtoParticleIterator endLoop = particlecollection->end();
414	for (pIter=startLoop;pIter!=endLoop;pIter++){
415	flowparticle = *pIter;
416	mQx += (cos(2flowparticle->FourMomentum().Phi()))(flowparticle->Track()->Pt());
417	mQy += (sin(2flowparticle->FourMomentum().Phi()))(flowparticle->Track()->Pt());
418	}
419
420	mQ.Set(mQx,mQy);
421	return mQ;
422	}
423
424	//__________________________________________________
425	double AliFemtoAnalysisAzimuthal::GetCurrentReactionPlane()
426	{
427	return fPsi;
428	}
429
430	//_________________________
431	TList* AliFemtoAnalysisAzimuthal::GetOutputList()
432	{
433	// Collect the list of output objects to be written
434
435	TList *tOutputList = new TList();
436	TList *p1Cut = fFemtoParticleCut->GetOutputList();
437
438	TListIter nextp1(p1Cut);
439	while (TObject *obj = nextp1.Next()) {
440	tOutputList->Add(obj);
441	}
442
443	TList *pairCut = fPairCut->GetOutputList();
444
445	TIter nextpair(pairCut);
446	while (TObject *obj = nextpair()) {
447	tOutputList->Add(obj);
448	}
449
450	TList *eventCut = fEventCut->GetOutputList();
451
452	TIter nextevent(eventCut);
453	while (TObject *obj = nextevent()) {
454	tOutputList->Add(obj);
455	}
456
457	AliFemtoCorrFctnIterator iter;
458	for (iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
459	TList tListCf = (iter)->GetOutputList();
460
461	TIter nextListCf(tListCf);
462	while (TObject *obj = nextListCf()) {
463	tOutputList->Add(obj);
464	}
465	}
466
467	return tOutputList;
468
469	}