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

////////////////////////////////////////////////////////////////////////////////
//345678901234567890123456789012345678901234567890123456789012345678901234567890
//       1         2         3         4         5         6         7         8
//
// Tool to study two-track effects in ALICE for femtoscopic analyses
// J. Mercado <mercado@physi.uni-heidelberg.de> Last modified: 20.01.2011
//
////////////////////////////////////////////////////////////////////////////////

#include "AliTwoTrackRes.h"
#include <iostream>
#include <fstream>
#include "TMath.h"
#include "TROOT.h"
#include "TFile.h"
#include "TChain.h"
#include "TLeaf.h"
#include "TNtuple.h"
#include "TRandom2.h"

ClassImp(AliTwoTrackRes)

//______________________________________________________________________________
// Constructor(s)

AliTwoTrackRes::AliTwoTrackRes(const char *name) : 
  AliAnalysisTask(name,""), fChain(0), fESDEvent(0), 
  fOutContainer(0), fTrackCuts(0), fNTuple1(0), 
  fNTuple2(0), fP1(), fP2(), fPb1(), fPb2(), fP(), fQ(), fTpcEnt1(), fTpcEnt2(), 
  fTpcDist(), fOutFilename()
{
  DefineInput(0, TChain::Class());     // Slot input 0 reads from a TChain  
  DefineOutput(0, TObjArray::Class()); // Slot output 0 writes into a TObjArray
}

AliTwoTrackRes::AliTwoTrackRes(const AliTwoTrackRes& aTwoTrackRes) : 
  AliAnalysisTask(aTwoTrackRes), fChain(0), fESDEvent(0), fOutContainer(0), 
  fTrackCuts(0), fNTuple1(0), fNTuple2(0), fP1(), fP2(), fPb1(), fPb2(), fP(), 
  fQ(), fTpcEnt1(), fTpcEnt2(), fTpcDist(), fOutFilename()
{
  //Copy constructor
  fChain = aTwoTrackRes.fChain;
  fESDEvent = aTwoTrackRes.fESDEvent;
  fOutContainer = aTwoTrackRes.fOutContainer;
  fTrackCuts = aTwoTrackRes.fTrackCuts;
  fNTuple1 = aTwoTrackRes.fNTuple1;
  fNTuple2 = aTwoTrackRes.fNTuple2;
  fP1 = aTwoTrackRes.fP1;
  fP2 = aTwoTrackRes.fP2;
  fPb1 = aTwoTrackRes.fPb1;
  fPb2 = aTwoTrackRes.fPb2;
  fP = aTwoTrackRes.fP;
  fQ = aTwoTrackRes.fQ;
  fTpcEnt1 = aTwoTrackRes.fTpcEnt1;
  fTpcEnt2 = aTwoTrackRes.fTpcEnt2;
  fTpcDist = aTwoTrackRes.fTpcDist;
  fOutFilename = aTwoTrackRes.fOutFilename;
}

AliTwoTrackRes& AliTwoTrackRes::operator=(const AliTwoTrackRes& aTwoTrackRes)
{
  // Assignment operator
  if (this == &aTwoTrackRes)
    return *this;
  fChain = aTwoTrackRes.fChain;
  fESDEvent = aTwoTrackRes.fESDEvent;
  fOutContainer = aTwoTrackRes.fOutContainer;
  fTrackCuts = aTwoTrackRes.fTrackCuts;
  fNTuple1 = aTwoTrackRes.fNTuple1;
  fNTuple2 = aTwoTrackRes.fNTuple2;
  fP1 = aTwoTrackRes.fP1;
  fP2 = aTwoTrackRes.fP2;
  fPb1 = aTwoTrackRes.fPb1;
  fPb2 = aTwoTrackRes.fPb2;
  fP = aTwoTrackRes.fP;
  fQ = aTwoTrackRes.fQ;
  fTpcEnt1 = aTwoTrackRes.fTpcEnt1;
  fTpcEnt2 = aTwoTrackRes.fTpcEnt2;
  fTpcDist = aTwoTrackRes.fTpcDist;
  fOutFilename = aTwoTrackRes.fOutFilename;
  return *this;
}

AliTwoTrackRes::~AliTwoTrackRes() {printf("AliTwoTrackRes destroyed\n");}

