1 ////////////////////////////////////////////////////////////////////////////////
2 //345678901234567890123456789012345678901234567890123456789012345678901234567890
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
8 ////////////////////////////////////////////////////////////////////////////////
10 #include "AliTwoTrackRes.h"
21 ClassImp(AliTwoTrackRes)
23 //______________________________________________________________________________
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()
32 DefineInput(0, TChain::Class()); // Slot input 0 reads from a TChain
33 DefineOutput(0, TObjArray::Class()); // Slot output 0 writes into a TObjArray
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()
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;
54 fTpcEnt1 = aTwoTrackRes.fTpcEnt1;
55 fTpcEnt2 = aTwoTrackRes.fTpcEnt2;
56 fTpcDist = aTwoTrackRes.fTpcDist;
57 fOutFilename = aTwoTrackRes.fOutFilename;
60 AliTwoTrackRes& AliTwoTrackRes::operator=(const AliTwoTrackRes& aTwoTrackRes)
62 // Assignment operator
63 if (this == &aTwoTrackRes)
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;
77 fTpcEnt1 = aTwoTrackRes.fTpcEnt1;
78 fTpcEnt2 = aTwoTrackRes.fTpcEnt2;
79 fTpcDist = aTwoTrackRes.fTpcDist;
80 fOutFilename = aTwoTrackRes.fOutFilename;
84 AliTwoTrackRes::~AliTwoTrackRes() {printf("AliTwoTrackRes destroyed\n");}
86 void AliTwoTrackRes::ConnectInputData(Option_t *) {
87 //______________________________________________________________________________
88 // Connect input data and initialize track cuts
90 fChain = (TChain*)GetInputData(0);
91 fESDEvent = new AliESDEvent();
92 fESDEvent->ReadFromTree(fChain);
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
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");}
114 void AliTwoTrackRes::CreateOutputObjects() {
115 //______________________________________________________________________________
116 // Create output objects
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");
123 fOutContainer = new TObjArray(2);
124 fOutContainer->AddAt(fNTuple1, c++);
125 fOutContainer->AddAt(fNTuple2, c++);
128 void AliTwoTrackRes::Exec(Option_t *) {
129 //______________________________________________________________________________
130 // Create true and mixed pairs keeping some track parameters
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]
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;
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());
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());
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();
184 Double_t mindist = 100000;
186 for (Int_t j = 0; j < tpcOut-tpcIn; j++) {
187 if (tdist[j] < mindist) {jmin=j; mindist = tdist[j]; }
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;
208 for (Int_t j = 0; j < tpcOut-tpcIn; j++) {
209 if (tdistrot[j] < mindistrot) {jmin=j; mindistrot = tdistrot[j]; }
211 if ((TMath::Abs(dphirot)<0.35)&&(deta<0.35)) {
212 if (rnd.Rndm() < 0.5) NoSwap();
214 FillNTuple2(mindistrot,distrot,corr,qfac,nsh1,nsh2);} // Mixed
217 PostData(0, fOutContainer);
221 void AliTwoTrackRes::Terminate(Option_t *) {
222 //______________________________________________________________________________
223 // Write output and clean up
225 fOutContainer = (TObjArray*)GetOutputData(0);
226 TFile *f1 = new TFile( fOutFilename, "RECREATE" );
227 fOutContainer->Write();
237 //______________________________________________________________________________
238 // Miscellaneous methods
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);}
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);}
252 double AliTwoTrackRes::MinDist(AliExternalTrackParam *trk1,
253 AliExternalTrackParam *trk2) {
254 // Calculate minimum track separation within the TPC
256 int tpcIn = 0; // [cm]
257 int tpcOut = 170; // [cm]
258 double tdist[170], pos[3];
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();
265 double maxdist = 0.0;
266 for (int j = 0; j < tpcOut-tpcIn; j++) {
267 if (tdist[j] > maxdist) { maxdist = tdist[j]; }
269 double mindist = maxdist;
270 for (int j = 0; j < tpcOut-tpcIn; j++) {
271 if (tdist[j] < mindist) { mindist = tdist[j]; }
275 int AliTwoTrackRes::GetNSha(TBits cl, TBits sh) {
276 // Get number of shared clusters
278 int ncl = cl.GetNbits();
280 for(int i = 0; i < ncl; i++) {
281 if (!cl.TestBitNumber(i)) continue;
282 int n = sh.TestBitNumber(i);
286 double AliTwoTrackRes::Corr(TBits cl1, TBits cl2, TBits sh1, TBits sh2) {
287 // Calculate correlation coefficient
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);
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 (sX*sY!=0) corr = (meanXY-meanX*meanY)/(sX*sY);
314 double AliTwoTrackRes::Qfac(TBits cl1, TBits cl2, TBits sh1, TBits sh2) {
315 // Quality factor from AliFemto
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
327 else {sumQ--; sumCls+=2;}
329 else if (cl1.TestBitNumber(i) || cl2.TestBitNumber(i)) { // Non shared
334 qfactor = sumQ*1.0/sumCls;
335 shfrac = sumSha*1.0/sumCls;
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()));
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);}
356 //______________________________________________________________________________