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

///////////////////////////////////////////////////////////////////////////
//                                                                       //
// AliFemtoQinvCorrFctn:                                                 //
// a simple Q-invariant correlation function                             //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

#include "AliFemtoQinvCorrFctn.h"
//#include "AliFemtoHisto.h"
#include <cstdio>

#ifdef __ROOT__
ClassImp(AliFemtoQinvCorrFctn)
#endif

//____________________________
AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
  fNumerator(0),
  fDenominator(0),
  fRatio(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)";
  char tTitNum[101] = "Num";
  strncat(tTitNum,title, 100);
  fNumerator = new TH1D(tTitNum,title,nbins,QinvLo,QinvHi);
  // set up denominator
  //title = "Den Qinv (MeV/c)";
  char tTitDen[101] = "Den";
  strncat(tTitDen,title, 100);
  fDenominator = new TH1D(tTitDen,title,nbins,QinvLo,QinvHi);
  // set up ratio
  //title = "Ratio Qinv (MeV/c)";
  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);
  //fDenominator->SetDirectory(0);
  //fRatio->SetDirectory(0);

  // to enable error bar calculation...
  fNumerator->Sumw2();
  fDenominator->Sumw2();
  fRatio->Sumw2();
  fkTMonitor->Sumw2();

  fNumDEtaDPhiS->Sumw2();
  fDenDEtaDPhiS->Sumw2();

}

//____________________________
AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(const AliFemtoQinvCorrFctn& aCorrFctn) :
  AliFemtoCorrFctn(),
  fNumerator(0),
  fDenominator(0),
  fRatio(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(){
  // destructor
  delete fNumerator;
  delete fDenominator;
  delete fRatio;
  delete fkTMonitor;
  delete fNumDEtaDPhiS;
  delete fDenDEtaDPhiS;
  delete PairReader;

}
//_________________________
AliFemtoQinvCorrFctn& AliFemtoQinvCorrFctn::operator=(const AliFemtoQinvCorrFctn& aCorrFctn)
{
  // assignment operator
  if (this == &aCorrFctn)
    return *this;

  if (fNumerator) delete fNumerator;
  fNumerator = new TH1D(*aCorrFctn.fNumerator);
  if (fDenominator) delete fDenominator;
  fDenominator = new TH1D(*aCorrFctn.fDenominator);
  if (fRatio) delete fRatio;
  fRatio = new TH1D(*aCorrFctn.fRatio);
  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;
}

//_________________________
void AliFemtoQinvCorrFctn::Finish(){
  // here is where we should normalize, fit, etc...
  // we should NOT Draw() the histos (as I had done it below),
  // since we want to insulate ourselves from root at this level
  // of the code.  Do it instead at root command line with browser.
  //  fNumerator->Draw();
  //fDenominator->Draw();
  //fRatio->Draw();
  fRatio->Divide(fNumerator,fDenominator,1.0,1.0);

}

//____________________________
AliFemtoString AliFemtoQinvCorrFctn::Report(){
  // construct report
  string stemp = "Qinv Correlation Function Report:\n";
  char ctemp[100];
  snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
  stemp += ctemp;
  snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
  stemp += ctemp;
  snprintf(ctemp , 100, "Number of entries in ratio:\t%E\n",fRatio->GetEntries());
  stemp += ctemp;
  //  stemp += mCoulombWeight->Report();
  AliFemtoString returnThis = stemp;
  return returnThis;
}
//____________________________
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();
  fkTMonitor->Write();
  if (fDetaDphiscal) {
    fNumDEtaDPhiS->Write();
    fDenDEtaDPhiS->Write();
  }
  if (fPairKinematics) {
    PairReader->Write();
  }
}
//______________________________
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(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;
}
Commit	Line	Data
76ce4b5b	1	///////////////////////////////////////////////////////////////////////////
	2	// //
	3	// AliFemtoQinvCorrFctn: //
59c3fe5c	4	// a simple Q-invariant correlation function //
76ce4b5b	5	// //
	6	///////////////////////////////////////////////////////////////////////////
	7
	8	#include "AliFemtoQinvCorrFctn.h"
	9	//#include "AliFemtoHisto.h"
	10	#include <cstdio>
	11