void AliTwoTrackRes::ConnectInputData(Option_t *) {  
//______________________________________________________________________________
// Connect input data and initialize track cuts

  fChain = (TChain*)GetInputData(0);
  fESDEvent = new AliESDEvent();
  fESDEvent->ReadFromTree(fChain);

  // Cuts to select primary tracks (ITS+TPC)
  fTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
  Double_t cov1, cov2, cov3, cov4, cov5; // diagonal cov. matrix elements
  Double_t nSigma;                       // max. DCA to primary vertex
  Int_t minNClustersTPC;                 // min. number of clusters per TPC tracks
  Double_t maxChi2PerClusterTPC;         // max. chi2 per cluster per TPC track  
  Int_t cutMode = 1;                     // select cut mode
  if (cutMode == 1) {
  cov1 = 2; cov2 = 2; cov3 = 0.5; cov4 = 0.5; cov5 = 2;
  nSigma = 3; minNClustersTPC = 75; maxChi2PerClusterTPC = 3.5;
  fTrackCuts->SetMaxCovDiagonalElements(cov1, cov2, cov3, cov4, cov5);
  fTrackCuts->SetMaxNsigmaToVertex(nSigma);
  fTrackCuts->SetRequireSigmaToVertex(kTRUE);
  fTrackCuts->SetRequireTPCRefit(kTRUE);
  fTrackCuts->SetAcceptKinkDaughters(kFALSE);
  fTrackCuts->SetMinNClustersTPC(minNClustersTPC);
  fTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);  
  TString tag("Global tracking");}
}

void AliTwoTrackRes::CreateOutputObjects() {
//______________________________________________________________________________
// Create output objects 

  fNTuple1 = new TNtuple("nt1","True pairs",
  "pt1:eta1:phi1:nsh1:pt2:eta2:phi2:nsh2:qinv:mindist:dist:corr:qfac");
  fNTuple2 = new TNtuple("nt2","Mixed pairs",
  "pt1:eta1:phi1:nsh1:pt2:eta2:phi2:nsh2:qinv:mindist:dist:corr:qfac");
  Int_t c = 0;
  fOutContainer = new TObjArray(2);
  fOutContainer->AddAt(fNTuple1, c++);
  fOutContainer->AddAt(fNTuple2, c++);
}

void AliTwoTrackRes::Exec(Option_t *) {
//______________________________________________________________________________
// Create true and mixed pairs keeping some track parameters

  double bfield = 5.0;
  static int nr=0;
  if (nr == 0) printf("\tStarting event loop...\n");
  printf("\rProcessing event %8d", nr);
  Double_t mpi = 0.13957; // [GeV/c^2]
  Double_t pidTrk1[AliPID::kSPECIES], pidTrk2[AliPID::kSPECIES];
  Int_t tpcIn = 80;   // [cm]
  Int_t tpcOut = 250; // [cm]
  Double_t tdist[170];
  Double_t tdistrot[170];
  Double_t tpcEnt1[3], tpcEnt2[3], pos1[3];
  TVector3 x1, x2, diff;
  TBits clu1, clu2, sha1, sha2;
  TRandom2 rnd;
  Int_t  ntracks = fESDEvent->GetNumberOfTracks();
  for(Int_t itrack = 0; itrack < ntracks; itrack++) {
    AliESDtrack *track1 = fESDEvent->GetTrack(itrack);
    AliExternalTrackParam *trp1 = const_cast<AliExternalTrackParam*> 
      (track1->GetTPCInnerParam());
    if (!trp1) continue;
    if (!track1->IsOn(AliESDtrack::kTPCpid)) continue;
    track1->GetTPCpid(pidTrk1);
    Int_t q1 = trp1->Charge();
    if (!((fTrackCuts->AcceptTrack(track1)) && (q1 == 1) &&
	  (pidTrk1[AliPID::kPion]+pidTrk1[AliPID::kMuon] > 0.5))) continue;
    if (!track1->GetInnerXYZ(tpcEnt1)) continue;
    clu1 = track1->GetTPCClusterMap();
    sha1 = track1->GetTPCSharedMap();
    SetTr1(track1->Pt(), track1->Eta(), track1->Phi(), mpi);
    SetTpcEnt1(tpcEnt1[0], tpcEnt1[1], tpcEnt1[2]);
    for(Int_t jtrack = 0; jtrack < itrack; jtrack++) {
      AliESDtrack *track2 = fESDEvent->GetTrack(jtrack);
      AliExternalTrackParam *trp2 = const_cast<AliExternalTrackParam*> 
	(track2->GetTPCInnerParam());
      if (!trp2) continue;
      if (!track2->IsOn(AliESDtrack::kTPCpid)) continue;
      track2->GetTPCpid(pidTrk2);
      Int_t q2 = trp2->Charge();
      if (!((fTrackCuts->AcceptTrack(track2)) && (q2 == 1) &&
	    (pidTrk2[AliPID::kPion]+pidTrk2[AliPID::kMuon] > 0.5))) continue;
      if (!track2->GetInnerXYZ(tpcEnt2)) continue;
      clu2 = track2->GetTPCClusterMap();
      sha2 = track2->GetTPCSharedMap();
      SetTr2(track2->Pt(), track2->Eta(), track2->Phi(), mpi);
      SetTpcEnt2(tpcEnt2[0], tpcEnt2[1], tpcEnt2[2]);
      for (Int_t i = tpcIn; i < tpcOut; i++) { // Minimum distance 
	trp1->GetDistance(trp2, (double) i, pos1, bfield);
	x1.SetXYZ(pos1[0], pos1[1], pos1[2]);
	tdist[i-tpcIn] = x1.Mag();
	x1.SetXYZ(-pos1[0], -pos1[1], pos1[2]);
	tdistrot[i-tpcIn] = x1.Mag();
      }
      Double_t mindist = 100000;
      Int_t jmin=0;
      for (Int_t j = 0; j < tpcOut-tpcIn; j++) {
	if (tdist[j] < mindist) {jmin=j;  mindist = tdist[j]; }
      }
      //      Double_t mindist = MinDist(track1, track2);
      Double_t dist = Dist();
      Double_t dphi = DPhi();
      Double_t deta = DEta();
      Int_t    nsh1 = GetNSha(clu1, sha1);
      Int_t    nsh2 = GetNSha(clu2, sha2);
      Double_t corr = Corr(clu1, clu2, sha1, sha2);
      Double_t qfac = Qfac(clu1, clu2, sha1, sha2);
      if ((TMath::Abs(track1->Eta())>0.8)&&(TMath::Abs(track2->Eta())>0.8)) continue;
      if ((TMath::Abs(dphi)<0.35)&&(deta<0.35)) {
      FillNTuple1(mindist,dist,corr,qfac,nsh1,nsh2);}    // True
      Double_t tr2rot = RotTr2Phi();                // Rotate trck2
      SetTr2(track2->Pt(), track2->Eta(), tr2rot, mpi);
      tpcEnt2[0] = -tpcEnt2[0];
      tpcEnt2[1] = -tpcEnt2[1];
      Double_t distrot = Dist();
      Double_t dphirot = DPhi();
      Double_t mindistrot = 100000;
      jmin=0;
      for (Int_t j = 0; j < tpcOut-tpcIn; j++) {
	if (tdistrot[j] < mindistrot) {jmin=j;  mindistrot = tdistrot[j]; }
      }
      if ((TMath::Abs(dphirot)<0.35)&&(deta<0.35)) {
	if (rnd.Rndm() < 0.5) NoSwap();
	else Swap(); 
	FillNTuple2(mindistrot,distrot,corr,qfac,nsh1,nsh2);} // Mixed
    }
  }
  PostData(0, fOutContainer);
  nr++;
}

