PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoAnalysisAzimuthalPbPb.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  *
  10  * Author: Johanna Gramling, University of Heidelberg, jgramlin@cern.ch
  11  *         Jorge Mercado, University of Heidelberg, jmercado@cern.ch
  12  *         Vera Loggins, Wayne State University, veraloggins@wayne.edu
  13  *
  14  **************************************************************************/
  15
  16 #include "AliFemtoAnalysisAzimuthalPbPb.h"
  17 #include <TMath.h>
  18 #include <string>
  19 #include <cstdio>
  20 #include "AliFemtoParticleCollection.h"
  21 #include "AliFemtoTrackCut.h"
  22 #include "AliFemtoV0Cut.h"
  23 #include "AliFemtoPairCut.h"
  24 #include "AliFemtoPairCutRadialDistanceLM.h"
  25 #include "AliVTrack.h"
  26 #include "TVector2.h"
  27 #include "TTreeFormula.h"
  28 #include "AliFemtoKinkCut.h"
  29 #include "AliEventplane.h"
  30 #include "AliAODEvent.h"
  31 #include "AliAODVertex.h"
  32 #include "AliAODTrack.h"
  33 #include "AliFemtoSimpleAnalysis.h"
  34 #include "AliFemtoPicoEventRP.h"
  35 #include "AliFemtoPicoEventCollectionVector.h"
  36 #include "AliFemtoPicoEventCollectionVectorHideAway.h"
  37
  38 #ifdef __ROOT__
  39 ClassImp(AliFemtoAnalysisAzimuthalPbPb)
  40 #endif
  41
  42 void FillHbtParticleCollection(       AliFemtoParticleCut*         partCut,
  43                                       AliFemtoEvent*               hbtEvent,
  44                                       AliFemtoPicoEventRP*         picoevent)
  45 {
  46 //      cout << "in filling" << endl;
  47       AliEventplane* evpl;
  48       evpl = picoevent->PicoEventplane();
  49       *evpl = *hbtEvent->EP();
  50 //      cout << "EP here " << evpl->GetEventplane("Q") << endl;
  51       AliFemtoTrackCut* pCut = (AliFemtoTrackCut*) partCut;
  52       AliFemtoTrack* pParticle;
  53       AliFemtoTrackIterator pIter;
  54       AliFemtoTrackIterator startLoop = hbtEvent->TrackCollection()->begin();
  55       AliFemtoTrackIterator endLoop   = hbtEvent->TrackCollection()->end();
  56       for (pIter=startLoop;pIter!=endLoop;pIter++){
  57         pParticle = *pIter;
  58         bool tmpPassParticle = pCut->Pass(pParticle);
  59         pCut->FillCutMonitor(pParticle, tmpPassParticle);
  60 //      cout << pCut->Report() << endl;
  61         if (tmpPassParticle){
  62           AliFemtoParticle* particle = new AliFemtoParticle(pParticle,pCut->Mass());
  63           picoevent->FirstParticleCollection()->push_back(particle);
  64         }
  65       }
  66 }
  67
  68 //____________________________
  69 AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(unsigned int binsVertex, double minVertex, double maxVertex,
  70                                                        unsigned int binsMult, double minMult, double maxMult, unsigned short binsRP)
  71   :
  72   fFirstParticleCut(0),
  73   fSecondParticleCut(0),
  74   fPairCutRD(0),
  75   fPicoEventRP(0),
  76   fVertexZBins(binsVertex),
  77   fOverFlowVertexZ(0),
  78   fUnderFlowVertexZ(0),
  79   fMultBins(binsMult),
  80   fOverFlowMult(0),
  81   fUnderFlowMult(0),
  82   fRPBins(binsRP),
  83   fRP(0),
  84   fphidist(0),
  85   fpairphi(0),
  86   fRPdist(0),
  87   fsubRPdist(0),
  88   frealpsi(0),
  89   fmixedpsi(0)
  90 {
  91   fCorrFctnCollection= 0;
  92   fCorrFctnCollection = new AliFemtoCorrFctnCollection;
  93   fVertexZ[0] = minVertex;
  94   fVertexZ[1] = maxVertex;
  95   fMult[0] = minMult;
  96   fMult[1] = maxMult;
  97   if (fMixingBuffer) delete fMixingBuffer;
  98   fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
  99                                                                                      fMultBins,fMult[0],fMult[1],
 100                                                                                      fRPBins,0.0,TMath::Pi());
 101
 102   fphidist = new TH1F("fphidist","fphidist; Phi Distribution",100,-TMath::Pi(),TMath::Pi());
 103   fpairphi = new TH1F("fpairphi","fpairphi; Pair Phi Distribution",100,0,TMath::TwoPi());
 104   fRPdist = new TH1F("fRPdist","fRPdist; RP Distribution",100,0,TMath::Pi());
 105   fsubRPdist = new TH1F("fsubRPdist","fsubRPdist; sub RP Distribution",200,-TMath::Pi(),TMath::Pi());
 106   frealpsi = new TH1F("frealpsi","frealpsi; real Psi Distribution",100,0,TMath::RadToDeg()*TMath::Pi());
 107   fmixedpsi = new TH1F("fmixedpsi","fmixedpsi; mixed Psi Distribution",100,0,TMath::RadToDeg()*TMath::Pi());
 108 }
 109
 110 //____________________________
 111
 112 AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) :
 113   AliFemtoSimpleAnalysis(),
 114   fFirstParticleCut(0),
 115   fSecondParticleCut(0),
 116   fPairCutRD(0),
 117   fPicoEventRP(0),
 118   fVertexZBins(a.fVertexZBins),
 119   fOverFlowVertexZ(0),
 120   fUnderFlowVertexZ(0),
 121   fMultBins(a.fMultBins) ,
 122   fOverFlowMult(0),
 123   fUnderFlowMult(0),
 124   fRPBins(a.fRPBins),
 125   fRP(0),
 126   fphidist(0),
 127   fpairphi(0),
 128   fRPdist(0),
 129   fsubRPdist(0),
 130   frealpsi(0),
 131   fmixedpsi(0)
 132
 133 {
 134   fCorrFctnCollection= 0;
 135   fCorrFctnCollection = new AliFemtoCorrFctnCollection;
 136   fVertexZ[0] = a.fVertexZ[0];
 137   fVertexZ[1] = a.fVertexZ[1];
 138   fMult[0] = a.fMult[0];
 139   fMult[1] = a.fMult[1];
 140   if (fMixingBuffer) delete fMixingBuffer;
 141   fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
 142                                                                                      fMultBins,fMult[0],fMult[1],
 143                                                                                      fRPBins,0.0,TMath::Pi());
 144   // find the right event cut
 145   fEventCut = a.fEventCut->Clone();
 146   // find the right femto particle cut
 147   fFirstParticleCut = a.fFirstParticleCut->Clone();
 148   // find the right flow particle cut
 149   fSecondParticleCut = a.fSecondParticleCut->Clone();
 150   // find the right pair cut
 151   fPairCut = a.fPairCut->Clone();
 152
 153   if ( fEventCut ) {
 154       SetEventCut(fEventCut); // this will set the myAnalysis pointer inside the cut
 155       cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - event cut set " << endl;
 156   }
 157   if ( fFirstParticleCut ) {
 158       SetFirstParticleCut(fFirstParticleCut); // this will set the myAnalysis pointer inside the cut
 159       cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - femto particle cut set " << endl;
 160   }
 161   if ( fSecondParticleCut ) {
 162       SetSecondParticleCut(fSecondParticleCut); // this will set the myAnalysis pointer inside the cut
 163       cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - flow particle cut set " << endl;
 164   }
 165   if ( fPairCut ) {
 166       SetPairCut(fPairCut); // this will set the myAnalysis pointer inside the cut
 167       cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - pair cut set " << endl;
 168   }
 169
 170   AliFemtoCorrFctnIterator iter;
 171   for (iter=a.fCorrFctnCollection->begin(); iter!=a.fCorrFctnCollection->end();iter++){
 172     cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - looking for correlation functions " << endl;
 173     AliFemtoCorrFctn* fctn = (*iter)->Clone();
 174     if (fctn) AddCorrFctn(fctn);
 175     else cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - correlation function not found " << endl;
 176   }
 177   fNumEventsToMix = a.fNumEventsToMix;
 178   cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - analysis copied " << endl;
 179 }
 180
 181 AliFemtoAnalysisAzimuthalPbPb& AliFemtoAnalysisAzimuthalPbPb::operator=(const AliFemtoAnalysisAzimuthalPbPb& a)
 182 {
 183   // Assignment operator
 184   if (this == &a)
 185     return *this;
 186
 187   fCorrFctnCollection= 0;
 188   fCorrFctnCollection = new AliFemtoCorrFctnCollection;
 189   fVertexZ[0] = a.fVertexZ[0];
 190   fVertexZ[1] = a.fVertexZ[1];
 191   fMult[0] = a.fMult[0];
 192   fMult[1] = a.fMult[1];
 193   if (fMixingBuffer) delete fMixingBuffer;
 194   fPicoEventCollectionVectorHideAway = new AliFemtoPicoEventCollectionVectorHideAway(fVertexZBins,fVertexZ[0],fVertexZ[1],
 195                                                                                      fMultBins,fMult[0],fMult[1],
 196                                                                                      fRPBins,0.0,TMath::Pi());
 197   // find the right event cut
 198   fEventCut = a.fEventCut->Clone();
 199   // find the right femto particle cut
 200   fFirstParticleCut = a.fFirstParticleCut->Clone();
 201   // find the right flow particle cut
 202   fSecondParticleCut = a.fSecondParticleCut->Clone();
 203   // find the right pair cut
 204   fPairCut = a.fPairCut->Clone();
 205
 206   if ( fEventCut ) {
 207       SetEventCut(fEventCut); // this will set the myAnalysis pointer inside the cut
 208       cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - event cut set " << endl;
 209   }
 210   if ( fFirstParticleCut ) {
 211       SetFirstParticleCut(fFirstParticleCut); // this will set the myAnalysis pointer inside the cut
 212       cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - femto particle cut set " << endl;
 213   }
 214   if ( fSecondParticleCut ) {
 215       SetSecondParticleCut(fSecondParticleCut); // this will set the myAnalysis pointer inside the cut
 216       cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - flow particle cut set " << endl;
 217   }
 218   if ( fPairCut ) {
 219       SetPairCut(fPairCut); // this will set the myAnalysis pointer inside the cut
 220       cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - pair cut set " << endl;
 221   }
 222
 223   AliFemtoCorrFctnIterator iter;
 224   for (iter=a.fCorrFctnCollection->begin(); iter!=a.fCorrFctnCollection->end();iter++){
 225     cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - looking for correlation functions " << endl;
 226     AliFemtoCorrFctn* fctn = (*iter)->Clone();
 227     if (fctn) AddCorrFctn(fctn);
 228     else cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - correlation function not found " << endl;
 229   }
 230   fNumEventsToMix = a.fNumEventsToMix;
 231   cout << " AliFemtoAnalysisAzimuthalPbPb::AliFemtoAnalysisAzimuthalPbPb(const AliFemtoAnalysisAzimuthalPbPb& a) - analysis copied " << endl;
 232
 233   return *this;
 234
 235 }
 236
 237 //____________________________
 238 AliFemtoAnalysisAzimuthalPbPb::~AliFemtoAnalysisAzimuthalPbPb(){
 239   // now delete every PicoEvent in the EventMixingBuffer and then the Buffer itself
 240   delete fPicoEventCollectionVectorHideAway;
 241 }
 242
 243 //_________________________
 244 void AliFemtoAnalysisAzimuthalPbPb::ProcessEvent(const AliFemtoEvent* hbtEvent) {
 245   // Perform event processing in bins of z vertex, multiplicity and Reaction Plane angle
 246   //****from AliFemtoSimpleAnalysis****
 247 // cout << "in processing event" << endl;
 248   fFirstParticleCut->EventBegin(hbtEvent);
 249 //   cout << "what is after that?" << endl;
 250   double vertexZ = hbtEvent->PrimVertPos().z();
 251 //   cout << "vertexZ" << vertexZ << endl;
 252   double mult = hbtEvent->UncorrectedNumberOfPrimaries();
 253 //   cout << "cent" << mult << endl;
 254   double RP = hbtEvent->ReactionPlaneAngle();
 255 //   cout << "RP " << RP << endl;
 256   fMixingBuffer = fPicoEventCollectionVectorHideAway->PicoEventCollection(vertexZ,mult,RP);
 257   if (!fMixingBuffer) {
 258 //     cout << "no mixing buffer!!!" << endl;
 259     if ( vertexZ < fVertexZ[0] ) fUnderFlowVertexZ++;
 260     if ( vertexZ > fVertexZ[1] ) fOverFlowVertexZ++;
 261     if ( mult < fMult[0] ) fUnderFlowMult++;
 262     if ( mult > fMult[1] ) fOverFlowMult++;
 263     return;
 264   }
 265
 266 //   cout << "in process event by simple analysis" << endl;
 267   // Add event to processed events
 268   fPicoEventRP=0; // we will get a new pico event, if not prevent corr. fctn to access old pico event
 269   fNeventsProcessed++;
 270
 271   fFirstParticleCut->EventBegin(hbtEvent);
 272   fSecondParticleCut->EventBegin(hbtEvent);
 273   fPairCut->EventBegin(hbtEvent);
 274   fPairCutRD->EventBegin(hbtEvent);
 275
 276   int magsign = 0;
 277   if(hbtEvent->MagneticField()>0) magsign = 1;
 278   else if(hbtEvent->MagneticField()<0) magsign = -1;
 279   fPairCutRD->SetMagneticFieldSign(magsign);
 280 //   cout << "magnetic field " << hbtEvent->MagneticField() << "magsign " << magsign << endl;
 281
 282   for (AliFemtoCorrFctnIterator iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
 283     (*iter)->EventBegin(hbtEvent);
 284   }
 285
 286   // event cut and event cut monitor
 287   bool tmpPassEvent = fEventCut->Pass(hbtEvent);
 288   if (!tmpPassEvent) {
 289     cout << "event not passed!!!!!!!!!!!!!!!!!!!!!" << endl;
 290     fEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
 291   }
 292   if (tmpPassEvent) {
 293
 294     if (RP>0){
 295      cout << "blaaaa " << hbtEvent->ReactionPlaneAngle() << endl;
 296 //       cout << " what is here?" << endl;
 297     fPicoEventRP = new AliFemtoPicoEventRP; // this is what we will make pairs from and put in Mixing Buffer
 298     // no memory leak. we will delete picoevents when they come out of the mixing buffer
 299 // cout << " what is here?" << endl;
 300     FillHbtParticleCollection(fFirstParticleCut,(AliFemtoEvent*)hbtEvent,fPicoEventRP);
 301
 302 //     cout << "after filling" << endl;
 303
 304 //     cout << "here" << endl;
 305    cout << "here " << fPicoEventRP->FirstParticleCollection()->size() << endl;
 306     if (fPicoEventRP->FirstParticleCollection()->size() >= fMinSizePartCollection) {
 307       fEventCut->FillCutMonitor(hbtEvent, tmpPassEvent);
 308 //       cout << "and here?" << endl;
 309       fRPdist->Fill(RP);
 310 //       fsubRPdist->Fill(fPicoEventRP->PicoEventplane()->GetQsubRes());
 311 //       cout << "before making real pairs" << endl;
 312
 313         MakePairs("real", fPicoEventRP);
 314
 315 //       cout << "AliFemtoSimpleAnalysis::ProcessEvent() - reals done ";
 316
 317       //---- Make pairs for mixed events, looping over events in mixingBuffer ----//
 318
 319       AliFemtoPicoEventRP* storedEvent;
 320       AliFemtoPicoEventIterator fPicoEventIter;
 321       for (fPicoEventIter=MixingBuffer()->begin();fPicoEventIter!=MixingBuffer()->end();fPicoEventIter++){
 322         storedEvent = (AliFemtoPicoEventRP*) *fPicoEventIter;
 323
 324 //      cout << "before making mixed pairs" << endl;
 325
 326           MakePairs("mixed",fPicoEventRP,
 327                             storedEvent );
 328
 329       }
 330 //       cout << " - mixed done   " << endl;
 331
 332       if ( MixingBufferFull() ) {
 333         delete MixingBuffer()->back();
 334         MixingBuffer()->pop_back();
 335       }
 336
 337       MixingBuffer()->push_front(fPicoEventRP);
 338
 339     }  // if ParticleCollections are big enough (mal jun2002)
 340     else{
 341 //       cout << "here down" << endl;
 342       fEventCut->FillCutMonitor(hbtEvent, !tmpPassEvent);
 343 //       cout << "and here?" << endl;
 344       delete fPicoEventRP;
 345 //       cout << "and here?" << endl;
 346     }
 347   }   // if currentEvent is accepted by currentAnalysis
 348 //   cout << "1" << endl;
 349   fFirstParticleCut->EventEnd(hbtEvent);
 350 //   cout << "2" << endl;
 351   fSecondParticleCut->EventEnd(hbtEvent);
 352 //   cout << "3" << endl;
 353   fPairCut->EventEnd(hbtEvent);
 354   fPairCutRD->EventEnd(hbtEvent);
 355 //   cout << "4" << endl;
 356   for (AliFemtoCorrFctnIterator iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
 357     (*iter)->EventEnd(hbtEvent);
 358   }
 359   }
 360 }
 361
 362 //_______________________________________________________________________________
 363 void AliFemtoAnalysisAzimuthalPbPb::MakePairs(const char* typeIn, AliFemtoPicoEventRP *picoevent1,
 364                                        AliFemtoPicoEventRP *picoevent2){
 365 //   cout << "In makepairs" << endl;
 366    string type = typeIn;
 367
 368   //  int swpart = ((long int) partCollection1) % 2;
 369   int swpart = fNeventsProcessed % 2;
 370
 371   AliFemtoParticleCollection* partCollection1 = picoevent1->FirstParticleCollection();
 372   AliFemtoParticleCollection* partCollection2 = 0;
 373   if (picoevent2)
 374     partCollection2 = picoevent2->FirstParticleCollection();
 375   AliFemtoPair* tPair = new AliFemtoPair;
 376
 377   AliFemtoCorrFctnIterator tCorrFctnIter;
 378
 379   AliFemtoParticleIterator tPartIter1, tPartIter2;
 380
 381   AliFemtoParticleIterator tStartOuterLoop = partCollection1->begin();  // always
 382   AliFemtoParticleIterator tEndOuterLoop   = partCollection1->end();    // will be one less if identical
 383   AliFemtoParticleIterator tStartInnerLoop;
 384   AliFemtoParticleIterator tEndInnerLoop;
 385   if (partCollection2) {                        // Two collections:
 386     tStartInnerLoop = partCollection2->begin();  //   Full inner & outer loops
 387     tEndInnerLoop   = partCollection2->end();    //
 388   }
 389   else {                                        // One collection:
 390     tEndOuterLoop--;                             //   Outer loop goes to next-to-last particle
 391     tEndInnerLoop = partCollection1->end() ;     //   Inner loop goes to last particle
 392   }
 393   for (tPartIter1=tStartOuterLoop;tPartIter1!=tEndOuterLoop;tPartIter1++) {
 394     if (!partCollection2){
 395       tStartInnerLoop = tPartIter1;
 396       tStartInnerLoop++;
 397     }
 398     tPair->SetTrack1(*tPartIter1);
 399     for (tPartIter2 = tStartInnerLoop; tPartIter2!=tEndInnerLoop;tPartIter2++) {
 400       tPair->SetTrack2(*tPartIter2);
 401
 402      if (!partCollection2) {
 403       if (swpart) {
 404           tPair->SetTrack1(*tPartIter2);
 405           tPair->SetTrack2(*tPartIter1);
 406           swpart = 0;
 407         }
 408         else {
 409           tPair->SetTrack1(*tPartIter1);
 410           tPair->SetTrack2(*tPartIter2);
 411           swpart = 1;
 412         }
 413       }
 414
 415
 416       //For getting the pair angle wrt EP
 417   if (type == "real"){
 418 //     cout << "real pairs " << endl;
 419         Double_t PairAngleEP=0;
 420         TVector2* q=0;
 421 //      float qx=0, qy=0;
 422         q = picoevent1->PicoEventplane()->GetQVector();
 423 //      cout << "real pairs 1" << endl;
 424
 425 //      cout << "track ID " << ((AliVTrack*)tPair->Track1()->Track())->GetID() << endl;
 426 //
 427 //      qx = q->X() - picoevent1->PicoEventplane()->GetQContributionX((AliVTrack*)tPair->Track1()->Track()) - picoevent1->PicoEventplane()->GetQContributionX((AliVTrack*)tPair->Track2()->Track());
 428 //      qy = q->Y() - picoevent1->PicoEventplane()->GetQContributionY((AliVTrack*)tPair->Track1()->Track()) - picoevent1->PicoEventplane()->GetQContributionY((AliVTrack*)tPair->Track2()->Track());
 429 //
 430 //      cout << "6" << endl;
 431 //      q->Set(qx,qy);
 432         float psi = q->Phi()/2;
 433 //      cout << "psi " << psi << endl;
 434         PairAngleEP = (tPair->EmissionAngle() - TMath::RadToDeg()*psi);
 435         while (PairAngleEP < 0) PairAngleEP += 180;
 436         while (PairAngleEP > 180) PairAngleEP -= 180;
 437         tPair->SetPairAngleEP(PairAngleEP);
 438         frealpsi->Fill(PairAngleEP);
 439         fphidist->Fill(tPair->Track1()->FourMomentum().Phi());
 440         fphidist->Fill(tPair->Track2()->FourMomentum().Phi());
 441         fpairphi->Fill(tPair->EmissionAngle()*TMath::DegToRad());
 442 //      cout << "tPair->EmissionAngle()*TMath::DegToRad()" << tPair->EmissionAngle()*TMath::DegToRad() << endl;
 443 //      cout << "real pairs 2" << endl;
 444   }
 445
 446   if (type == "mixed"){
 447 //     cout << "mixed pairs " << endl;
 448 //      float qx1=0, qy1=0, qx2=0, qy2=0;
 449         TVector2* q1=0;
 450         TVector2* q2=0;
 451         q1 = picoevent1->PicoEventplane()->GetQVector();
 452         q2 = picoevent2->PicoEventplane()->GetQVector();
 453         Double_t PairAngleEP=0;
 454 //      cout << "mixed pairs 1" << endl;
 455 //      qx1 = q1->X() - picoevent1->PicoEventplane()->GetQContributionX((AliVTrack*)tPair->Track1()->Track());
 456 //      qy1 = q1->Y() - picoevent1->PicoEventplane()->GetQContributionY((AliVTrack*)tPair->Track1()->Track());
 457 //      qx2 = q2->X() - picoevent2->PicoEventplane()->GetQContributionX((AliVTrack*)tPair->Track2()->Track());
 458 //      qy2 = q2->Y() - picoevent2->PicoEventplane()->GetQContributionY((AliVTrack*)tPair->Track2()->Track());
 459 //
 460 //      q1->Set(qx1,qy1);
 461 //      q2->Set(qx2,qy2);
 462
 463         float psi1 = q1->Phi()/2;
 464         float psi2 = q2->Phi()/2;
 465 //      cout << "mixed pairs 2" << endl;
 466         PairAngleEP = TMath::RadToDeg()*(TMath::ATan2(((tPair->Track1()->Track()->Pt()*TMath::Sin(tPair->Track1()->FourMomentum().Phi() - psi1))+(tPair->Track2()->Track()->Pt()*TMath::Sin(tPair->Track2()->FourMomentum().Phi() - psi2))),((tPair->Track1()->Track()->Pt()*TMath::Cos(tPair->Track1()->FourMomentum().Phi() - psi1))+(tPair->Track2()->Track()->Pt()*TMath::Cos(tPair->Track2()->FourMomentum().Phi() - psi2)))));
 467         while (PairAngleEP < 0) PairAngleEP += 180;
 468         while (PairAngleEP > 180) PairAngleEP -= 180;
 469 //      cout << "PairAngleEP " << PairAngleEP << endl;
 470         fmixedpsi->Fill(PairAngleEP);
 471         tPair->SetPairAngleEP(PairAngleEP);
 472 //      cout << "mixed pairs 3" << endl;
 473 //      double eta1 = tPair->Track1()->Track()->P().PseudoRapidity();
 474 //      double eta2 = tPair->Track2()->Track()->P().PseudoRapidity();
 475 //      cout << "mixed pair: deta = " << fabs(eta2 - eta1) << endl;
 476 //
 477 //      double phi1 = tPair->Track1()->Track()->P().Phi();
 478 //      double phi2 = tPair->Track2()->Track()->P().Phi();
 479 //      double chg1 = tPair->Track1()->Track()->Charge();
 480 //      double chg2 = tPair->Track2()->Track()->Charge();
 481 //      double ptv1 = tPair->Track1()->Track()->Pt();
 482 //      double ptv2 = tPair->Track2()->Track()->Pt();
 483 //
 484 //      Double_t dps = (phi1-phi2+(TMath::ASin(+0.075*chg1*1.2/ptv1))-(TMath::ASin(+0.075*chg2*1.2/ptv2)));
 485 //      cout << "mixed pair: dphi = " << fabs(dps) << endl;
 486   }
 487
 488       if (fPairCutRD->Pass(tPair)){
 489         for (tCorrFctnIter=fCorrFctnCollection->begin();
 490              tCorrFctnIter!=fCorrFctnCollection->end();tCorrFctnIter++){
 491           AliFemtoCorrFctn* tCorrFctn = *tCorrFctnIter;
 492           if(type == "real")
 493             tCorrFctn->AddRealPair(tPair);
 494           else if(type == "mixed") {
 495 //          cout << "mixed pair accepted" << endl;
 496             tCorrFctn->AddMixedPair(tPair);
 497           }
 498           else
 499             cout << "Problem with pair type, type = " << type.c_str() << endl;
 500         }
 501       }
 502 //       else if (!fPairCutRD->Pass(tPair)) {
 503 //      if (type == "mixed") cout << "mixed pair rejected!!!!!!!!!"<< endl;
 504 //      if (type == "real") cout << "real pair rejected!!!!!!!!!"<< endl;
 505 //       }
 506
 507
 508     }    // loop over second particle
 509
 510   }      // loop over first particle
 511
 512   delete tPair;
 513
 514 }
 515
 516 //_____________________________________________
 517 TVector2 AliFemtoAnalysisAzimuthalPbPb::GetQVector(AliFemtoParticleCollection* particlecollection){
 518
 519   TVector2 mQ;
 520   float mQx=0, mQy=0;
 521
 522   if (!particlecollection) {
 523     mQ.Set(0.0, 0.0);
 524     return mQ;
 525   }
 526
 527   AliFemtoParticle* flowparticle;
 528   AliFemtoParticleIterator pIter;
 529   AliFemtoParticleIterator startLoop = particlecollection->begin();
 530   AliFemtoParticleIterator endLoop   = particlecollection->end();
 531   for (pIter=startLoop;pIter!=endLoop;pIter++){
 532     flowparticle = *pIter;
 533     mQx += (cos(2*flowparticle->FourMomentum().Phi()))*(flowparticle->Track()->Pt());
 534     mQy += (sin(2*flowparticle->FourMomentum().Phi()))*(flowparticle->Track()->Pt());
 535   }
 536
 537   mQ.Set(mQx,mQy);
 538   return mQ;
 539 }
 540
 541 //__________________________________________________
 542 double AliFemtoAnalysisAzimuthalPbPb::GetCurrentReactionPlane()
 543 {
 544   return fRP;
 545 }
 546
 547 //______________________________________________________________________________
 548 void AliFemtoAnalysisAzimuthalPbPb::SetEPhistname(char* histname)
 549 {
 550   fphidist->SetName(Form("fphidist%s",histname));
 551   fpairphi->SetName(Form("fpairphi%s",histname));
 552   fRPdist->SetName(Form("fRPdist%s",histname));
 553   fsubRPdist->SetName(Form("fsubRPdist%s",histname));
 554   frealpsi->SetName(Form("frealpsi%s",histname));
 555   fmixedpsi->SetName(Form("fmixedpsi%s",histname));
 556 }
 557
 558 //_________________________
 559 TList* AliFemtoAnalysisAzimuthalPbPb::GetOutputList()
 560 {
 561   // Collect the list of output objects to be written
 562
 563   TList *tOutputList = new TList();
 564 //   tOutputList->SetOwner();
 565 //   TList *p1Cut = fFirstParticleCut->GetOutputList();
 566 //
 567 //   TListIter nextp1(p1Cut);
 568 //   while (TObject *obj = nextp1.Next()) {
 569 //     tOutputList->Add(obj);
 570 //   }
 571 //
 572 //   TList *pairCut = fPairCut->GetOutputList();
 573 //
 574 //   TIter nextpair(pairCut);
 575 //   while (TObject *obj = nextpair()) {
 576 //     tOutputList->Add(obj);
 577 //   }
 578 //
 579 //   TList *eventCut = fEventCut->GetOutputList();
 580 //
 581 //   TIter nextevent(eventCut);
 582 //   while (TObject *obj = nextevent()) {
 583 //     tOutputList->Add(obj);
 584 //   }
 585
 586   AliFemtoCorrFctnIterator iter;
 587   for (iter=fCorrFctnCollection->begin(); iter!=fCorrFctnCollection->end();iter++){
 588     TList *tListCf = (*iter)->GetOutputList();
 589
 590     TIter nextListCf(tListCf);
 591     while (TObject *obj = nextListCf()) {
 592       tOutputList->Add(obj);
 593     }
 594   }
 595
 596   tOutputList->Add(fphidist);
 597   tOutputList->Add(fpairphi);
 598   tOutputList->Add(fRPdist);
 599   tOutputList->Add(fsubRPdist);
 600   tOutputList->Add(frealpsi);
 601   tOutputList->Add(fmixedpsi);
 602
 603   return tOutputList;
 604
 605 }