59c3fe5c	12	#ifdef __ROOT__
76ce4b5b	13	ClassImp(AliFemtoQinvCorrFctn)
	14	#endif
	15
	16	//____________________________
	17	AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(char* title, const int& nbins, const float& QinvLo, const float& QinvHi):
	18	fNumerator(0),
	19	fDenominator(0),
	20	fRatio(0),
59c3fe5c	21	fkTMonitor(0),
59c3fe5c	22	fDetaDphiscal(kFALSE),
070bfe14	23	fPairKinematics(kFALSE),
59c3fe5c	24	fRaddedps(1.2),
59c3fe5c	25	fNumDEtaDPhiS(0),
070bfe14	26	fDenDEtaDPhiS(0),
	27	PairReader(0)// ,
	28	// fTrack1(NULL),
	29	// fTrack2(NULL)
	30
76ce4b5b	31	{
	32	// set up numerator
	33	// title = "Num Qinv (MeV/c)";
	34	char tTitNum[101] = "Num";
	35	strncat(tTitNum,title, 100);
	36	fNumerator = new TH1D(tTitNum,title,nbins,QinvLo,QinvHi);
	37	// set up denominator
	38	//title = "Den Qinv (MeV/c)";
	39	char tTitDen[101] = "Den";
	40	strncat(tTitDen,title, 100);
	41	fDenominator = new TH1D(tTitDen,title,nbins,QinvLo,QinvHi);
	42	// set up ratio
	43	//title = "Ratio Qinv (MeV/c)";
	44	char tTitRat[101] = "Rat";
	45	strncat(tTitRat,title, 100);
	46	fRatio = new TH1D(tTitRat,title,nbins,QinvLo,QinvHi);
59c3fe5c	47
76ce4b5b	48	char tTitkT[101] = "kTDep";
76ce4b5b	49	strncat(tTitkT,title, 100);
5e2038a9	50	fkTMonitor = new TH1D(tTitkT,title,250,0.0,5.0);
59c3fe5c	51
	52	char tTitNumDeDp[101] = "NumDEtaDPhiS";
	53	strncat(tTitNumDeDp,title, 100);
	54	fNumDEtaDPhiS = new TH2D(tTitNumDeDp,title,500,-0.2TMath::Pi(),0.2TMath::Pi(),500,-0.5,0.5);
	55
	56	char tTitDenDeDp[101] = "DenDEtaDPhiS";
	57	strncat(tTitDenDeDp,title, 100);
	58	fDenDEtaDPhiS = new TH2D(tTitDenDeDp,title,500,-0.2TMath::Pi(),0.2TMath::Pi(),500,-0.5,0.5);
	59
070bfe14	60	char tTitPair[101] = "Pair";
	61	strncat(tTitPair,title, 100);
	62	PairReader = new TNtuple(tTitPair,title, "px1:py1:pz1:e1:px2:py2:pz2:e2");
	63
59c3fe5c	64
76ce4b5b	65	// this next bit is unfortunately needed so that we can have many histos of same "title"
	66	// it is neccessary if we typedef TH1D to TH1d (which we do)
	67	//fNumerator->SetDirectory(0);
	68	//fDenominator->SetDirectory(0);
	69	//fRatio->SetDirectory(0);
	70
	71	// to enable error bar calculation...
	72	fNumerator->Sumw2();
	73	fDenominator->Sumw2();
	74	fRatio->Sumw2();
	75	fkTMonitor->Sumw2();
59c3fe5c	76
	77	fNumDEtaDPhiS->Sumw2();
	78	fDenDEtaDPhiS->Sumw2();
	79
76ce4b5b	80	}
	81
	82	//____________________________
	83	AliFemtoQinvCorrFctn::AliFemtoQinvCorrFctn(const AliFemtoQinvCorrFctn& aCorrFctn) :
	84	AliFemtoCorrFctn(),
	85	fNumerator(0),
	86	fDenominator(0),
	87	fRatio(0),
59c3fe5c	88	fkTMonitor(0),
59c3fe5c	89	fDetaDphiscal(kFALSE),
070bfe14	90	fPairKinematics(kFALSE),
59c3fe5c	91	fRaddedps(1.2),
59c3fe5c	92	fNumDEtaDPhiS(0),
070bfe14	93	fDenDEtaDPhiS(0),
	94	PairReader(0)// ,
	95	// fTrack1(NULL),
	96	// fTrack2(NULL)
	97
