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

////////////////////////////////////////////////////////////////////////////////
//345678901234567890123456789012345678901234567890123456789012345678901234567890
//       1         2         3         4         5         6         7         8
//
// Tool to study two-track effects in ALICE for femtoscopic studies
// J. Mercado. Last modified: 28.12.2009
//
////////////////////////////////////////////////////////////////////////////////

#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()
{
  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()
{
  //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;
}

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