void AliTwoTrackRes::Terminate(Option_t *) {
//______________________________________________________________________________
// Write output and clean up

  fOutContainer = (TObjArray*)GetOutputData(0);
  TFile *f1  = new TFile( fOutFilename, "RECREATE" );
  fOutContainer->Write();
  f1->Flush();
  f1->Close();
  delete f1;
  delete fChain;
  delete fNTuple1;
  delete fNTuple2;
  printf("\n");
}

//______________________________________________________________________________
// Miscellaneous methods

// Set tracks 
void AliTwoTrackRes::SetTr1(double pt1, double eta1, double phi1, double m) {
  fP1.SetPtEtaPhiM(pt1, eta1, phi1, m);}
void AliTwoTrackRes::SetTr2(double pt2, double eta2, double phi2, double m) {
  fP2.SetPtEtaPhiM(pt2, eta2, phi2, m);}

// Set nominal TPC entrance coordinates
void AliTwoTrackRes::SetTpcEnt1(double x1, double y1, double z1) {
  fTpcEnt1.SetX(x1); fTpcEnt1.SetY(y1); fTpcEnt1.SetZ(z1);}
void AliTwoTrackRes::SetTpcEnt2(double x2, double y2, double z2) {
  fTpcEnt2.SetX(x2); fTpcEnt2.SetY(y2); fTpcEnt2.SetZ(z2);}

double AliTwoTrackRes::MinDist(AliExternalTrackParam *trk1,
			       AliExternalTrackParam *trk2) {
// Calculate minimum track separation within the TPC

  int tpcIn = 0;   // [cm]
  int tpcOut = 170; // [cm]
  double tdist[170], pos[3];
  TVector3 x;
  for (int i = tpcIn; i < tpcOut; i++) {
    trk1->GetDistance(trk2, i, pos, 5000);
    x.SetXYZ(pos[0], pos[1], pos[2]);
    tdist[i-tpcIn] = x.Mag();
  }
  double maxdist = 0.0;
  for (int j = 0; j < tpcOut-tpcIn; j++) {
    if (tdist[j] > maxdist) { maxdist = tdist[j]; }
  }
  double mindist = maxdist;
  for (int j = 0; j < tpcOut-tpcIn; j++) {
    if (tdist[j] < mindist) { mindist = tdist[j]; }
  }
  return mindist;}

