]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/FEMTOSCOPY/AliFemto/AliFemtoQinvCorrFctn.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGCF / FEMTOSCOPY / AliFemto / AliFemtoQinvCorrFctn.cxx
index a7301a6dc44a1380b2b5dbb998208cfce0c98f3c..ac29d09bf449150069121b5011c891d24f5e803b 100644 (file)
@@ -1,7 +1,7 @@
 ///////////////////////////////////////////////////////////////////////////
 //                                                                       //
 // AliFemtoQinvCorrFctn:                                                 //
-// a simple Q-invariant correlation function                             // 
+// a simple Q-invariant correlation function                             //
 //                                                                       //
 ///////////////////////////////////////////////////////////////////////////
 
@@ -9,7 +9,7 @@
 //#include "AliFemtoHisto.h"
 #include <cstdio>
 
-#ifdef __ROOT__ 
+#ifdef __ROOT__
 ClassImp(AliFemtoQinvCorrFctn)
 #endif
 
@@ -18,7 +18,16 @@ AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(char* title, const int& nbins, const
   fNumerator(0),
   fDenominator(0),
   fRatio(0),
-  fkTMonitor(0)
+  fkTMonitor(0),
+  fDetaDphiscal(kFALSE),
+  fPairKinematics(kFALSE),
+  fRaddedps(1.2),
+  fNumDEtaDPhiS(0),
+  fDenDEtaDPhiS(0),
+  PairReader(0)// ,
+  // fTrack1(NULL),
+  // fTrack2(NULL)
+
 {
   // set up numerator
   //  title = "Num Qinv (MeV/c)";
@@ -35,9 +44,24 @@ AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(char* title, const int& nbins, const
   char tTitRat[101] = "Rat";
   strncat(tTitRat,title, 100);
   fRatio = new TH1D(tTitRat,title,nbins,QinvLo,QinvHi);
+
   char tTitkT[101] = "kTDep";
   strncat(tTitkT,title, 100);
   fkTMonitor = new TH1D(tTitkT,title,250,0.0,5.0);
+
+  char tTitNumDeDp[101] = "NumDEtaDPhiS";
+  strncat(tTitNumDeDp,title, 100);
+  fNumDEtaDPhiS = new TH2D(tTitNumDeDp,title,500,-0.2*TMath::Pi(),0.2*TMath::Pi(),500,-0.5,0.5);
+
+  char tTitDenDeDp[101] = "DenDEtaDPhiS";
+  strncat(tTitDenDeDp,title, 100);
+  fDenDEtaDPhiS = new TH2D(tTitDenDeDp,title,500,-0.2*TMath::Pi(),0.2*TMath::Pi(),500,-0.5,0.5);
+
+  char tTitPair[101] = "Pair";
+  strncat(tTitPair,title, 100);
+  PairReader = new TNtuple(tTitPair,title,  "px1:py1:pz1:e1:px2:py2:pz2:e2");
+
+
   // this next bit is unfortunately needed so that we can have many histos of same "title"
   // it is neccessary if we typedef TH1D to TH1d (which we do)
   //fNumerator->SetDirectory(0);
@@ -49,6 +73,10 @@ AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(char* title, const int& nbins, const
   fDenominator->Sumw2();
   fRatio->Sumw2();
   fkTMonitor->Sumw2();
+
+  fNumDEtaDPhiS->Sumw2();
+  fDenDEtaDPhiS->Sumw2();
+
 }
 
 //____________________________
@@ -57,13 +85,34 @@ AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(const AliFemtoQinvCorrFctn& aCorrFctn
   fNumerator(0),
   fDenominator(0),
   fRatio(0),
-  fkTMonitor(0)
+  fkTMonitor(0),
+  fDetaDphiscal(kFALSE),
+  fPairKinematics(kFALSE),
+  fRaddedps(1.2),
+  fNumDEtaDPhiS(0),
+  fDenDEtaDPhiS(0),
+  PairReader(0)// ,
+  // fTrack1(NULL),
+  // fTrack2(NULL)
+
 {
   // copy constructor
   fNumerator = new TH1D(*aCorrFctn.fNumerator);
   fDenominator = new TH1D(*aCorrFctn.fDenominator);
   fRatio = new TH1D(*aCorrFctn.fRatio);
   fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
+
+  fNumDEtaDPhiS = new TH2D(*aCorrFctn.fNumDEtaDPhiS);
+  fDenDEtaDPhiS = new TH2D(*aCorrFctn.fDenDEtaDPhiS);
+
+  fDetaDphiscal = aCorrFctn.fDetaDphiscal;
+  fRaddedps = aCorrFctn.fRaddedps;
+
+  fPairKinematics = aCorrFctn.fPairKinematics;
+
+  if (aCorrFctn.PairReader)
+    PairReader = (TNtuple*)aCorrFctn.PairReader;
+
 }
 //____________________________
 AliFemtoQinvCorrFctn::~AliFemtoQinvCorrFctn(){
@@ -72,6 +121,10 @@ AliFemtoQinvCorrFctn::~AliFemtoQinvCorrFctn(){
   delete fDenominator;
   delete fRatio;
   delete fkTMonitor;
+  delete fNumDEtaDPhiS;
+  delete fDenDEtaDPhiS;
+  delete PairReader;
+
 }
 //_________________________
 AliFemtoQinvCorrFctn& AliFemtoQinvCorrFctn::operator=(const AliFemtoQinvCorrFctn& aCorrFctn)
@@ -89,6 +142,19 @@ AliFemtoQinvCorrFctn& AliFemtoQinvCorrFctn::operator=(const AliFemtoQinvCorrFctn
   if (fkTMonitor) delete fkTMonitor;
   fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
 
+  if (fNumDEtaDPhiS) delete fNumDEtaDPhiS;
+  fNumDEtaDPhiS = new TH2D(*aCorrFctn.fNumDEtaDPhiS);
+  if (fDenDEtaDPhiS) delete fDenDEtaDPhiS;
+  fDenDEtaDPhiS = new TH2D(*aCorrFctn.fDenDEtaDPhiS);
+
+  fDetaDphiscal = aCorrFctn.fDetaDphiscal;
+  fRaddedps = aCorrFctn.fRaddedps;
+
+  fPairKinematics = aCorrFctn.fPairKinematics;
+
+  if (aCorrFctn.PairReader)
+    PairReader = (TNtuple*)aCorrFctn.PairReader;
+
   return *this;
 }
 
@@ -125,29 +191,153 @@ void AliFemtoQinvCorrFctn::AddRealPair(AliFemtoPair* pair){
   // add true pair
   if (fPairCut)
     if (!fPairCut->Pass(pair)) return;
-  
+
   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
+
   fNumerator->Fill(tQinv);
   fkTMonitor->Fill(pair->KT());
+
+
+//_______________________________________
+  if (fDetaDphiscal) {
+
+    double phi1 = pair->Track1()->Track()->P().Phi();
+    double phi2 = pair->Track2()->Track()->P().Phi();
+    double chg1 = pair->Track1()->Track()->Charge();
+    double chg2 = pair->Track2()->Track()->Charge();
+    double ptv1 = pair->Track1()->Track()->Pt();
+    double ptv2 = pair->Track2()->Track()->Pt();
+    double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
+    double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
+
+    AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    Int_t magsign = 0;
+
+    if (!aodH) {
+      //AliWarning("Could not get AODInputHandler");
+    }
+    else {
+      AliAODEvent *fAOD; // = new AliAODEvent()
+      fAOD = aodH->GetEvent();
+      magsign = fAOD->GetMagneticField();
+    }
+
+
+    Int_t fMagSign;
+
+    if (magsign > 1)
+      fMagSign = 1;
+    else if ( magsign < 1)
+      fMagSign = -1;
+    else
+      fMagSign = magsign;
+
+    Double_t rad = fRaddedps;
+
+    double afsi0b = 0.07510020733*chg1*fMagSign*rad/ptv1;
+    double afsi1b = 0.07510020733*chg2*fMagSign*rad/ptv2;
+    Double_t dps6 =  phi2 - phi1 + TMath::ASin(afsi1b) - TMath::ASin(afsi0b);
+    dps6 = TVector2::Phi_mpi_pi(dps6);
+
+    // Double_t dps = (phi1-phi2+(TMath::ASin(-0.075*chg1*fMagSign*rad/ptv1))-(TMath::ASin(-0.075*chg2*fMagSign*rad/ptv2)));
+    // dps = TVector2::Phi_mpi_pi(dps);
+    double etad = eta2 - eta1;
+
+    fNumDEtaDPhiS->Fill(dps6,etad);
+  }
+//_______________________________________________________________
+
   //  cout << "AliFemtoQinvCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv <<
   //" " << pair->track1().FourMomentum() << " " << pair->track2().FourMomentum() << endl;
 }
+
 //____________________________
 void AliFemtoQinvCorrFctn::AddMixedPair(AliFemtoPair* pair){
   // add mixed (background) pair
   if (fPairCut)
     if (!fPairCut->Pass(pair)) return;
-  
+
   double weight = 1.0;
   double tQinv = fabs(pair->QInv());   // note - qInv() will be negative for identical pairs...
   fDenominator->Fill(tQinv,weight);
+
+  if (fPairKinematics) {
+    AliFemtoParticle* fTrack1 = pair->Track1();
+    AliFemtoParticle* fTrack2 = pair->Track2();
+
+    double px1 = fTrack1->FourMomentum().vect().x();
+    double py1 = fTrack1->FourMomentum().vect().y();
+    double pz1 = fTrack1->FourMomentum().vect().z();
+    double e1 = fTrack1->FourMomentum().e();
+
+    double px2 = fTrack2->FourMomentum().vect().x();
+    double py2 = fTrack2->FourMomentum().vect().y();
+    double pz2 = fTrack2->FourMomentum().vect().z();
+    double e2 = fTrack2->FourMomentum().e();
+    PairReader->Fill(px1, py1, pz1, e1, px2, py2, pz2, e2);
+  }
+
+//_______________________________________
+  if (fDetaDphiscal) {
+
+    double phi1 = pair->Track1()->Track()->P().Phi();
+    double phi2 = pair->Track2()->Track()->P().Phi();
+    double chg1 = pair->Track1()->Track()->Charge();
+    double chg2 = pair->Track2()->Track()->Charge();
+    double ptv1 = pair->Track1()->Track()->Pt();
+    double ptv2 = pair->Track2()->Track()->Pt();
+    double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
+    double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
+
+    AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
+    Int_t magsign = 0;
+
+
+    if (!aodH) {
+      //AliWarning("Could not get AODInputHandler");
+      cout << "Could not get AODInputHandler" << endl;
+    }
+    else {
+      AliAODEvent *fAOD;
+      fAOD = aodH->GetEvent();
+      magsign = fAOD->GetMagneticField();
+    }
+
+
+    Int_t fMagSign;
+    if (magsign > 1)
+      fMagSign = 1;
+    else if ( magsign < 1)
+      fMagSign = -1;
+    else
+      fMagSign = magsign;
+
+    Double_t rad = fRaddedps;
+
+    double afsi0b = 0.07510020733*chg1*fMagSign*rad/ptv1;
+    double afsi1b = 0.07510020733*chg2*fMagSign*rad/ptv2;
+    Double_t dps6 =  phi2 - phi1 + TMath::ASin(afsi1b) - TMath::ASin(afsi0b);
+    dps6 = TVector2::Phi_mpi_pi(dps6);
+    double etad = eta2 - eta1;
+
+    fDenDEtaDPhiS->Fill(dps6,etad);
+  }
+//_______________________________________________________________
+
 }
 //____________________________
 void AliFemtoQinvCorrFctn::Write(){
   // Write out neccessary objects
-  fNumerator->Write(); 
-  fDenominator->Write();  
+  fNumerator->Write();
+  fDenominator->Write();
   fkTMonitor->Write();
+  if (fDetaDphiscal) {
+    fNumDEtaDPhiS->Write();
+    fDenDEtaDPhiS->Write();
+  }
+  if (fPairKinematics) {
+    PairReader->Write();
+  }
 }
 //______________________________
 TList* AliFemtoQinvCorrFctn::GetOutputList()
@@ -155,11 +345,24 @@ TList* AliFemtoQinvCorrFctn::GetOutputList()
   // Prepare the list of objects to be written to the output
   TList *tOutputList = new TList();
 
-  tOutputList->Add(fNumerator); 
-  tOutputList->Add(fDenominator);  
+  tOutputList->Add(fNumerator);
+  tOutputList->Add(fDenominator);
   tOutputList->Add(fkTMonitor);
-
+  if (fDetaDphiscal) {
+    tOutputList->Add(fNumDEtaDPhiS);
+    tOutputList->Add(fDenDEtaDPhiS);
+  }
+  if (fPairKinematics) {
+    tOutputList->Add(PairReader);
+  }
   return tOutputList;
 }
 
+void AliFemtoQinvCorrFctn::CalculateDetaDphis(Bool_t dedpsc, Double_t rad) {
+  fDetaDphiscal = dedpsc;
+  fRaddedps = rad;
+}
 
+void AliFemtoQinvCorrFctn::CalculatePairKinematics(Bool_t pk) {
+  fPairKinematics = pk;
+}