1 /////////////////////////////////////////////////////////////////////////////
3 // AliFemtoShareQualityPairCut - a pair cut which checks for some pair //
4 // qualities that attempt to identify slit/doubly reconstructed tracks //
6 /////////////////////////////////////////////////////////////////////////////
7 /***************************************************************************
9 * $Id: AliFemtoShareQualityPairCut.cxx 50722 2011-07-21 15:18:38Z akisiel $
11 * Author: Adam Kisiel, Ohio State, kisiel@mps.ohio-state.edu
12 ***************************************************************************
14 * Description: part of STAR HBT Framework: AliFemtoMaker package
15 * a cut to remove "shared" and "split" pairs
17 ***************************************************************************
20 **************************************************************************/
22 #include "AliFemtoV0TrackPairCut.h"
27 ClassImp(AliFemtoV0TrackPairCut)
31 AliFemtoV0TrackPairCut::AliFemtoV0TrackPairCut():
35 fShareQualityMax(1.0),
36 fShareFractionMax(1.0),
43 fFirstParticleType(kLambda),
44 fSecondParticleType(kProton),
45 fMinAvgSepTrackPos(0),
48 // Default constructor
52 AliFemtoV0TrackPairCut::~AliFemtoV0TrackPairCut(){
55 AliFemtoV0TrackPairCut& AliFemtoV0TrackPairCut::operator=(const AliFemtoV0TrackPairCut& cut)
59 AliFemtoPairCut::operator=(cut);
63 fShareQualityMax = 1.0;
64 fShareFractionMax = 1.0;
70 fMinAvgSepTrackPos = 0;
71 fMinAvgSepTrackNeg = 0;
77 bool AliFemtoV0TrackPairCut::Pass(const AliFemtoPair* pair){
78 // Check for pairs that are possibly shared/double reconstruction
84 if(!(pair->Track1()->V0() && pair->Track2()->Track()))
90 if( (-(pair->Track1()->V0()->IdNeg()+1)) ==pair->Track2()->TrackId() || (-(pair->Track1()->V0()->IdPos()+1)) ==pair->Track2()->TrackId())
97 if(pair->Track1()->V0()->IdNeg()==pair->Track2()->TrackId() || pair->Track1()->V0()->IdPos()==pair->Track2()->TrackId())
103 bool tempTPCEntrancePos = true;
104 bool tempTPCEntranceNeg = true;
105 bool tempTPCExitPos = true;
106 bool tempTPCExitNeg = true;
107 if(fDataType==kESD || fDataType==kAOD)
109 double distx = pair->Track1()->V0()->NominalTpcEntrancePointPos().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
110 double disty = pair->Track1()->V0()->NominalTpcEntrancePointPos().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
111 double distz = pair->Track1()->V0()->NominalTpcEntrancePointPos().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
112 double distPos = sqrt(distx*distx + disty*disty + distz*distz);
114 distx = pair->Track1()->V0()->NominalTpcEntrancePointNeg().x() - pair->Track2()->Track()->NominalTpcEntrancePoint().x();
115 disty = pair->Track1()->V0()->NominalTpcEntrancePointNeg().y() - pair->Track2()->Track()->NominalTpcEntrancePoint().y();
116 distz = pair->Track1()->V0()->NominalTpcEntrancePointNeg().z() - pair->Track2()->Track()->NominalTpcEntrancePoint().z();
117 double distNeg = sqrt(distx*distx + disty*disty + distz*distz);
119 double distExitX = pair->Track1()->V0()->NominalTpcExitPointPos().x() - pair->Track2()->Track()->NominalTpcExitPoint().x();
120 double distExitY = pair->Track1()->V0()->NominalTpcExitPointPos().y() - pair->Track2()->Track()->NominalTpcExitPoint().y();
121 double distExitZ = pair->Track1()->V0()->NominalTpcExitPointPos().z() - pair->Track2()->Track()->NominalTpcExitPoint().z();
122 double distExitPos = sqrt(distExitX*distExitX + distExitY*distExitY + distExitZ*distExitZ);
124 distExitX = pair->Track1()->V0()->NominalTpcExitPointNeg().x() - pair->Track2()->Track()->NominalTpcExitPoint().x();
125 distExitY = pair->Track1()->V0()->NominalTpcExitPointNeg().y() - pair->Track2()->Track()->NominalTpcExitPoint().y();
126 distExitZ = pair->Track1()->V0()->NominalTpcExitPointNeg().z() - pair->Track2()->Track()->NominalTpcExitPoint().z();
127 double distExitNeg = sqrt(distExitX*distExitX + distExitY*distExitY + distExitZ*distExitZ);
129 tempTPCEntrancePos = distPos > fDTPCMin;
130 tempTPCEntranceNeg = distNeg > fDTPCMin;
131 tempTPCExitPos = distExitPos > fDTPCExitMin;
132 tempTPCExitNeg = distExitNeg > fDTPCExitMin;
135 if(!tempTPCEntrancePos || !tempTPCEntranceNeg || !tempTPCExitPos || !tempTPCExitNeg) return false;
137 //reject merged trakcs in TPC
140 //temp = dist > fDTPCMin;
145 //kopia z AliFemtoShareQualityPairCut.cxx
152 if ((fShareFractionMax < 1.0) && ( fShareQualityMax < 1.0)) {
153 for (unsigned int imap=0; imap<pair->Track1()->V0()->TPCclustersPos().GetNbits(); imap++) {
154 // If both have clusters in the same row
155 if (pair->Track1()->V0()->TPCclustersPos().TestBitNumber(imap) &&
156 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
157 // Do they share it ?
158 if (pair->Track1()->V0()->TPCsharingPos().TestBitNumber(imap) &&
159 pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
161 // cout << "A shared cluster !!!" << endl;
162 // cout << "imap idx1 idx2 "
164 // << tP1idx[imap] << " " << tP2idx[imap] << endl;
170 // Different hits on the same padrow
176 else if (pair->Track1()->V0()->TPCclustersPos().TestBitNumber(imap) ||
177 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
178 // One track has a hit, the other does not
184 Float_t hsmval = 0.0;
185 Float_t hsfval = 0.0;
191 // if (hsmval > -0.4) {
192 // cout << "Pair quality: " << hsmval << " " << an << " " << nh << " "
193 // << (pair->Track1()->Track()) << " "
194 // << (pair->Track2()->Track()) << endl;
195 // cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl;
197 // if (hsfval > 0.0) {
198 // cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << " " << hsmval << " " << an << " " << nh << endl;
201 temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax);
202 if(!temp) return false;
208 for (unsigned int imap=0; imap<pair->Track1()->V0()->TPCclustersNeg().GetNbits(); imap++) {
209 // If both have clusters in the same row
210 if (pair->Track1()->V0()->TPCclustersNeg().TestBitNumber(imap) &&
211 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
212 // Do they share it ?
213 if (pair->Track1()->V0()->TPCsharingNeg().TestBitNumber(imap) &&
214 pair->Track2()->Track()->TPCsharing().TestBitNumber(imap))
216 // cout << "A shared cluster !!!" << endl;
217 // cout << "imap idx1 idx2 "
219 // << tP1idx[imap] << " " << tP2idx[imap] << endl;
225 // Different hits on the same padrow
231 else if (pair->Track1()->V0()->TPCclustersNeg().TestBitNumber(imap) ||
232 pair->Track2()->Track()->TPCclusters().TestBitNumber(imap)) {
233 // One track has a hit, the other does not
246 // if (hsmval > -0.4) {
247 // cout << "Pair quality: " << hsmval << " " << an << " " << nh << " "
248 // << (pair->Track1()->Track()) << " "
249 // << (pair->Track2()->Track()) << endl;
250 // cout << "Bits: " << pair->Track1()->Track()->TPCclusters().GetNbits() << endl;
252 // if (hsfval > 0.0) {
253 // cout << "Pair sharity: " << hsfval << " " << ns << " " << nh << " " << hsmval << " " << an << " " << nh << endl;
256 temp = (hsmval < fShareQualityMax) && (hsfval < fShareFractionMax);
267 AliFemtoThreeVector first, second, tmp;
268 for(int i=0; i<8 ;i++)
270 tmp = pair->Track1()->V0()->NominalTpcPointPos(i);
271 //cout<<"X pos: "<<tmp.x()<<endl;
272 first.SetX((double)(tmp.x()));
273 first.SetY((double)tmp.y());
274 first.SetZ((double)tmp.z());
276 tmp = pair->Track2()->Track()->NominalTpcPoint(i);
277 second.SetX((double)tmp.x());
278 second.SetY((double)tmp.y());
279 second.SetZ((double)tmp.z());
281 avgSep += TMath::Sqrt(((double)first.x()-(double)second.x())*((double)first.x()-(double)second.x())+((double)first.y()-(double)second.y())*((double)first.y()-second.y())+((double)first.z()-(double)second.z())*((double)first.z()-(double)second.z()));
285 if(avgSep<fMinAvgSepTrackPos) return false;
289 for(int i=0; i<8 ;i++)
291 tmp = pair->Track1()->V0()->NominalTpcPointNeg(i);
292 //cout<<"X pos: "<<tmp.x()<<endl;
293 first.SetX((double)(tmp.x()));
294 first.SetY((double)tmp.y());
295 first.SetZ((double)tmp.z());
297 tmp = pair->Track2()->Track()->NominalTpcPoint(i);
298 second.SetX((double)tmp.x());
299 second.SetY((double)tmp.y());
300 second.SetZ((double)tmp.z());
302 avgSep += TMath::Sqrt(((double)first.x()-(double)second.x())*((double)first.x()-(double)second.x())+((double)first.y()-(double)second.y())*((double)first.y()-second.y())+((double)first.z()-(double)second.z())*((double)first.z()-(double)second.z()));
306 if(avgSep<fMinAvgSepTrackNeg) return false;
311 //Qinv cut (we are trying to get rid of antiresidual correlation between primary protons)
315 //2 daughters of first V0
316 double mom1PosX = pair->Track1()->V0()->MomPosX();
317 double mom1PosY = pair->Track1()->V0()->MomPosY();
318 double mom1PosZ = pair->Track1()->V0()->MomPosZ();
319 double mom1NegX = pair->Track1()->V0()->MomNegX();
320 double mom1NegY = pair->Track1()->V0()->MomNegY();
321 double mom1NegZ = pair->Track1()->V0()->MomNegZ();
324 //double PionMass = 0.13956995;
325 //double KaonMass = 0.493677;
326 double ProtonMass = 0.938272;
327 //double LambdaMass = 1.115683;
329 AliFemtoLorentzVector fFourMomentum1; // Particle momentum
330 AliFemtoLorentzVector fFourMomentum2; // Particle momentum
332 AliFemtoThreeVector temp1;
334 if(fFirstParticleType == 0) //lambda
336 if(fSecondParticleType == 2) //proton
338 //temp1 = ::sqrt(mom1PosX*mom1PosX+mom1PosY*mom1PosY+mom1PosZ*mom1PosZ);
339 temp1.SetX(mom1PosX);
340 temp1.SetY(mom1PosY);
341 temp1.SetZ(mom1PosZ);
342 ener1 = ::sqrt(temp1.Mag2()+ProtonMass*ProtonMass);
345 else if(fFirstParticleType == 1) //antilambda
347 if(fSecondParticleType == 3) //antiproton
349 //temp1 = ::sqrt(mom1NegX*mom1NegX+mom1NegY*mom1NegY+mom1NegZ*mom1NegZ);
350 temp1.SetX(mom1NegX);
351 temp1.SetY(mom1NegY);
352 temp1.SetZ(mom1NegZ);
353 ener1 = ::sqrt(temp1.Mag2()+ProtonMass*ProtonMass);
356 fFourMomentum1.SetVect(temp1);
357 fFourMomentum1.SetE(ener1);
359 //AliFemtoLorentzVector fFourMomentum2; // Particle momentum
360 AliFemtoThreeVector temp2;
363 //2 daughters of second V0
364 temp2 = pair->Track2()->Track()->P();
367 if(fSecondParticleType == 2 || fSecondParticleType == 3) //proton
369 ener2 = ::sqrt(temp2.Mag2()+ProtonMass*ProtonMass);
372 fFourMomentum2.SetVect(temp2);
373 fFourMomentum2.SetE(ener2);
376 AliFemtoLorentzVector tDiff = (fFourMomentum1-fFourMomentum2);
378 double tQinv = fabs(-1.*tDiff.m()); // note - qInv() will be negative for identical pairs...
380 //cout<<"tQinv/2: "<<tQinv/2<<endl;
381 if(tQinv/2 < fKstarCut)
390 AliFemtoString AliFemtoV0TrackPairCut::Report(){
391 // Prepare the report from the execution
392 string stemp = "AliFemtoV0 Pair Cut - remove shared and split pairs\n"; char ctemp[100];
393 snprintf(ctemp , 100, "Number of pairs which passed:\t%ld Number which failed:\t%ld\n",fNPairsPassed,fNPairsFailed);
395 AliFemtoString returnThis = stemp;
399 void AliFemtoV0TrackPairCut::SetV0Max(Double_t aV0Max) {
403 Double_t AliFemtoV0TrackPairCut::GetAliFemtoV0Max() const {
408 TList *AliFemtoV0TrackPairCut::ListSettings()
410 // return a list of settings in a writable form
411 TList *tListSetttings = new TList();
413 snprintf(buf, 200, "AliFemtoV0TrackPairCut.sharequalitymax=%f", fV0Max);
414 snprintf(buf, 200, "AliFemtoV0TrackPairCut.sharefractionmax=%f", fShareFractionMax);
415 tListSetttings->AddLast(new TObjString(buf));
417 return tListSetttings;
420 void AliFemtoV0TrackPairCut::SetRemoveSameLabel(Bool_t aRemove)
422 fRemoveSameLabel = aRemove;
425 void AliFemtoV0TrackPairCut::SetTPCOnly(Bool_t tpconly)
427 fTrackTPCOnly = tpconly;
430 void AliFemtoV0TrackPairCut::SetShareQualityMax(Double_t aShareQualityMax) {
431 fShareQualityMax = aShareQualityMax;
434 void AliFemtoV0TrackPairCut::SetShareFractionMax(Double_t aShareFractionMax) {
435 fShareFractionMax = aShareFractionMax;
438 void AliFemtoV0TrackPairCut::SetDataType(AliFemtoDataType type)
443 void AliFemtoV0TrackPairCut::SetTPCEntranceSepMinimum(double dtpc)
448 void AliFemtoV0TrackPairCut::SetTPCExitSepMinimum(double dtpc)
453 void AliFemtoV0TrackPairCut::SetKstarCut(double kstar, AliFemtoParticleType firstParticle, AliFemtoParticleType secondParticle)
456 fFirstParticleType = firstParticle; //for kstar - first particle type (V0 type)
457 fSecondParticleType = secondParticle;
460 void AliFemtoV0TrackPairCut::SetMinAvgSeparation(int type, double minSep)
462 if(type == 0) //Track-Pos
463 fMinAvgSepTrackPos = minSep;
464 else if(type == 1) //Track-Neg
465 fMinAvgSepTrackNeg = minSep;