int AliTwoTrackRes::GetNSha(TBits cl, TBits sh) {
// Get number of shared clusters

  int ncl = cl.GetNbits();
  int sum = 0; 
  for(int i = 0; i < ncl; i++) {
    if (!cl.TestBitNumber(i)) continue;
    int n = sh.TestBitNumber(i);
    sum += n;}
  return sum;}

double AliTwoTrackRes::Corr(TBits cl1,  TBits cl2, TBits sh1, TBits sh2) {
// Calculate correlation coefficient

  int ncl1 = cl1.GetNbits();
  int ncl2 = cl2.GetNbits();
  double sumN = 0;  double sumX = 0;  double sumY = 0;
  double sumXX = 0; double sumYY = 0; double sumXY = 0; double corr = -2.0;
  for(int i = 0; i < ncl1 && i < ncl2; i++) {
    if (!(cl1.TestBitNumber(i)&&cl2.TestBitNumber(i))) continue;
    int x = sh1.TestBitNumber(i);
    int y = sh2.TestBitNumber(i);
    sumN += 1.0;
    sumX += x;
    sumY += y;
    sumXX += x*x;
    sumYY += y*y;
    sumXY += x*y;
  }
  double meanX = sumX/sumN;
  double meanY = sumY/sumN;
  double meanXX = sumXX/sumN;
  double meanYY = sumYY/sumN;
  double meanXY = sumXY/sumN;
  double sX = TMath::Sqrt(TMath::Abs(meanXX-meanX*meanX));
  double sY = TMath::Sqrt(TMath::Abs(meanYY-meanY*meanY));
  if (sX*sY!=0) corr = (meanXY-meanX*meanY)/(sX*sY);
  return corr;}

double AliTwoTrackRes::Qfac(TBits cl1,  TBits cl2, TBits sh1, TBits sh2) {
// Quality factor from AliFemto 

  int ncl1 = cl1.GetNbits();
  int ncl2 = cl2.GetNbits();
  int sumCls = 0; int sumSha = 0; int sumQ = 0;
  double shfrac = 0; double qfactor = 0;
  for(int i = 0; i < ncl1 && i < ncl2; i++) {
    if (cl1.TestBitNumber(i) && cl2.TestBitNumber(i)) { // Both clusters
      if (sh1.TestBitNumber(i) && sh2.TestBitNumber(i)) { // Shared
	sumQ++;
	sumCls+=2;
	sumSha+=2;}
      else {sumQ--; sumCls+=2;}
    }
    else if (cl1.TestBitNumber(i) || cl2.TestBitNumber(i)) { // Non shared
      sumQ++;
      sumCls++;}
  }
  if (sumCls>0) {
    qfactor = sumQ*1.0/sumCls;
    shfrac = sumSha*1.0/sumCls;
  }
  return qfactor;
}

// Rotate second track for mixed pairs
double AliTwoTrackRes::RotTr2Phi() {
  double rot = TVector2::Phi_mpi_pi(fP2.Phi()+TMath::Pi());
  fTpcEnt2.SetPhi(TVector2::Phi_mpi_pi(fTpcEnt2.Phi()+TMath::Pi())); 
  return rot;}

// Fill NTuples
void AliTwoTrackRes::FillNTuple1(double minsep, double sep, double corr, 
				 double qf, int ns1, int ns2) {
  fNTuple1->Fill(fP1.Pt(),fP1.Eta(),fP1.Phi(),ns1,fP2.Pt(),fP2.Eta(),
		 fP2.Phi(),ns2,Qinv(),minsep,sep,corr,qf);}
void AliTwoTrackRes::FillNTuple2(double minsep, double sep, double corr, 
				 double qf, int ns1, int ns2) {
  fNTuple2->Fill(fPb1.Pt(),fPb1.Eta(),fPb1.Phi(),ns1,fPb2.Pt(),fPb2.Eta(),
		 fPb2.Phi(),ns2,Qinv(),minsep,sep,corr,qf);}