76ce4b5b	98	{
	99	// copy constructor
	100	fNumerator = new TH1D(*aCorrFctn.fNumerator);
	101	fDenominator = new TH1D(*aCorrFctn.fDenominator);
	102	fRatio = new TH1D(*aCorrFctn.fRatio);
	103	fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
59c3fe5c	104
	105	fNumDEtaDPhiS = new TH2D(*aCorrFctn.fNumDEtaDPhiS);
	106	fDenDEtaDPhiS = new TH2D(*aCorrFctn.fDenDEtaDPhiS);
	107
	108	fDetaDphiscal = aCorrFctn.fDetaDphiscal;
	109	fRaddedps = aCorrFctn.fRaddedps;
	110
070bfe14	111	fPairKinematics = aCorrFctn.fPairKinematics;
	112
	113	if (aCorrFctn.PairReader)
	114	PairReader = (TNtuple*)aCorrFctn.PairReader;
	115
76ce4b5b	116	}
	117	//____________________________
	118	AliFemtoQinvCorrFctn::~AliFemtoQinvCorrFctn(){
	119	// destructor
	120	delete fNumerator;
	121	delete fDenominator;
	122	delete fRatio;
	123	delete fkTMonitor;
59c3fe5c	124	delete fNumDEtaDPhiS;
59c3fe5c	125	delete fDenDEtaDPhiS;
070bfe14	126	delete PairReader;
59c3fe5c	127
76ce4b5b	128	}
	129	//_________________________
	130	AliFemtoQinvCorrFctn& AliFemtoQinvCorrFctn::operator=(const AliFemtoQinvCorrFctn& aCorrFctn)
	131	{
	132	// assignment operator
	133	if (this == &aCorrFctn)
	134	return *this;
	135
	136	if (fNumerator) delete fNumerator;
	137	fNumerator = new TH1D(*aCorrFctn.fNumerator);
	138	if (fDenominator) delete fDenominator;
	139	fDenominator = new TH1D(*aCorrFctn.fDenominator);
	140	if (fRatio) delete fRatio;
	141	fRatio = new TH1D(*aCorrFctn.fRatio);
	142	if (fkTMonitor) delete fkTMonitor;
	143	fkTMonitor = new TH1D(*aCorrFctn.fkTMonitor);
	144
59c3fe5c	145	if (fNumDEtaDPhiS) delete fNumDEtaDPhiS;
	146	fNumDEtaDPhiS = new TH2D(*aCorrFctn.fNumDEtaDPhiS);
	147	if (fDenDEtaDPhiS) delete fDenDEtaDPhiS;
	148	fDenDEtaDPhiS = new TH2D(*aCorrFctn.fDenDEtaDPhiS);
	149
	150	fDetaDphiscal = aCorrFctn.fDetaDphiscal;
	151	fRaddedps = aCorrFctn.fRaddedps;
	152
070bfe14	153	fPairKinematics = aCorrFctn.fPairKinematics;
	154
	155	if (aCorrFctn.PairReader)
	156	PairReader = (TNtuple*)aCorrFctn.PairReader;
	157
76ce4b5b	158	return *this;
	159	}
	160
	161	//_________________________
	162	void AliFemtoQinvCorrFctn::Finish(){
	163	// here is where we should normalize, fit, etc...
	164	// we should NOT Draw() the histos (as I had done it below),
	165	// since we want to insulate ourselves from root at this level
	166	// of the code. Do it instead at root command line with browser.
	167	// fNumerator->Draw();
	168	//fDenominator->Draw();
	169	//fRatio->Draw();
	170	fRatio->Divide(fNumerator,fDenominator,1.0,1.0);
	171
	172	}
	173
	174	//____________________________
	175	AliFemtoString AliFemtoQinvCorrFctn::Report(){
	176	// construct report
	177	string stemp = "Qinv Correlation Function Report:\n";
	178	char ctemp[100];
	179	snprintf(ctemp , 100, "Number of entries in numerator:\t%E\n",fNumerator->GetEntries());
	180	stemp += ctemp;
	181	snprintf(ctemp , 100, "Number of entries in denominator:\t%E\n",fDenominator->GetEntries());
	182	stemp += ctemp;
	183	snprintf(ctemp , 100, "Number of entries in ratio:\t%E\n",fRatio->GetEntries());
	184	stemp += ctemp;
	185	// stemp += mCoulombWeight->Report();
	186	AliFemtoString returnThis = stemp;
	187	return returnThis;
	188	}
	189	//____________________________
	190	void AliFemtoQinvCorrFctn::AddRealPair(AliFemtoPair* pair){
	191	// add true pair
	192	if (fPairCut)
	193	if (!fPairCut->Pass(pair)) return;
59c3fe5c	194
76ce4b5b	195	double tQinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs...
59c3fe5c	196
76ce4b5b	197	fNumerator->Fill(tQinv);
76ce4b5b	198	fkTMonitor->Fill(pair->KT());
59c3fe5c	199
	200
	201	//_______________________________________
	202	if (fDetaDphiscal) {
	203
	204	double phi1 = pair->Track1()->Track()->P().Phi();
	205	double phi2 = pair->Track2()->Track()->P().Phi();
	206	double chg1 = pair->Track1()->Track()->Charge();
	207	double chg2 = pair->Track2()->Track()->Charge();
	208	double ptv1 = pair->Track1()->Track()->Pt();
	209	double ptv2 = pair->Track2()->Track()->Pt();
	210	double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
	211	double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
	212
	213	AliAODInputHandler aodH = dynamic_cast<AliAODInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
0056f672	214	Int_t magsign = 0;
59c3fe5c	215
	216	if (!aodH) {
	217	//AliWarning("Could not get AODInputHandler");
	218	}
	219	else {
7d3e2025	220	AliAODEvent *fAOD; // = new AliAODEvent()
59c3fe5c	221	fAOD = aodH->GetEvent();
7d3e2025	222	magsign = fAOD->GetMagneticField();
59c3fe5c	223	}
59c3fe5c	224
0056f672	225
59c3fe5c	226	Int_t fMagSign;
	227
	228	if (magsign > 1)
	229	fMagSign = 1;
	230	else if ( magsign < 1)
	231	fMagSign = -1;
	232	else
	233	fMagSign = magsign;
	234
	235	Double_t rad = fRaddedps;
	236
	237	double afsi0b = 0.07510020733chg1fMagSign*rad/ptv1;
	238	double afsi1b = 0.07510020733chg2fMagSign*rad/ptv2;
	239	Double_t dps6 = phi2 - phi1 + TMath::ASin(afsi1b) - TMath::ASin(afsi0b);
	240	dps6 = TVector2::Phi_mpi_pi(dps6);
	241
	242	// Double_t dps = (phi1-phi2+(TMath::ASin(-0.075chg1fMagSignrad/ptv1))-(TMath::ASin(-0.075chg2fMagSignrad/ptv2)));
	243	// dps = TVector2::Phi_mpi_pi(dps);
	244	double etad = eta2 - eta1;
	245
	246	fNumDEtaDPhiS->Fill(dps6,etad);
	247	}
	248	//_______________________________________________________________
	249
76ce4b5b	250	// cout << "AliFemtoQinvCorrFctn::AddRealPair : " << pair->qInv() << " " << tQinv <<
	251	//" " << pair->track1().FourMomentum() << " " << pair->track2().FourMomentum() << endl;
	252	}
