Add Reaction Plane aware analysis
authorakisiel <akisiel@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 9 Sep 2009 13:11:42 +0000 (13:11 +0000)
committerakisiel <akisiel@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 9 Sep 2009 13:11:42 +0000 (13:11 +0000)
PWG2/FEMTOSCOPY/AliFemto/AliFemtoAnalysisReactionPlane.cxx
PWG2/FEMTOSCOPY/AliFemto/AliFemtoAnalysisReactionPlane.h
PWG2/FEMTOSCOPY/AliFemto/AliFemtoBPLCMS3DCorrFctn.cxx
PWG2/FEMTOSCOPY/AliFemto/AliFemtoBPLCMS3DCorrFctn.h
PWG2/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx
PWG2/FEMTOSCOPY/AliFemto/AliFemtoKTPairCut.cxx
PWG2/FEMTOSCOPY/AliFemto/AliFemtoKTPairCut.h
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.cxx
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelBPLCMSCorrFctn.h
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnSource.cxx
PWG2/FEMTOSCOPY/AliFemtoUser/AliFemtoModelCorrFctnSource.h

index af893de..2812e83 100644 (file)
@@ -171,7 +171,7 @@ void AliFemtoAnalysisReactionPlane::ProcessEvent(const AliFemtoEvent* hbtEvent)
   // get right mixing buffer
   double vertexZ = hbtEvent->PrimVertPos().z();
   double mult = hbtEvent->UncorrectedNumberOfPrimaries();
-  fCurrentRP = hbtEvent->GetReactionPlaneAngle();
+  fCurrentRP = hbtEvent->ReactionPlaneAngle();
   
   fMixingBuffer = fPicoEventCollectionVectorHideAway->PicoEventCollection(vertexZ,mult,fCurrentRP); 
   if (!fMixingBuffer) {
@@ -185,7 +185,7 @@ void AliFemtoAnalysisReactionPlane::ProcessEvent(const AliFemtoEvent* hbtEvent)
   AliFemtoSimpleAnalysis::ProcessEvent(hbtEvent);
 }
 
-double AliFemtoSimpleAnalysis::GetCurrentReactionPlane()
+double AliFemtoAnalysisReactionPlane::GetCurrentReactionPlane()
 {
   return fCurrentRP;
 }
index ae8cf50..94cbca9 100644 (file)
@@ -6,8 +6,8 @@
 //                                                                        //
 ////////////////////////////////////////////////////////////////////////////
 
-#ifndef ALIFEMTOVERTEXMULTANALYSIS_H
-#define ALIFEMTOVERTEXMULTANALYSIS_H
+#ifndef ALIFEMTOANALYSISREACTIONPLANE_H
+#define ALIFEMTOANALYSISREACTIONPLANE_H
 
 #include "AliFemtoSimpleAnalysis.h"        // base analysis class
 
index 8d384c6..9c4a7b8 100644 (file)
@@ -11,6 +11,8 @@
 ///////////////////////////////////////////////////////////////////////////
 
 #include "AliFemtoBPLCMS3DCorrFctn.h"
+#include "AliFemtoKTPairCut.h"
+#include "AliFemtoAnalysisReactionPlane.h"
 //#include "AliFemtoHisto.h"
 #include <cstdio>
 
@@ -41,7 +43,8 @@ AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins
   fQinvNormLo(0),
   fQinvNormHi(0),
   fNumRealsNorm(0),