//______________________________________________________________________________
// EOF
Commit	Line	Data
76ce4b5b	1	////////////////////////////////////////////////////////////////////////////////
	2	//345678901234567890123456789012345678901234567890123456789012345678901234567890
	3	// 1 2 3 4 5 6 7 8
	4	//
	5	// Tool to study two-track effects in ALICE for femtoscopic analyses
	6	// J. Mercado <mercado@physi.uni-heidelberg.de> Last modified: 20.01.2011
	7	//
	8	////////////////////////////////////////////////////////////////////////////////
	9
	10	#include "AliTwoTrackRes.h"
	11	#include <iostream>
	12	#include <fstream>
	13	#include "TMath.h"
	14	#include "TROOT.h"
	15	#include "TFile.h"
	16	#include "TChain.h"
	17	#include "TLeaf.h"
	18	#include "TNtuple.h"
	19	#include "TRandom2.h"
	20
	21	ClassImp(AliTwoTrackRes)
	22
	23	//______________________________________________________________________________
	24	// Constructor(s)
	25
	26	AliTwoTrackRes::AliTwoTrackRes(const char *name) :
	27	AliAnalysisTask(name,""), fChain(0), fESDEvent(0),
	28	fOutContainer(0), fTrackCuts(0), fNTuple1(0),
	29	fNTuple2(0), fP1(), fP2(), fPb1(), fPb2(), fP(), fQ(), fTpcEnt1(), fTpcEnt2(),
	30	fTpcDist(), fOutFilename()
	31	{
	32	DefineInput(0, TChain::Class()); // Slot input 0 reads from a TChain
	33	DefineOutput(0, TObjArray::Class()); // Slot output 0 writes into a TObjArray
	34	}
	35
	36	AliTwoTrackRes::AliTwoTrackRes(const AliTwoTrackRes& aTwoTrackRes) :
	37	AliAnalysisTask(aTwoTrackRes), fChain(0), fESDEvent(0), fOutContainer(0),
	38	fTrackCuts(0), fNTuple1(0), fNTuple2(0), fP1(), fP2(), fPb1(), fPb2(), fP(),
	39	fQ(), fTpcEnt1(), fTpcEnt2(), fTpcDist(), fOutFilename()
	40	{
	41	//Copy constructor
	42	fChain = aTwoTrackRes.fChain;
	43	fESDEvent = aTwoTrackRes.fESDEvent;
	44	fOutContainer = aTwoTrackRes.fOutContainer;
	45	fTrackCuts = aTwoTrackRes.fTrackCuts;
	46	fNTuple1 = aTwoTrackRes.fNTuple1;
	47	fNTuple2 = aTwoTrackRes.fNTuple2;
	48	fP1 = aTwoTrackRes.fP1;
	49	fP2 = aTwoTrackRes.fP2;
	50	fPb1 = aTwoTrackRes.fPb1;
	51	fPb2 = aTwoTrackRes.fPb2;
	52	fP = aTwoTrackRes.fP;
	53	fQ = aTwoTrackRes.fQ;
	54	fTpcEnt1 = aTwoTrackRes.fTpcEnt1;
	55	fTpcEnt2 = aTwoTrackRes.fTpcEnt2;
	56	fTpcDist = aTwoTrackRes.fTpcDist;
	57	fOutFilename = aTwoTrackRes.fOutFilename;
	58	}
	59
	60	AliTwoTrackRes& AliTwoTrackRes::operator=(const AliTwoTrackRes& aTwoTrackRes)
	61	{
	62	// Assignment operator
	63	if (this == &aTwoTrackRes)
	64	return *this;
65	fChain = aTwoTrackRes.fChain;
66	fESDEvent = aTwoTrackRes.fESDEvent;
67	fOutContainer = aTwoTrackRes.fOutContainer;
68	fTrackCuts = aTwoTrackRes.fTrackCuts;
69	fNTuple1 = aTwoTrackRes.fNTuple1;
70	fNTuple2 = aTwoTrackRes.fNTuple2;
71	fP1 = aTwoTrackRes.fP1;
72	fP2 = aTwoTrackRes.fP2;
73	fPb1 = aTwoTrackRes.fPb1;
74	fPb2 = aTwoTrackRes.fPb2;
75	fP = aTwoTrackRes.fP;
76	fQ = aTwoTrackRes.fQ;
77	fTpcEnt1 = aTwoTrackRes.fTpcEnt1;
78	fTpcEnt2 = aTwoTrackRes.fTpcEnt2;
79	fTpcDist = aTwoTrackRes.fTpcDist;
80	fOutFilename = aTwoTrackRes.fOutFilename;
81	return *this;
82	}
83
84	AliTwoTrackRes::~AliTwoTrackRes() {printf("AliTwoTrackRes destroyed\n");}
85
86	void AliTwoTrackRes::ConnectInputData(Option_t *) {
87	//______________________________________________________________________________
88	// Connect input data and initialize track cuts
89
90	fChain = (TChain*)GetInputData(0);
91	fESDEvent = new AliESDEvent();
92	fESDEvent->ReadFromTree(fChain);
93
94	// Cuts to select primary tracks (ITS+TPC)
95	fTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
96	Double_t cov1, cov2, cov3, cov4, cov5; // diagonal cov. matrix elements
97	Double_t nSigma; // max. DCA to primary vertex
98	Int_t minNClustersTPC; // min. number of clusters per TPC tracks
99	Double_t maxChi2PerClusterTPC; // max. chi2 per cluster per TPC track
100	Int_t cutMode = 1; // select cut mode
101	if (cutMode == 1) {
102	cov1 = 2; cov2 = 2; cov3 = 0.5; cov4 = 0.5; cov5 = 2;
103	nSigma = 3; minNClustersTPC = 75; maxChi2PerClusterTPC = 3.5;
104	fTrackCuts->SetMaxCovDiagonalElements(cov1, cov2, cov3, cov4, cov5);
105	fTrackCuts->SetMaxNsigmaToVertex(nSigma);
106	fTrackCuts->SetRequireSigmaToVertex(kTRUE);
107	fTrackCuts->SetRequireTPCRefit(kTRUE);
108	fTrackCuts->SetAcceptKinkDaughters(kFALSE);
109	fTrackCuts->SetMinNClustersTPC(minNClustersTPC);
110	fTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
111	TString tag("Global tracking");}
112	}
113
114	void AliTwoTrackRes::CreateOutputObjects() {
115	//______________________________________________________________________________
116	// Create output objects
117
118	fNTuple1 = new TNtuple("nt1","True pairs",
119	"pt1:eta1:phi1:nsh1:pt2:eta2:phi2:nsh2:qinv:mindist:dist:corr:qfac");
120	fNTuple2 = new TNtuple("nt2","Mixed pairs",
121	"pt1:eta1:phi1:nsh1:pt2:eta2:phi2:nsh2:qinv:mindist:dist:corr:qfac");
122	Int_t c = 0;
123	fOutContainer = new TObjArray(2);
124	fOutContainer->AddAt(fNTuple1, c++);
125	fOutContainer->AddAt(fNTuple2, c++);
126	}
127
128	void AliTwoTrackRes::Exec(Option_t *) {
129	//______________________________________________________________________________
130	// Create true and mixed pairs keeping some track parameters
131
132	double bfield = 5.0;
133	static int nr=0;
134	if (nr == 0) printf("\tStarting event loop...\n");
135	printf("\rProcessing event %8d", nr);
136	Double_t mpi = 0.13957; // [GeV/c^2]
137	Double_t pidTrk1[AliPID::kSPECIES], pidTrk2[AliPID::kSPECIES];
138	Int_t tpcIn = 80; // [cm]
139	Int_t tpcOut = 250; // [cm]
140	Double_t tdist[170];
141	Double_t tdistrot[170];
142	Double_t tpcEnt1[3], tpcEnt2[3], pos1[3];
143	TVector3 x1, x2, diff;
144	TBits clu1, clu2, sha1, sha2;
145	TRandom2 rnd;
146	Int_t ntracks = fESDEvent->GetNumberOfTracks();
147	for(Int_t itrack = 0; itrack < ntracks; itrack++) {
148	AliESDtrack *track1 = fESDEvent->GetTrack(itrack);
149	AliExternalTrackParam trp1 = const_cast<AliExternalTrackParam>
150	(track1->GetTPCInnerParam());
151	if (!trp1) continue;
152	if (!track1->IsOn(AliESDtrack::kTPCpid)) continue;
153	track1->GetTPCpid(pidTrk1);
154	Int_t q1 = trp1->Charge();
155	if (!((fTrackCuts->AcceptTrack(track1)) && (q1 == 1) &&
156	(pidTrk1[AliPID::kPion]+pidTrk1[AliPID::kMuon] > 0.5))) continue;
157	if (!track1->GetInnerXYZ(tpcEnt1)) continue;
158	clu1 = track1->GetTPCClusterMap();
159	sha1 = track1->GetTPCSharedMap();
160	SetTr1(track1->Pt(), track1->Eta(), track1->Phi(), mpi);
161	SetTpcEnt1(tpcEnt1[0], tpcEnt1[1], tpcEnt1[2]);
162	for(Int_t jtrack = 0; jtrack < itrack; jtrack++) {
163	AliESDtrack *track2 = fESDEvent->GetTrack(jtrack);
164	AliExternalTrackParam trp2 = const_cast<AliExternalTrackParam>
165	(track2->GetTPCInnerParam());
166	if (!trp2) continue;
167	if (!track2->IsOn(AliESDtrack::kTPCpid)) continue;
168	track2->GetTPCpid(pidTrk2);
169	Int_t q2 = trp2->Charge();
170	if (!((fTrackCuts->AcceptTrack(track2)) && (q2 == 1) &&
171	(pidTrk2[AliPID::kPion]+pidTrk2[AliPID::kMuon] > 0.5))) continue;
172	if (!track2->GetInnerXYZ(tpcEnt2)) continue;
173	clu2 = track2->GetTPCClusterMap();
174	sha2 = track2->GetTPCSharedMap();
175	SetTr2(track2->Pt(), track2->Eta(), track2->Phi(), mpi);
176	SetTpcEnt2(tpcEnt2[0], tpcEnt2[1], tpcEnt2[2]);
177	for (Int_t i = tpcIn; i < tpcOut; i++) { // Minimum distance
178	trp1->GetDistance(trp2, (double) i, pos1, bfield);
179	x1.SetXYZ(pos1[0], pos1[1], pos1[2]);
180	tdist[i-tpcIn] = x1.Mag();
181	x1.SetXYZ(-pos1[0], -pos1[1], pos1[2]);
182	tdistrot[i-tpcIn] = x1.Mag();
183	}
184	Double_t mindist = 100000;
185	Int_t jmin=0;
186	for (Int_t j = 0; j < tpcOut-tpcIn; j++) {
187	if (tdist[j] < mindist) {jmin=j; mindist = tdist[j]; }
188	}
189	// Double_t mindist = MinDist(track1, track2);
190	Double_t dist = Dist();
191	Double_t dphi = DPhi();
192	Double_t deta = DEta();
193	Int_t nsh1 = GetNSha(clu1, sha1);
194	Int_t nsh2 = GetNSha(clu2, sha2);
195	Double_t corr = Corr(clu1, clu2, sha1, sha2);
196	Double_t qfac = Qfac(clu1, clu2, sha1, sha2);
197	if ((TMath::Abs(track1->Eta())>0.8)&&(TMath::Abs(track2->Eta())>0.8)) continue;
198	if ((TMath::Abs(dphi)<0.35)&&(deta<0.35)) {
199	FillNTuple1(mindist,dist,corr,qfac,nsh1,nsh2);} // True
200	Double_t tr2rot = RotTr2Phi(); // Rotate trck2
201	SetTr2(track2->Pt(), track2->Eta(), tr2rot, mpi);
202	tpcEnt2[0] = -tpcEnt2[0];
203	tpcEnt2[1] = -tpcEnt2[1];
204	Double_t distrot = Dist();
205	Double_t dphirot = DPhi();
206	Double_t mindistrot = 100000;
207	jmin=0;
208	for (Int_t j = 0; j < tpcOut-tpcIn; j++) {
209	if (tdistrot[j] < mindistrot) {jmin=j; mindistrot = tdistrot[j]; }
210	}
211	if ((TMath::Abs(dphirot)<0.35)&&(deta<0.35)) {
212	if (rnd.Rndm() < 0.5) NoSwap();
213	else Swap();
214	FillNTuple2(mindistrot,distrot,corr,qfac,nsh1,nsh2);} // Mixed
215	}
216	}
217	PostData(0, fOutContainer);
218	nr++;
219	}
220
221	void AliTwoTrackRes::Terminate(Option_t *) {
222	//______________________________________________________________________________
223	// Write output and clean up
224
225	fOutContainer = (TObjArray*)GetOutputData(0);
226	TFile *f1 = new TFile( fOutFilename, "RECREATE" );
227	fOutContainer->Write();
228	f1->Flush();
229	f1->Close();
230	delete f1;
231	delete fChain;
232	delete fNTuple1;
233	delete fNTuple2;
234	printf("\n");
235	}
236
237	//______________________________________________________________________________
238	// Miscellaneous methods
239
240	// Set tracks
241	void AliTwoTrackRes::SetTr1(double pt1, double eta1, double phi1, double m) {
242	fP1.SetPtEtaPhiM(pt1, eta1, phi1, m);}
243	void AliTwoTrackRes::SetTr2(double pt2, double eta2, double phi2, double m) {
244	fP2.SetPtEtaPhiM(pt2, eta2, phi2, m);}
245
246	// Set nominal TPC entrance coordinates
247	void AliTwoTrackRes::SetTpcEnt1(double x1, double y1, double z1) {
248	fTpcEnt1.SetX(x1); fTpcEnt1.SetY(y1); fTpcEnt1.SetZ(z1);}
249	void AliTwoTrackRes::SetTpcEnt2(double x2, double y2, double z2) {
250	fTpcEnt2.SetX(x2); fTpcEnt2.SetY(y2); fTpcEnt2.SetZ(z2);}
251
252	double AliTwoTrackRes::MinDist(AliExternalTrackParam *trk1,
253	AliExternalTrackParam *trk2) {
254	// Calculate minimum track separation within the TPC
255
256	int tpcIn = 0; // [cm]
257	int tpcOut = 170; // [cm]
258	double tdist[170], pos[3];
259	TVector3 x;
260	for (int i = tpcIn; i < tpcOut; i++) {
261	trk1->GetDistance(trk2, i, pos, 5000);
262	x.SetXYZ(pos[0], pos[1], pos[2]);
263	tdist[i-tpcIn] = x.Mag();
264	}
265	double maxdist = 0.0;
266	for (int j = 0; j < tpcOut-tpcIn; j++) {
267	if (tdist[j] > maxdist) { maxdist = tdist[j]; }
268	}
269	double mindist = maxdist;
270	for (int j = 0; j < tpcOut-tpcIn; j++) {
271	if (tdist[j] < mindist) { mindist = tdist[j]; }
272	}
273	return mindist;}
274
275	int AliTwoTrackRes::GetNSha(TBits cl, TBits sh) {
276	// Get number of shared clusters
277
278	int ncl = cl.GetNbits();
279	int sum = 0;
280	for(int i = 0; i < ncl; i++) {
281	if (!cl.TestBitNumber(i)) continue;
282	int n = sh.TestBitNumber(i);
283	sum += n;}
284	return sum;}
285
286	double AliTwoTrackRes::Corr(TBits cl1, TBits cl2, TBits sh1, TBits sh2) {
287	// Calculate correlation coefficient
288
289	int ncl1 = cl1.GetNbits();
290	int ncl2 = cl2.GetNbits();
291	double sumN = 0; double sumX = 0; double sumY = 0;
292	double sumXX = 0; double sumYY = 0; double sumXY = 0; double corr = -2.0;
293	for(int i = 0; i < ncl1 && i < ncl2; i++) {
294	if (!(cl1.TestBitNumber(i)&&cl2.TestBitNumber(i))) continue;
295	int x = sh1.TestBitNumber(i);
296	int y = sh2.TestBitNumber(i);
297	sumN += 1.0;
298	sumX += x;
299	sumY += y;
300	sumXX += x*x;
301	sumYY += y*y;
302	sumXY += x*y;
303	}
304	double meanX = sumX/sumN;
305	double meanY = sumY/sumN;
306	double meanXX = sumXX/sumN;
307	double meanYY = sumYY/sumN;
308	double meanXY = sumXY/sumN;
309	double sX = TMath::Sqrt(TMath::Abs(meanXX-meanX*meanX));
310	double sY = TMath::Sqrt(TMath::Abs(meanYY-meanY*meanY));
311	if (sXsY!=0) corr = (meanXY-meanXmeanY)/(sX*sY);
312	return corr;}
313
314	double AliTwoTrackRes::Qfac(TBits cl1, TBits cl2, TBits sh1, TBits sh2) {
315	// Quality factor from AliFemto
316
317	int ncl1 = cl1.GetNbits();
318	int ncl2 = cl2.GetNbits();
319	int sumCls = 0; int sumSha = 0; int sumQ = 0;
320	double shfrac = 0; double qfactor = 0;
321	for(int i = 0; i < ncl1 && i < ncl2; i++) {
322	if (cl1.TestBitNumber(i) && cl2.TestBitNumber(i)) { // Both clusters
323	if (sh1.TestBitNumber(i) && sh2.TestBitNumber(i)) { // Shared
324	sumQ++;
325	sumCls+=2;
326	sumSha+=2;}
327	else {sumQ--; sumCls+=2;}
328	}
329	else if (cl1.TestBitNumber(i) \|\| cl2.TestBitNumber(i)) { // Non shared
330	sumQ++;
331	sumCls++;}
332	}
333	if (sumCls>0) {
334	qfactor = sumQ*1.0/sumCls;
335	shfrac = sumSha*1.0/sumCls;
336	}
337	return qfactor;
338	}
339
340	// Rotate second track for mixed pairs
341	double AliTwoTrackRes::RotTr2Phi() {
342	double rot = TVector2::Phi_mpi_pi(fP2.Phi()+TMath::Pi());
343	fTpcEnt2.SetPhi(TVector2::Phi_mpi_pi(fTpcEnt2.Phi()+TMath::Pi()));
344	return rot;}
345
346	// Fill NTuples
347	void AliTwoTrackRes::FillNTuple1(double minsep, double sep, double corr,
348	double qf, int ns1, int ns2) {
349	fNTuple1->Fill(fP1.Pt(),fP1.Eta(),fP1.Phi(),ns1,fP2.Pt(),fP2.Eta(),
350	fP2.Phi(),ns2,Qinv(),minsep,sep,corr,qf);}
351	void AliTwoTrackRes::FillNTuple2(double minsep, double sep, double corr,
352	double qf, int ns1, int ns2) {
353	fNTuple2->Fill(fPb1.Pt(),fPb1.Eta(),fPb1.Phi(),ns1,fPb2.Pt(),fPb2.Eta(),
354	fPb2.Phi(),ns2,Qinv(),minsep,sep,corr,qf);}
355
356	//______________________________________________________________________________
357	// EOF
358
359