59c3fe5c	253
76ce4b5b	254	//____________________________
	255	void AliFemtoQinvCorrFctn::AddMixedPair(AliFemtoPair* pair){
	256	// add mixed (background) pair
	257	if (fPairCut)
	258	if (!fPairCut->Pass(pair)) return;
59c3fe5c	259
76ce4b5b	260	double weight = 1.0;
	261	double tQinv = fabs(pair->QInv()); // note - qInv() will be negative for identical pairs...
	262	fDenominator->Fill(tQinv,weight);
59c3fe5c	263
070bfe14	264	if (fPairKinematics) {
	265	AliFemtoParticle* fTrack1 = pair->Track1();
	266	AliFemtoParticle* fTrack2 = pair->Track2();
	267
	268	double px1 = fTrack1->FourMomentum().vect().x();
	269	double py1 = fTrack1->FourMomentum().vect().y();
	270	double pz1 = fTrack1->FourMomentum().vect().z();
	271	double e1 = fTrack1->FourMomentum().e();
	272
	273	double px2 = fTrack2->FourMomentum().vect().x();
	274	double py2 = fTrack2->FourMomentum().vect().y();
	275	double pz2 = fTrack2->FourMomentum().vect().z();
	276	double e2 = fTrack2->FourMomentum().e();
	277	PairReader->Fill(px1, py1, pz1, e1, px2, py2, pz2, e2);
	278	}
	279
59c3fe5c	280	//_______________________________________
	281	if (fDetaDphiscal) {
	282
	283	double phi1 = pair->Track1()->Track()->P().Phi();
	284	double phi2 = pair->Track2()->Track()->P().Phi();
	285	double chg1 = pair->Track1()->Track()->Charge();
	286	double chg2 = pair->Track2()->Track()->Charge();
	287	double ptv1 = pair->Track1()->Track()->Pt();
	288	double ptv2 = pair->Track2()->Track()->Pt();
	289	double eta1 = pair->Track1()->Track()->P().PseudoRapidity();
	290	double eta2 = pair->Track2()->Track()->P().PseudoRapidity();
	291
	292	AliAODInputHandler aodH = dynamic_cast<AliAODInputHandler> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
0056f672	293	Int_t magsign = 0;
0056f672	294
59c3fe5c	295
	296	if (!aodH) {
	297	//AliWarning("Could not get AODInputHandler");
	298	cout << "Could not get AODInputHandler" << endl;
	299	}
	300	else {
0056f672	301	AliAODEvent *fAOD;
59c3fe5c	302	fAOD = aodH->GetEvent();
0056f672	303	magsign = fAOD->GetMagneticField();
59c3fe5c	304	}
59c3fe5c	305
0056f672	306
59c3fe5c	307	Int_t fMagSign;
	308	if (magsign > 1)
	309	fMagSign = 1;
	310	else if ( magsign < 1)
	311	fMagSign = -1;
	312	else
	313	fMagSign = magsign;
	314
	315	Double_t rad = fRaddedps;
	316
	317	double afsi0b = 0.07510020733chg1fMagSign*rad/ptv1;
	318	double afsi1b = 0.07510020733chg2fMagSign*rad/ptv2;
	319	Double_t dps6 = phi2 - phi1 + TMath::ASin(afsi1b) - TMath::ASin(afsi0b);
	320	dps6 = TVector2::Phi_mpi_pi(dps6);
	321	double etad = eta2 - eta1;
	322
	323	fDenDEtaDPhiS->Fill(dps6,etad);
	324	}
	325	//_______________________________________________________________
	326
76ce4b5b	327	}
	328	//____________________________
	329	void AliFemtoQinvCorrFctn::Write(){
	330	// Write out neccessary objects
59c3fe5c	331	fNumerator->Write();
59c3fe5c	332	fDenominator->Write();
76ce4b5b	333	fkTMonitor->Write();
59c3fe5c	334	if (fDetaDphiscal) {
	335	fNumDEtaDPhiS->Write();
	336	fDenDEtaDPhiS->Write();
	337	}
070bfe14	338	if (fPairKinematics) {
	339	PairReader->Write();
	340	}
76ce4b5b	341	}
	342	//______________________________
	343	TList* AliFemtoQinvCorrFctn::GetOutputList()
	344	{
	345	// Prepare the list of objects to be written to the output
	346	TList *tOutputList = new TList();
	347
59c3fe5c	348	tOutputList->Add(fNumerator);
59c3fe5c	349	tOutputList->Add(fDenominator);
76ce4b5b	350	tOutputList->Add(fkTMonitor);
59c3fe5c	351	if (fDetaDphiscal) {
	352	tOutputList->Add(fNumDEtaDPhiS);
	353	tOutputList->Add(fDenDEtaDPhiS);
	354	}
070bfe14	355	if (fPairKinematics) {
	356	tOutputList->Add(PairReader);
	357	}
76ce4b5b	358	return tOutputList;
	359	}
	360
59c3fe5c	361	void AliFemtoQinvCorrFctn::CalculateDetaDphis(Bool_t dedpsc, Double_t rad) {
	362	fDetaDphiscal = dedpsc;
	363	fRaddedps = rad;
	364	}
070bfe14	365
	366	void AliFemtoQinvCorrFctn::CalculatePairKinematics(Bool_t pk) {
	367	fPairKinematics = pk;
	368	}