-  fNumMixedNorm(0)
+  fNumMixedNorm(0),
+  fUseRPSelection(0)
 {
   // Basic constructor
   // set some stuff...
@@ -58,23 +61,23 @@ AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins
   // set up numerator
   char tTitNum[100] = "Num";
   strcat(tTitNum,title);
-  fNumerator = new TH3D(tTitNum,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fNumerator = new TH3D(tTitNum,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
   // set up denominator
   char tTitDen[100] = "Den";
   strcat(tTitDen,title);
-  fDenominator = new TH3D(tTitDen,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fDenominator = new TH3D(tTitDen,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
   // set up uncorrected denominator
   char tTitDenUncoul[100] = "DenNoCoul";
   strcat(tTitDenUncoul,title);
-  //  fUncorrectedDenominator = new TH3D(tTitDenUncoul,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  //  fUncorrectedDenominator = new TH3D(tTitDenUncoul,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
   // set up ratio
   char tTitRat[100] = "Rat";
   strcat(tTitRat,title);
-  fRatio = new TH3D(tTitRat,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fRatio = new TH3D(tTitRat,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
   // set up ave qInv
   char tTitQinv[100] = "Qinv";
   strcat(tTitQinv,title);
-  fQinvHisto = new TH3D(tTitQinv,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fQinvHisto = new TH3D(tTitQinv,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
 
   // to enable error bar calculation...
   fNumerator->Sumw2();
@@ -87,13 +90,13 @@ AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins
   // here comes the "idea" numerator and denominator and ratio...
   char tTitNumID[100] = "IDNum";
   strcat(tTitNumID,title);
-  fIDNumHisto = new TH3D(tTitNumID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fIDNumHisto = new TH3D(tTitNumID,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
   char tTitDenID[100] = "IDDen";
   strcat(tTitDenID,title);
-  fIDDenHisto = new TH3D(tTitDenID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fIDDenHisto = new TH3D(tTitDenID,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
   char tTitRatID[100] = "IDRat";
   strcat(tTitRatID,title);
-  fIDRatHisto = new TH3D(tTitRatID,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fIDRatHisto = new TH3D(tTitRatID,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
 
   fIDNumHisto->Sumw2();
   fIDDenHisto->Sumw2();
@@ -103,13 +106,13 @@ AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins
   // here comes the "smeared" numerator and denominator...
   char tTitNumSM[100] = "SMNum";
   strcat(tTitNumSM,title);
-  fSMNumHisto = new TH3D(tTitNumSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fSMNumHisto = new TH3D(tTitNumSM,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
   char tTitDenSM[100] = "SMDen";
   strcat(tTitDenSM,title);
-  fSMDenHisto = new TH3D(tTitDenSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fSMDenHisto = new TH3D(tTitDenSM,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
   char tTitRatSM[100] = "SMRat";
   strcat(tTitRatSM,title);
-  fSMRatHisto = new TH3D(tTitRatSM,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fSMRatHisto = new TH3D(tTitRatSM,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
   //
   fSMNumHisto->Sumw2();
   fSMDenHisto->Sumw2();
@@ -118,12 +121,12 @@ AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(char* title, const int& nbins
   // here comes the correction factor (which is just ratio of ideal ratio to smeared ratio)
   char tTitCorrection[100] = "CorrectionFactor";
   strcat(tTitCorrection,title);
-  fCorrectionHisto = new TH3D(tTitCorrection,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);  
+  fCorrectionHisto = new TH3D(tTitCorrection,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);  
   fCorrectionHisto->Sumw2();
   // here comes the fully corrected correlation function
   char tTitCorrCF[100] = "CorrectedCF";
   strcat(tTitCorrCF,title);
-  fCorrCFHisto = new TH3D(tTitCorrCF,title,nbins,QLo,QHi,nbins,QLo,QHi,nbins,QLo,QHi);
+  fCorrCFHisto = new TH3D(tTitCorrCF,title,nbins,-QHi,QHi,nbins,-QHi,QHi,nbins,-QHi,QHi);
   fCorrCFHisto->Sumw2();
 
   // user can (and should) override these defaults...
@@ -156,7 +159,8 @@ AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(const AliFemtoBPLCMS3DCorrFct
   fQinvNormLo(0),
   fQinvNormHi(0),
   fNumRealsNorm(0),
-  fNumMixedNorm(0)
+  fNumMixedNorm(0),
+  fUseRPSelection(0)
 {
   // Copy constructor
   fIDNumHisto = new TH3D(*aCorrFctn.fIDNumHisto);
@@ -180,6 +184,7 @@ AliFemtoBPLCMS3DCorrFctn::AliFemtoBPLCMS3DCorrFctn(const AliFemtoBPLCMS3DCorrFct
   fQinvNormHi = aCorrFctn.fQinvNormHi;
   fNumRealsNorm = aCorrFctn.fNumRealsNorm;
   fNumMixedNorm = aCorrFctn.fNumMixedNorm;
+  fUseRPSelection = aCorrFctn.fUseRPSelection;
 }
 //____________________________
 AliFemtoBPLCMS3DCorrFctn::~AliFemtoBPLCMS3DCorrFctn(){
@@ -238,7 +243,8 @@ AliFemtoBPLCMS3DCorrFctn& AliFemtoBPLCMS3DCorrFctn::operator=(const AliFemtoBPLC
   fQinvNormHi = aCorrFctn.fQinvNormHi;
   fNumRealsNorm = aCorrFctn.fNumRealsNorm;
   fNumMixedNorm = aCorrFctn.fNumMixedNorm;
-  
+  fUseRPSelection = aCorrFctn.fUseRPSelection;
+
   return *this;
 }
 
@@ -361,22 +367,57 @@ AliFemtoString AliFemtoBPLCMS3DCorrFctn::Report(){
 void AliFemtoBPLCMS3DCorrFctn::AddRealPair( AliFemtoPair* pair){
   // perform operations on real pairs
   if (fPairCut){
-    if (!(fPairCut->Pass(pair))) return;
+    if (fUseRPSelection) {
+      AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
+      if (!ktc) { 
+       cout << "RP aware cut requested, but not connected to the CF" << endl;
+       if (!(fPairCut->Pass(pair))) return;
+      }
+      else {
+       AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
+       if (!arp) {
+         cout << "RP aware cut requested, but not connected to the CF" << endl;
+         if (!(fPairCut->Pass(pair))) return;
+       }
+       if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
+      }
+    }
+    else
+      if (!(fPairCut->Pass(pair))) return;
   }
 
   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
   if ((tQinv < fQinvNormHi) && (tQinv > fQinvNormLo)) fNumRealsNorm++;
-  double qOut = fabs(pair->QOutCMS());
-  double qSide = fabs(pair->QSideCMS());
-  double qLong = fabs(pair->QLongCMS());
+  double qOut = (pair->QOutCMS());
+  double qSide = (pair->QSideCMS());
+  double qLong = (pair->QLongCMS());
 
   fNumerator->Fill(qOut,qSide,qLong);
 }
 //____________________________
 void AliFemtoBPLCMS3DCorrFctn::AddMixedPair( AliFemtoPair* pair){
   // perform operations on mixed pairs
+//   if (fPairCut){
+//     if (!(fPairCut->Pass(pair))) return;
+//   }
   if (fPairCut){
-    if (!(fPairCut->Pass(pair))) return;
+    if (fUseRPSelection) {
+      AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
+      if (!ktc) { 
+       cout << "RP aware cut requested, but not connected to the CF" << endl;
+       if (!(fPairCut->Pass(pair))) return;
+      }
+      else {
+       AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
+       if (!arp) {
+         cout << "RP aware cut requested, but not connected to the CF" << endl;
+         if (!(fPairCut->Pass(pair))) return;
+       }
+       if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
+      }
+    }
+    else
+      if (!(fPairCut->Pass(pair))) return;
   }
 
   //  double CoulombWeight = (fCorrection ? fCorrection->CoulombCorrect(pair) : 1.0);
@@ -384,9 +425,9 @@ void AliFemtoBPLCMS3DCorrFctn::AddMixedPair( AliFemtoPair* pair){
 
   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
   if ((tQinv < fQinvNormHi) && (tQinv > fQinvNormLo)) fNumMixedNorm++;
-  double qOut = fabs(pair->QOutCMS());
-  double qSide = fabs(pair->QSideCMS());
-  double qLong = fabs(pair->QLongCMS());
+  double qOut = (pair->QOutCMS());
+  double qSide = (pair->QSideCMS());
+  double qLong = (pair->QLongCMS());
 
   fDenominator->Fill(qOut,qSide,qLong,tCoulombWeight);
   //  fUncorrectedDenominator->Fill(qOut,qSide,qLong,1.0);
@@ -419,3 +460,7 @@ void AliFemtoBPLCMS3DCorrFctn::AddMixedPair( AliFemtoPair* pair){
 }
 
 
+void AliFemtoBPLCMS3DCorrFctn::SetUseRPSelection(unsigned short aRPSel)
+{
+  fUseRPSelection = aRPSel;
+}
index d1e5deb..4340fc2 100644 (file)
@@ -48,6 +48,7 @@ public:
   //  void SetCoulombCorrection(AliFemtoCoulomb* Correction);
 
   void SetSpecificPairCut(AliFemtoPairCut* aCut);
+  void SetUseRPSelection(unsigned short aRPSel);
 
   //  void SetSmearPair(AliFemtoSmearPair*);
   void SetRout(double guess);
@@ -91,8 +92,10 @@ private:
   unsigned long int fNumRealsNorm; // pairs in numerator in Qinv normalization range
   unsigned long int fNumMixedNorm; // pairs in denominator in Qinv normalization range
 
-  //  AliFemtoCoulomb* fCorrection; //!
+  unsigned short fUseRPSelection;  // The pair cut uses RP selection
 
+  //  AliFemtoCoulomb* fCorrection; //!
+  
 
 #ifdef __ROOT__
   ClassDef(AliFemtoBPLCMS3DCorrFctn, 1)
index e5635e6..2ca9b1d 100644 (file)
@@ -239,6 +239,8 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoEvent(AliFemtoEvent *tEvent)
     }
   }
 
+  tEvent->SetReactionPlaneAngle(fEvent->GetHeader()->GetQTheta(0)/2.0);
+
   Int_t *motherids;
   if (mcP) {
     motherids = new Int_t[((AliAODMCParticle *) mcP->At(mcP->GetEntries()-1))->GetLabel()];
index 41531f2..5db33ba 100644 (file)
@@ -22,6 +22,7 @@
 #include "AliFemtoKTPairCut.h"
 #include <string>
 #include <cstdio>
+#include <TMath.h>
 
 #ifdef __ROOT__
 ClassImp(AliFemtoKTPairCut)
@@ -31,7 +32,9 @@ ClassImp(AliFemtoKTPairCut)
 AliFemtoKTPairCut::AliFemtoKTPairCut():
   AliFemtoPairCut(),
   fKTMin(0),
-  fKTMax(1.0e6)
+  fKTMax(1.0e6),
+  fPhiMin(0),
+  fPhiMax(360.0)
 {
   fKTMin = 0;
    fKTMax = 1.0e6;
@@ -40,17 +43,23 @@ AliFemtoKTPairCut::AliFemtoKTPairCut():
 AliFemtoKTPairCut::AliFemtoKTPairCut(double lo, double hi) :
   AliFemtoPairCut(),
   fKTMin(lo),
-  fKTMax(hi)
+  fKTMax(hi),
+  fPhiMin(0),
+  fPhiMax(360)
 {
 }
 //__________________
 AliFemtoKTPairCut::AliFemtoKTPairCut(const AliFemtoKTPairCut& c) : 
   AliFemtoPairCut(c),
   fKTMin(0),
-  fKTMax(1.0e6)
+  fKTMax(1.0e6),
+  fPhiMin(0),
+  fPhiMax(360)
 { 
   fKTMin = c.fKTMin;
   fKTMax = c.fKTMax;
+  fPhiMin = c.fPhiMin;
+  fPhiMax = c.fPhiMax;
 }
 
 //__________________
@@ -74,6 +83,7 @@ AliFemtoString AliFemtoKTPairCut::Report(){
   // Prepare a report from the execution
   string stemp = "AliFemtoKT Pair Cut \n";  char ctemp[100];
   sprintf(ctemp,"Accept pair with kT in range %f , %f",fKTMin,fKTMax);
+  sprintf(ctemp,"Accept pair with angle in range %f , %f",fPhiMin,fPhiMax);
   stemp += ctemp;
   AliFemtoString returnThis = stemp;
   return returnThis;}
@@ -88,6 +98,10 @@ TList *AliFemtoKTPairCut::ListSettings()
   tListSetttings->AddLast(new TObjString(buf));
   snprintf(buf, 200, "AliFemtoKTPairCut.ktmin=%f", fKTMin);
   tListSetttings->AddLast(new TObjString(buf));
+  snprintf(buf, 200, "AliFemtoKTPairCut.phimax=%f", fPhiMax);
+  tListSetttings->AddLast(new TObjString(buf));
+  snprintf(buf, 200, "AliFemtoKTPairCut.phimin=%f", fPhiMin);
+  tListSetttings->AddLast(new TObjString(buf));
 
   return tListSetttings;
 }
@@ -97,3 +111,38 @@ void AliFemtoKTPairCut::SetKTRange(double ktmin, double ktmax)
   fKTMin = ktmin;
   fKTMax = ktmax;
 }
+
+void AliFemtoKTPairCut::SetPhiRange(double phimin, double phimax)
+{
+  fPhiMin = phimin;
+  fPhiMax = phimax;
+}
+
+bool AliFemtoKTPairCut::Pass(const AliFemtoPair* pair, double aRPAngle)
+{
+  if (!(Pass(pair))) return false;
+
+  //  cout << "Got pair angle RP " << pair->EmissionAngle() << "   " << aRPAngle << endl;
+
+  bool temp = true;
+  double rpangle = pair->EmissionAngle();
+  if (rpangle > 180.0) rpangle -= 180.0;
+  rpangle -= aRPAngle*180/TMath::Pi();
+  if (rpangle > 180.0) rpangle -= 180.0;
+  if (rpangle < 0.0) rpangle += 180.0;
+
+  //  cout << "Got difference " << rpangle << endl;
+
+  if (fPhiMin < 0) {
+    if ((rpangle > fPhiMax) && (rpangle < 180.0+fPhiMin)) 
+      temp = false;
+  }
+  else {
+    if ((rpangle < fPhiMin) || (rpangle > fPhiMax))
+      temp = false;
+  }
+      
+  //  if (temp) cout << "Accepted !" << endl;
+
+  return temp;
+}
index affb8d6..d6aed9a 100644 (file)
@@ -42,10 +42,14 @@ public:
   virtual TList *ListSettings();
   AliFemtoPairCut* Clone();
   void SetKTRange(double ktmin, double ktmax);
-  
+  void SetPhiRange(double phimin, double phimax);
+  bool Pass(const AliFemtoPair* pair, double aRPAngle);
+
  protected:
   Double_t fKTMin;          // Minimum allowed pair transverse momentum
   Double_t fKTMax;          // Maximum allowed pair transverse momentum 
+  Double_t fPhiMin;         // Minimum angle vs. reaction plane 
+  Double_t fPhiMax;         // Maximum angle vs. reaction plane
 
 #ifdef __ROOT__
   ClassDef(AliFemtoKTPairCut, 0)
index e5503f1..83972a0 100644 (file)
@@ -9,6 +9,8 @@
 #include "AliFemtoModelBPLCMSCorrFctn.h"
 #include "AliFemtoPair.h"
 #include "AliFemtoModelManager.h"
+#include "AliFemtoKTPairCut.h"
+#include "AliFemtoAnalysisReactionPlane.h"
 #include <cstdio>
 
 #ifdef __ROOT__ 
@@ -23,7 +25,8 @@ AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(char* title, const int&
   fNumerator3DFake(0),
   fDenominator3D(0),
   fQinvHisto(0),
-  fPairCut(0)
+  fPairCut(0),
+  fUseRPSelection(0)
 {
   // set up true numerator
   char tTitNumT[100] = "Num3DTrue";
@@ -54,7 +57,8 @@ AliFemtoModelBPLCMSCorrFctn::AliFemtoModelBPLCMSCorrFctn(const AliFemtoModelBPLC
   fNumerator3DFake(0),
   fDenominator3D(0),
   fQinvHisto(0),
-  fPairCut(0)
+  fPairCut(0),
+  fUseRPSelection(0)
 {
   // Copy constructor
   fNumerator3DTrue = new TH3D(*aCorrFctn.fNumerator3DTrue);
@@ -153,8 +157,27 @@ void AliFemtoModelBPLCMSCorrFctn::AddRealPair( AliFemtoPair* pair)
 {
   // Store a real pair in numerator
   if (fPairCut){
-    if (!(fPairCut->Pass(pair))) return;
+    if (fUseRPSelection) {
+      AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
+      if (!ktc) { 
+       cout << "RP aware cut requested, but not connected to the CF" << endl;
+       if (!(fPairCut->Pass(pair))) return;
+      }
+      else {
+       AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
+       if (!arp) {
+         cout << "RP aware cut requested, but not connected to the CF" << endl;
+         if (!(fPairCut->Pass(pair))) return;
+       }
+       if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
+      }
+    }
+    else
+      if (!(fPairCut->Pass(pair))) return;
   }
+//   if (fPairCut){
+//     if (!(fPairCut->Pass(pair))) return;
+//   }
   
   Double_t weight = fManager->GetWeight(pair);
 
@@ -169,8 +192,27 @@ void AliFemtoModelBPLCMSCorrFctn::AddRealPair( AliFemtoPair* pair)
 void AliFemtoModelBPLCMSCorrFctn::AddMixedPair( AliFemtoPair* pair){
   // store mixed pair in denominator
   if (fPairCut){
-    if (!(fPairCut->Pass(pair))) return;
+    if (fUseRPSelection) {
+      AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
+      if (!ktc) { 
+       cout << "RP aware cut requested, but not connected to the CF" << endl;
+       if (!(fPairCut->Pass(pair))) return;
+      }
+      else {
+       AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
+       if (!arp) {
+         cout << "RP aware cut requested, but not connected to the CF" << endl;
+         if (!(fPairCut->Pass(pair))) return;
+       }
+       if (!(ktc->Pass(pair, arp->GetCurrentReactionPlane()))) return;
+      }
+    }
+    else
+      if (!(fPairCut->Pass(pair))) return;
   }
+//   if (fPairCut){
+//     if (!(fPairCut->Pass(pair))) return;
+//   }
 
   Double_t weight = fManager->GetWeight(pair);
 
@@ -197,3 +239,8 @@ void AliFemtoModelBPLCMSCorrFctn::SetSpecificPairCut(AliFemtoPairCut* aCut)
 {
   fPairCut = aCut;
 }
+
+void AliFemtoModelBPLCMSCorrFctn::SetUseRPSelection(unsigned short aRPSel)
+{
+  fUseRPSelection = aRPSel;
+}
index 2cee22c..ca92456 100644 (file)
@@ -22,7 +22,8 @@ class AliFemtoModelBPLCMSCorrFctn : public AliFemtoModelCorrFctn {
     fNumerator3DFake(0),
     fDenominator3D(0),
     fQinvHisto(0),
-    fPairCut(0){}
+    fPairCut(0),
+    fUseRPSelection(0){}
   AliFemtoModelBPLCMSCorrFctn(char* title, const int& nbins, const float& QLo, const float& QHi);
   AliFemtoModelBPLCMSCorrFctn(const AliFemtoModelBPLCMSCorrFctn& aCorrFctn);
   virtual ~AliFemtoModelBPLCMSCorrFctn();
@@ -39,6 +40,7 @@ class AliFemtoModelBPLCMSCorrFctn : public AliFemtoModelCorrFctn {
   virtual TList* GetOutputList();
 
   void SetSpecificPairCut(AliFemtoPairCut* aCut);
+  void SetUseRPSelection(unsigned short aRPSel);
 
   virtual AliFemtoModelCorrFctn* Clone();
 
@@ -51,6 +53,7 @@ protected:
 
   AliFemtoPairCut* fPairCut;    //! this is a PairCut specific to THIS CorrFctn, not the Analysis
 
+  unsigned short fUseRPSelection;  // The pair cut uses RP selection
 #ifdef __ROOT__
   ClassDef(AliFemtoModelBPLCMSCorrFctn, 1)
 #endif
index 1b7ecdb..d36d43e 100644 (file)
@@ -13,6 +13,8 @@
 #include "AliFemtoModelGausLCMSFreezeOutGenerator.h"
 #include "AliFemtoModelHiddenInfo.h"
 #include "AliFemtoModelCorrFctnSource.h"
+#include "AliFemtoKTPairCut.h"
+#include "AliFemtoAnalysisReactionPlane.h"
     
 //_______________________
 AliFemtoModelCorrFctnSource::AliFemtoModelCorrFctnSource(): 
@@ -21,7 +23,10 @@ AliFemtoModelCorrFctnSource::AliFemtoModelCorrFctnSource():
   fHistRSide(0),
   fHistRLong(0),
   fHistRStar(0),
-  fHistdNdR(0)
+  fHistdNdR(0),
+  fHistNumWS(0),
+  fHistDenWS(0),
+  fUseRPSelection(0)
 {
   // default constructor
   char buf[100];
@@ -37,6 +42,11 @@ AliFemtoModelCorrFctnSource::AliFemtoModelCorrFctnSource():
   sprintf(buf, "%sdNdR", title);
   fHistdNdR = new TH1D(buf,buf,100,-50.0,50.0);
 
+  sprintf(buf, "%sNWS", title);
+  fHistNumWS = new TH2D(buf,buf,50,0.0,0.5,100,0.0,2.0);
+  sprintf(buf, "%sDWS", title);
+  fHistDenWS = new TH2D(buf,buf,50,0.0,0.5,100,0.0,2.0);
+
   fHistROut->Sumw2();
   fHistRSide->Sumw2();
   fHistRLong->Sumw2();
@@ -50,7 +60,10 @@ AliFemtoModelCorrFctnSource::AliFemtoModelCorrFctnSource(const char *title, Int_
   fHistRSide(0),
   fHistRLong(0),
   fHistRStar(0),
-  fHistdNdR(0)
+  fHistdNdR(0),
+  fHistNumWS(0),
+  fHistDenWS(0),
+  fUseRPSelection(0)
 {
   // basic constructor
   char buf[100];
@@ -65,6 +78,11 @@ AliFemtoModelCorrFctnSource::AliFemtoModelCorrFctnSource(const char *title, Int_
   sprintf(buf, "%sdNdR", title);
   fHistdNdR = new TH1D(buf,buf,100,-50.0,50.0);
 
+  sprintf(buf, "%sNWS", title);
+  fHistNumWS = new TH2D(buf,buf,50,0.0,0.5,100,0.0,2.0);
+  sprintf(buf, "%sDWS", title);
+  fHistDenWS = new TH2D(buf,buf,50,0.0,0.5,100,0.0,2.0);
+
   fHistROut->Sumw2();
   fHistRSide->Sumw2();
   fHistRLong->Sumw2();
@@ -78,7 +96,10 @@ AliFemtoModelCorrFctnSource::AliFemtoModelCorrFctnSource(const AliFemtoModelCorr
   fHistRSide(0),
   fHistRLong(0),
   fHistRStar(0),
-  fHistdNdR(0)
+  fHistdNdR(0),
+  fHistNumWS(0),
+  fHistDenWS(0),
+  fUseRPSelection(0)
 {
   // copy constructor
   fHistROut = new TH1D (*aCorrFctn.fHistROut);
@@ -86,6 +107,10 @@ AliFemtoModelCorrFctnSource::AliFemtoModelCorrFctnSource(const AliFemtoModelCorr
   fHistRLong = new TH1D(*aCorrFctn.fHistRLong);
   fHistRStar = new TH1D(*aCorrFctn.fHistRStar);
   fHistdNdR = new TH1D(*aCorrFctn.fHistdNdR);
+  fHistNumWS = new TH2D(*aCorrFctn.fHistNumWS);
+  fHistDenWS = new TH2D(*aCorrFctn.fHistDenWS);
+
+  fUseRPSelection = aCorrFctn.fUseRPSelection;
 }
 //_______________________
 AliFemtoModelCorrFctnSource::~AliFemtoModelCorrFctnSource()
@@ -96,6 +121,8 @@ AliFemtoModelCorrFctnSource::~AliFemtoModelCorrFctnSource()
   if (fHistRLong) delete fHistRLong;
   if (fHistRStar) delete fHistRStar;
   if (fHistdNdR) delete fHistdNdR;
+  if (fHistNumWS) delete fHistNumWS;
+  if (fHistDenWS) delete fHistDenWS;
   if (fNumeratorTrue) delete fNumeratorTrue;
   if (fNumeratorFake) delete fNumeratorFake;
   if (fDenominator) delete fDenominator;
@@ -122,6 +149,14 @@ AliFemtoModelCorrFctnSource& AliFemtoModelCorrFctnSource::operator=(const AliFem
   if (aCorrFctn.fHistdNdR)
     fHistdNdR = new TH1D(*aCorrFctn.fHistdNdR);
   else fHistdNdR = 0;
+  if (aCorrFctn.fHistNumWS)
+    fHistNumWS = new TH2D(*aCorrFctn.fHistNumWS);
+  else fHistNumWS = 0;
+  if (aCorrFctn.fHistDenWS)
+    fHistDenWS = new TH2D(*aCorrFctn.fHistDenWS);
+  else fHistDenWS = 0;
+
+  fUseRPSelection = aCorrFctn.fUseRPSelection;
 
   return *this;
 }
@@ -138,19 +173,72 @@ AliFemtoString AliFemtoModelCorrFctnSource::Report()
 void AliFemtoModelCorrFctnSource::AddRealPair(AliFemtoPair* aPair)
 {
   // add real (effect) pair
+//   if (fPairCut){
+//     if (!(fPairCut->Pass(aPair))) return;
+//   }
+  if (fPairCut){
+    if (fUseRPSelection) {
+      AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
+      if (!ktc) { 
+       cout << "RP aware cut requested, but not connected to the CF" << endl;
+       if (!(fPairCut->Pass(aPair))) return;
+      }
+      else {
+       AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
+       if (!arp) {
+         cout << "RP aware cut requested, but not connected to the CF" << endl;
+         if (!(fPairCut->Pass(aPair))) return;
+       }
+       if (!(ktc->Pass(aPair, arp->GetCurrentReactionPlane()))) return;
+      }
+    }
+    else
+      if (!(fPairCut->Pass(aPair))) return;
+  }
+  
   AliFemtoModelCorrFctn::AddRealPair(aPair);
+
 }
 //_______________________
 void AliFemtoModelCorrFctnSource::AddMixedPair(AliFemtoPair* aPair)
 {
   // add mixed (background) pair
+//   if (fPairCut){
+//     if (!(fPairCut->Pass(aPair))) return;
+//   }
+  if (fPairCut){
+    if (fUseRPSelection) {
+      AliFemtoKTPairCut *ktc = dynamic_cast<AliFemtoKTPairCut *>(fPairCut);
+      if (!ktc) { 
+       cout << "RP aware cut requested, but not connected to the CF" << endl;
+       if (!(fPairCut->Pass(aPair))) return;
+      }
+      else {
+       AliFemtoAnalysisReactionPlane *arp = dynamic_cast<AliFemtoAnalysisReactionPlane *> (HbtAnalysis());
+       if (!arp) {
+         cout << "RP aware cut requested, but not connected to the CF" << endl;
+         if (!(fPairCut->Pass(aPair))) return;
+       }
+       if (!(ktc->Pass(aPair, arp->GetCurrentReactionPlane()))) return;
+      }
+    }
+    else
+      if (!(fPairCut->Pass(aPair))) return;
+  }
+  
   AliFemtoModelCorrFctn::AddMixedPair(aPair);
   // save the generated positions
-  fHistROut->Fill (fManager->GetWeightGenerator()->GetRStarOut());
-  fHistRSide->Fill(fManager->GetWeightGenerator()->GetRStarSide());
-  fHistRLong->Fill(fManager->GetWeightGenerator()->GetRStarLong());
-  fHistRStar->Fill(fManager->GetWeightGenerator()->GetRStar());
-  fHistdNdR->Fill (fManager->GetWeightGenerator()->GetRStar(),1.0/(fManager->GetWeightGenerator()->GetRStar()*fManager->GetWeightGenerator()->GetRStar()));
+  if (aPair->KStar() < 0.2) {
+    fHistROut->Fill (fManager->GetWeightGenerator()->GetRStarOut());
+    fHistRSide->Fill(fManager->GetWeightGenerator()->GetRStarSide());
+    fHistRLong->Fill(fManager->GetWeightGenerator()->GetRStarLong());
+    fHistRStar->Fill(fManager->GetWeightGenerator()->GetRStar());
+    fHistdNdR->Fill (fManager->GetWeightGenerator()->GetRStar(),1.0/(fManager->GetWeightGenerator()->GetRStar()*fManager->GetWeightGenerator()->GetRStar()));
+  }
+
+  fHistDenWS->Fill(aPair->QInv(), 1.0);
+  Double_t weight = fManager->GetWeight(aPair);
+  fHistNumWS->Fill(aPair->QInv(), weight);
 }
 //_______________________
 void AliFemtoModelCorrFctnSource::Write()
@@ -161,7 +249,9 @@ void AliFemtoModelCorrFctnSource::Write()
   fHistRLong->Write();
   fHistRStar->Write();
   fHistdNdR->Write();
-  
+  fHistNumWS->Write();
+  fHistDenWS->Write();
+
   AliFemtoModelCorrFctn::Write();
 }
 //________________________
@@ -175,6 +265,8 @@ TList* AliFemtoModelCorrFctnSource::GetOutputList()
   tOutputList->Add(fHistRLong);  
   tOutputList->Add(fHistRStar);  
   tOutputList->Add(fHistdNdR);  
+  tOutputList->Add(fHistDenWS);
+  tOutputList->Add(fHistNumWS);
 
   return tOutputList;
 }
@@ -187,3 +279,7 @@ AliFemtoModelCorrFctn* AliFemtoModelCorrFctnSource::Clone()
   return tCopy;
 }
 
+void AliFemtoModelCorrFctnSource::SetUseRPSelection(unsigned short aRPSel)
+{
+  fUseRPSelection = aRPSel;
+}
index f4265cd..f5ecc71 100644 (file)
@@ -9,6 +9,7 @@
 #ifndef ALIFEMTOMODELCORRFCTNSOURCE_H
 #define ALIFEMTOMODELCORRFCTNSOURCE_H
 
+#include "TH2D.h"
 #include "AliFemtoCorrFctn.h"
 #include "AliFemtoPair.h"
 #include "AliFemtoModelManager.h"
@@ -34,6 +35,7 @@ public:
 
   virtual AliFemtoModelCorrFctn* Clone();
 
+  void SetUseRPSelection(unsigned short aRPSel);
 protected:
 
   TH1D *fHistROut;     // Distribution of Rout
@@ -41,9 +43,13 @@ protected:
   TH1D *fHistRLong;    // Distribution of Rlong
   TH1D *fHistRStar;    // Distribution of RStar
   TH1D *fHistdNdR;     // Distribution of RStar weighted by Jacobian 
+  TH2D *fHistNumWS;    // Weight spread for numerator
+  TH2D *fHistDenWS;    // Weight spread for denominator
 
 private:
 
+  unsigned short fUseRPSelection;  // The pair cut uses RP selection
+
 #ifdef __ROOT__
   ClassDef(AliFemtoModelCorrFctnSource, 1)
 #endif