StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructPairCuts.h
1 /**********************************************************************
2  *
3  * $Id: StEStructPairCuts.h,v 1.21 2012/11/16 21:22:27 prindle Exp $
4  *
5  * Author: Jeff Porter
6  *
7  **********************************************************************
8  *
9  * Description: Cut class for track-pair level quantities
10  *
11  *
12  ***********************************************************************/
13 #ifndef __STEBYEPAIRCUTS_H
14 #define __STEBYEPAIRCUTS_H
15 
16 #include "StEStructPairLUT.h"
17 #include "StEStructPool/AnalysisMaker/StEStructCuts.h"
18 #include "StEStructPool/EventMaker/StEStructTrack.h"
19 
20 #include "Stiostream.h"
21 
22 
23 class StEStructPairCuts : public StEStructCuts {
24 
25 protected:
26 
27 
28  CutName mdphiName;
29  CutName mdetaName;
30  CutName mgooddzdxyName;
31  CutName mdmtName;
32  CutName mqInvName;
33  CutName mEntSepName;
34  CutName mExitSepName;
35  CutName mQualityName;
36  CutName mMidTpcSepLSName;
37  CutName mMidTpcSepUSName;
38  CutName mHBTName;
39  CutName mCoulombName;
40  CutName mMergingName;
41  CutName mMergingName2;
42  CutName mCrossingName;
43  CutName mCrossingName2;
44  CutName mLUTName;
45  CutName mpionMomentumName;
46  CutName mKaonMomentumName;
47  CutName mprotonMomentumName;
48 
49  CutName mpionOtherMassName;
50  CutName mpionpionMassName;
51  CutName mpionKaonMassName;
52  CutName mpionprotonMassName;
53  CutName mKaonOtherMassName;
54  CutName mKaonKaonMassName;
55  CutName mKaonprotonMassName;
56  CutName mprotonOtherMassName;
57  CutName mprotonprotonMassName;
58  CutName mOtherOtherMassName;
59 
60  float mdphi[2];
61  float mdeta[2];
62  float mgooddzdxy[2];
63  float mdmt[2];
64  float mqInv[2];
65  float mEntSep[2];
66  float mExitSep[2];
67  float mQuality[2];
68  float mMidTpcSepLS[2];
69  float mMidTpcSepUS[2];
70  float mHBT[4];
71  float mCoulomb[4];
72  float mMerging[2];
73  float mMerging2[2];
74  float mCrossing[2];
75  float mCrossing2[2];
76  float mLUTParams[2];
77  float mdEdxMomentumCut[4][2];
78  float mToFMomentumCut[4][2];
79  float mpionOtherMass[2];
80  float mpionpionMass[2];
81  float mpionKaonMass[2];
82  float mpionprotonMass[2];
83  float mKaonOtherMass[2];
84  float mKaonKaonMass[2];
85  float mKaonprotonMass[2];
86  float mprotonOtherMass[2];
87  float mprotonprotonMass[2];
88  float mOtherOtherMass[2];
89 
90  bool mdeltaPhiCut, mdeltaEtaCut, mGooddeltaZdeltaXYCut, mdeltaMtCut;
91  bool mqInvCut, mEntSepCut, mExitSepCut, mQualityCut;
92  bool mMidTpcSepLSCut, mMidTpcSepUSCut;
93  bool mHBTCut, mCoulombCut, mMergingCut, mCrossingCut, mMergingCut2, mCrossingCut2, mLUTCut;
94  bool mpionMomentumCut, mKaonMomentumCut, mprotonMomentumCut;
95  bool mpionOtherMassCut, mpionpionMassCut, mpionKaonMassCut, mpionprotonMassCut;
96  bool mKaonOtherMassCut, mKaonKaonMassCut, mKaonprotonMassCut, mprotonOtherMassCut;
97  bool mprotonprotonMassCut, mOtherOtherMassCut;
98 
99  int mdphiCounter[4], mdetaCounter[4], mgooddzdxyCounter[4], mdmtCounter[4];
100  int mqInvCounter[4], mEntSepCounter[4], mExitSepCounter[4], mQualityCounter[4];
101  int msplitLSCounter[4], msplitUSCounter[4];
102  int mHBTCounter[4], mCoulombCounter[4], mMergingCounter[4], mCrossingCounter[4], mMergingCounter2[4], mCrossingCounter2[4], mLUTCounter[4];
103  int mpionMomentumCounter[4], mKaonMomentumCounter[4], mprotonMomentumCounter[4];
104  int mpionOtherMassCounter[4], mpionpionMassCounter[4], mpionKaonMassCounter[4], mpionprotonMassCounter[4];
105  int mKaonOtherMassCounter[4], mKaonKaonMassCounter[4], mKaonprotonMassCounter[4], mprotonOtherMassCounter[4];
106  int mprotonprotonMassCounter[4], mOtherOtherMassCounter[4];
107 
108  // if data stored for subsequent use.... e.g. next cut, histogramming
109  float mdeltaPhi, mdeltaEta, mdeltaMt;
110  float mqInvariant, mEntranceSeparation, mExitSeparation, mQualityVal;
111  float mMidTpcSeparationLS, mMidTpcSeparationUS;
112  float mBField; // magnetic field (kilogauss as MuDst)
113 
114  int mType;
115  unsigned long mapMask0;
116  unsigned long mapMask1;
117  unsigned long bitI[32];
118 
119  int mretVal;
120 
121  int mcutMode; // cut space definitions
122 
123  float mZoffset; // Z offset in primary vertex; correction for mixed event pair separations
124 
125  void init();
126  void initCuts();
127  void initNames();
128 
129 public:
130 
131  StEStructTrack* mTrack1;
132  StEStructTrack* mTrack2;
133 
134  StEStructPairCuts();
135  StEStructPairCuts(const char* cutFileName);
136  virtual ~StEStructPairCuts();
137 
138  // Look Up Table object.
139  StEStructPairLUT *mLUT;
140 
141  virtual bool loadBaseCuts(const char* name, const char** vals, int nvals);
142  virtual void loadUserCuts(const char* name, const char** vals, int nvals);
143  virtual void printCutStats(ostream& ofs);
144  virtual void printCutCounts(ostream& ofs, const char* cutType,int c1, int c2);
145 
146  float BField() const { return mBField; };
147  void SetBField(const float bfield){ mBField=bfield; };
148 
149  // StEStructPairCuts stuff
150 
151  const StEStructTrack* Track1() const;
152  const StEStructTrack* Track2() const;
153  void SetTrack1(const StEStructTrack* trkPtr);
154  void SetTrack2(const StEStructTrack* trkPtr);
155 
156 
157 //-----For Auto Correlationn Study
158  float DeltaMt() const;
159  float DeltaXt() const;
160  float DeltaYt() const;
161  float DeltaYt(float mass1,float mass2) const;
162  float DeltaEta() const;
163  float DeltaEta(float mass1,float mass2) const;
164  float DeltaPhi() const;
165  float DeltaPt() const;
166  float SigmaMt() const;
167  float SigmaXt() const;
168  float SigmaYt() const;
169  float SigmaYt(float mass1,float mass2) const;
170  float SigmaEta() const;
171  float SigmaEta(float mass1,float mass2) const;
172  float SigmaPhi() const;
173  float SigmaPt() const;
174 //--- For Cut Bin Definitions (found these in more than 1 place ...)
175  bool awaySide();
176  bool sameSide();
177 
178 
179 //-----For HBT Study
180  StLorentzVectorF fourMomentumDiff() const;
181  StLorentzVectorF fourMomentumSum() const;
182  float qInv() const;
183  float kT() const;
184  float mInv() const;
185 
186 //-----For Pair Cut
187  double quality() const;
188  double OpeningAngle() const;
189  double NominalTpcEntranceSeparation() const;
190  double NominalTpcXYEntranceSeparation() const;
191  double NominalTpcZEntranceSeparation() const;
192  double MidTpcXYSeparation() const;
193  double MidTpcZSeparation() const;
194  double MidTpcSeparation() const;
195  double OuterMidTpcXYSeparation() const;
196  double OuterMidTpcZSeparation() const;
197  double OuterMidTpcSeparation() const;
198  double NominalTpcExitSeparation() const;
199  double NominalTpcXYExitSeparation() const;
200  double NominalTpcZExitSeparation() const;
201  double TpcZEndcapExitSeparation() const;
202  double TpcRadialEndCapSeparation() const;
203  double NominalTpcAvgXYSeparation() const;
204  double NominalTpcAvgZSeparation() const;
205  bool TracksCrossInPhi() const;
206 
207  void SetZoffset(float Zoffset) { mZoffset = Zoffset; };
208  float GetZoffset() { return mZoffset; };
209 
210  void setPairType(int type);
211  int cutPair();
212  int cutPair(bool doHistograms);
213  int cutPairHistograms();
214 
215  // explicit cut coding...
216  // for Pair cuts we have to break the model of
217  // non-duplicate codes simply because even checking a
218  // bool for each cut is expensive when done N1*N2*Ncut times
219 
220  /* because of the expense 1st check wether we need to do further */
221 
222  int goodDeltaXY();
223  int goodDeltaZ();
224  int goodDeltaMt();
225 
226  int cutDeltaPhi();
227  int cutDeltaEta();
228  int cutDeltaMt();
229  int cutqInvORNominalEntranceSep();
230  int cutqInv();
231  int cutEntranceSep();
232  int cutMidTpcSepUS();
233  int cutMidTpcSepLS();
234  int cutExitSep();
235  int cutQuality();
236  int cutHBT();
237  int cutCoulomb();
238  int cutMerging();
239  int cutCrossing();
240  int cutMerging2();
241  int cutCrossing2();
242  int cutLUT();
243  int cutMass();
244 
245  // calls above set but fills histogramming variables
246  int cutDeltaPhiH();
247  int cutDeltaEtaH();
248  int cutDeltaMtH();
249  int cutqInvH();
250  int cutEntranceSepH();
251  int cutMidTpcSepUSH();
252  int cutMidTpcSepLSH();
253  int cutExitSepH();
254  int cutQualityH();
255 
256  // for new GoodPair study
257 
258  int correlationDepth();
259 
260 
261 
262  ClassDef(StEStructPairCuts,1)
263 };
264 
265 inline void StEStructPairCuts::loadUserCuts(const char* name, const char** vals, int nvals){ }
266 
267 inline void StEStructPairCuts::setPairType(int type) { mType=type; }
268 inline int StEStructPairCuts::cutPair(bool doHistograms){
269  if(!doHistograms) return cutPair();
270  return cutPairHistograms();
271 }
272 
273 
274 inline void StEStructPairCuts::SetTrack1(const StEStructTrack* trkPtr){
275  mTrack1=(StEStructTrack*)trkPtr;
276 }
277 
278 inline void StEStructPairCuts::SetTrack2(const StEStructTrack* trkPtr){
279  mTrack2=(StEStructTrack*)trkPtr;
280 }
281 
282 inline const StEStructTrack* StEStructPairCuts::Track1() const {return mTrack1;}
283 inline const StEStructTrack* StEStructPairCuts::Track2() const {return mTrack2;}
284 
285 inline float StEStructPairCuts::SigmaMt() const {
286  return mTrack1->FourMomentum().mt()+mTrack2->FourMomentum().mt();
287 }
288 inline float StEStructPairCuts::SigmaXt() const {
289  return mTrack1->Xt()+mTrack2->Xt();
290 }
291 inline float StEStructPairCuts::SigmaYt() const {
292  return mTrack1->Yt()+mTrack2->Yt();
293 }
294 inline float StEStructPairCuts::SigmaYt(float mass1, float mass2) const {
295  return mTrack1->Yt(mass1)+mTrack2->Yt(mass2);
296 }
297 inline float StEStructPairCuts::SigmaEta() const {
298  return mTrack1->Eta()+mTrack2->Eta();
299 }
300 inline float StEStructPairCuts::SigmaEta(float mass1, float mass2) const {
301  return mTrack1->Eta(mass1)+mTrack2->Eta(mass2);
302 }
303 inline float StEStructPairCuts::SigmaPhi() const {
304  return mTrack1->Phi()+mTrack2->Phi();
305 }
306 
307 inline float StEStructPairCuts::SigmaPt() const {
308  return mTrack1->Pt()+mTrack2->Pt();
309 }
310 
311 inline float StEStructPairCuts::DeltaMt() const {
312  return mTrack1->FourMomentum().mt()-mTrack2->FourMomentum().mt();
313 }
314 inline float StEStructPairCuts::DeltaXt() const {
315  return mTrack1->Xt()-mTrack2->Xt();
316 }
317 inline float StEStructPairCuts::DeltaYt() const {
318  return mTrack1->Yt()-mTrack2->Yt();
319 }
320 inline float StEStructPairCuts::DeltaYt(float mass1, float mass2) const {
321  return mTrack1->Yt(mass1)-mTrack2->Yt(mass2);
322 }
323 inline float StEStructPairCuts::DeltaEta() const {
324  return mTrack1->Eta()-mTrack2->Eta();
325 }
326 inline float StEStructPairCuts::DeltaEta(float mass1, float mass2) const {
327  return mTrack1->Eta(mass1)-mTrack2->Eta(mass2);
328 }
329 inline float StEStructPairCuts::DeltaPhi() const {
330  return mTrack1->Phi()-mTrack2->Phi();
331 }
332 
333 inline float StEStructPairCuts::DeltaPt() const {
334  return mTrack1->Pt()-mTrack2->Pt();
335 }
336 inline float StEStructPairCuts::mInv() const {
337  return abs(mTrack1->FourMomentum()+mTrack2->FourMomentum());
338 }
339 
340 inline float StEStructPairCuts::kT() const {
341  return 0.5*(mTrack1->FourMomentum()+mTrack2->FourMomentum()).perp();
342 }
343 
344 inline float StEStructPairCuts::qInv() const {
345  return -1.0*(mTrack1->FourMomentum()-mTrack2->FourMomentum()).m();
346 }
347 
348 inline bool StEStructPairCuts::awaySide() {
349  mdeltaPhi=fabs(DeltaPhi());
350  if((mdeltaPhi>M_PI/2.0) && (mdeltaPhi<1.5*M_PI)) return true;
351  return false;
352 }
353 
354 inline bool StEStructPairCuts::sameSide() {
355  mdeltaPhi=fabs(DeltaPhi());
356  if((mdeltaPhi < M_PI/2.0) || (mdeltaPhi> (1.5*M_PI))) return true;
357  return false;
358 }
359 
360 // Both goodDeltaPhi and goodDeltaEta were looking for their deltas to
361 // be smaller than the cut value. To speed up the pairCuts we want
362 // to quickly identify pairs that will not be cut. Also, cuts are actually done
363 // on difference in Z and XY directions. If either of these are too big we don't
364 // want to bother trying rest of pair cuts that we know won't reject pair.
365 inline int StEStructPairCuts::goodDeltaXY() {
366  if ( mGooddeltaZdeltaXYCut &&
367  MidTpcXYSeparation() > mgooddzdxy[1] &&
368  OuterMidTpcXYSeparation() > mgooddzdxy[1] ) {
369  return 1;
370  }
371  return 0;
372 }
373 
374 // For OuterMidTPC we have the complication of tracks exiting via endcap.
375 // OuterMidTpcZSeparation() will return a negative number in this case and
376 // the pair will not be flagged to skip the pair cuts.
377 // I expect a small fraction of pairs are like this so we keep this code
378 // simple and hopefully fast, passing a few more pairs on to detailed
379 // pair cuts than are really necessary.
380 //
381 // I find that the quality cuts can be unexpectedly high (therefore rejecting
382 // the pair) when one track exits via the endcap while the other does not.
383 // I guess we need to do this cut properly after all.
384 inline int StEStructPairCuts::goodDeltaZ() {
385  if ( mGooddeltaZdeltaXYCut &&
386  MidTpcZSeparation() > mgooddzdxy[0] &&
387  OuterMidTpcZSeparation() > mgooddzdxy[0] ) {
388  return 1;
389  }
390  double outerDelZ;
391  if ( mGooddeltaZdeltaXYCut ) {
392  if ((outerDelZ = OuterMidTpcZSeparation()) < 0) {
393  return 1;
394  } else if ((MidTpcZSeparation() > mgooddzdxy[0]) &&
395  (outerDelZ > mgooddzdxy[0] )) {
396  return 1;
397  }
398  }
399  return 0;
400 }
401 
402 inline int StEStructPairCuts::goodDeltaMt(){
403  if(mdeltaMtCut &&
404  ( (mdeltaMt=fabs(DeltaMt())) <mdmt[1] )) return 0;
405  return 1;
406 }
407 
408 
409 inline int StEStructPairCuts::cutDeltaPhi(){
410  if ( mdeltaPhiCut &&
411  ( (mdeltaPhi=DeltaPhi()) <mdphi[0] || mdeltaPhi>mdphi[1])
412  ) return ++(mdphiCounter[mType]);
413  return 0;
414 }
415 
416 inline int StEStructPairCuts::cutDeltaEta(){
417  if(mdeltaEtaCut &&
418  ( (mdeltaEta=DeltaEta()) <mdeta[0] || mdeltaEta>mdeta[1])
419  ) return ++(mdetaCounter[mType]);
420  return 0;
421 }
422 
423 
424 inline int StEStructPairCuts::cutDeltaMt(){
425  if(mdeltaMtCut &&
426  ( (mdeltaMt=DeltaMt()) <mdmt[0] || mdeltaMt>mdmt[1])
427  ) return ++(mdmtCounter[mType]);
428  return 0;
429 }
430 
431 inline int StEStructPairCuts::cutqInvORNominalEntranceSep(){
432  /* small qInv and entrance cut */
433 
434  if( mqInvCut && mEntSepCut &&
435  ( (mqInvariant=qInv()) <mqInv[0] &&
436  ((mEntranceSeparation=NominalTpcEntranceSeparation())<mEntSep[0]))
437  ) return ++(mqInvCounter[mType]);
438  return 0;
439 }
440 
441 inline int StEStructPairCuts::cutqInv(){
442  if( mqInvCut &&
443  ( (mqInvariant=qInv())<mqInv[0]
444  || mqInvariant>mqInv[1]
445  )
446  ) return ++(mqInvCounter[mType]);
447  return 0;
448 }
449 
450 inline int StEStructPairCuts::cutEntranceSep(){
451  if( mEntSepCut &&
452  ( (mEntranceSeparation=NominalTpcEntranceSeparation())<mEntSep[0]
453  || mEntranceSeparation>mEntSep[1]
454  )
455  ) return ++(mEntSepCounter[mType]);
456  return 0;
457 }
458 
459 
460 inline int StEStructPairCuts::cutMidTpcSepUS(){
461  if( mMidTpcSepUSCut ){
462  double x1=mMidTpcSeparationUS=MidTpcSeparation();
463  double x2=NominalTpcEntranceSeparation();
464  double x3=NominalTpcExitSeparation();
465  if(((x1+x2+x3)/3.)>mMidTpcSepUS[1])return 0; // ok, average is large
466  if( x1<x2 && x1<x3 && x1<mMidTpcSepUS[0]) return ++(msplitUSCounter[mType]);
467  }
468  return 0;
469 }
470 
471 inline int StEStructPairCuts::cutMidTpcSepLS(){
472  if( mMidTpcSepLSCut ){
473  double x1=mMidTpcSeparationLS=MidTpcSeparation();
474  double x2=NominalTpcEntranceSeparation();
475  double x3=NominalTpcExitSeparation();
476  if(((x1+x2+x3)/3.)>mMidTpcSepLS[1])return 0; //ok, average is large
477  if( x1<x2 && x1<x3 && x1<mMidTpcSepLS[0]) return ++(msplitLSCounter[mType]); }
478  return 0;
479 }
480 
481 
482 inline int StEStructPairCuts::cutExitSep(){
483  if( mExitSepCut &&
484  ( (mExitSeparation=NominalTpcExitSeparation()) <mExitSep[0])
485  // || mExitSeparation>mExitSep[1] )
486  ) return ++(mExitSepCounter[mType]);
487  return 0;
488 }
489 
490 inline int StEStructPairCuts::cutQuality() {
491  if( mQualityCut &&
492  ((mQualityVal=quality()) <mQuality[0] || mQualityVal>mQuality[1] )) {
493  return ++(mQualityCounter[mType]);
494  }
495  return 0;
496 }
497 
498 inline int StEStructPairCuts::cutHBT(){
499  if(!mHBTCut || mType==1 || mType==3) return 0; // HBT applies only to LS pairs
500  float dpt = fabs(DeltaPt()); // DeltaPt is signed,
501  float deta = fabs(DeltaEta()); // now DeltaEta and DeltaPhi are signed...
502  float dphi = fabs(DeltaPhi()); //
503  if ( deta<mHBT[0] && dphi<mHBT[1] && dpt<mHBT[2]
504  && mTrack1->Pt()<mHBT[3] && mTrack2->Pt()<mHBT[3] )
505  return ++(mHBTCounter[mType]);
506  return 0;
507 }
508 
509 inline int StEStructPairCuts::cutCoulomb(){
510  if(!mCoulombCut) return 0;
511  float dpt = fabs(DeltaPt()); // DeltaPt is signed,
512  float deta = fabs(DeltaEta()); // now DeltaEta and DeltaPhi are signed...
513  float dphi = fabs(DeltaPhi()); //
514  if ( deta<mCoulomb[0] && dphi<mCoulomb[1] && dpt<mCoulomb[2]
515  && mTrack1->Pt()<mCoulomb[3] && mTrack2->Pt()<mCoulomb[3] )
516  return ++(mCoulombCounter[mType]);
517  return 0;
518 }
519 
520 inline int StEStructPairCuts::cutMerging(){
521  if(!mMergingCut) return 0;
522  if ( NominalTpcAvgXYSeparation()<mMerging[0] &&
523  NominalTpcAvgZSeparation() <mMerging[1])
524  return ++(mMergingCounter[mType]);
525  return 0;
526 }
527 
528 inline int StEStructPairCuts::cutCrossing(){
529  // Crossing geometries are a little complicated...
530  // Pair Cut if:
531  // + + dPt & dPhi opp. sign
532  // - - dPt & dPhi same sign
533  // + - dPhi>0
534  // - + dPhi<0
535  // Above table is for positive field, for reversed field switch charge signs
536  // I'm going to explicity list each case rather then try to be clever here.
537 
538  if(!mCrossingCut) return 0; // make the common case fast
539  mretVal = 0;
540  if ( MidTpcXYSeparation()<mCrossing[0] && MidTpcZSeparation()<mCrossing[1]) {
541  float dphi = mTrack1->Phi()-mTrack2->Phi(); // signed DeltaPhi
542  float dpt = mTrack1->Pt()- mTrack2->Pt(); // signed DeltaPt
543  if (mType==1 || mType==3) { // US pair
544  if(mBField>=0) { // pos field
545  if (mTrack1->Charge()>0 && dphi>0) mretVal = 1; // + -
546  if (mTrack1->Charge()<0 && dphi<0) mretVal = 1; // - +
547  } else { // rev field
548  if (mTrack1->Charge()>0 && dphi<0) mretVal = 1; // rev +- : same as -+ above
549  if (mTrack1->Charge()<0 && dphi>0) mretVal = 1; // rev -+ : same as +- above
550  }
551  } else { // LS pair
552  if(mBField>=0) { // pos field
553  if (mTrack1->Charge()>0 && dphi*dpt<0) mretVal = 1; // + +
554  if (mTrack1->Charge()<0 && dphi*dpt>0) mretVal = 1; // - -
555  } else { // rev field
556  if (mTrack1->Charge()>0 && dphi*dpt>0) mretVal = 1; // rev ++ : same as --
557  if (mTrack1->Charge()<0 && dphi*dpt<0) mretVal = 1; // rev -- : same as ++
558  }
559  }
560  }
561 
562  // This should do what that does. Might be faster?
563  // Problem is that although signbit is not zero if the bit is set it is not actually
564  // guaranteed to be 1.
565 // mretVal = 0;
566 // if ( MidTpcXYSeparation()<mCrossing[0] && MidTpcZSeparation()<mCrossing[1]) {
567 // int dphi = 1-2*signbit(mTrack1->Phi()-mTrack2->Phi()); // sign of DeltaPhi
568 // int chrg = 1-2*signbit(mBField*mTrack1->Charge());
569 // if (mType==1 || mType==3) { // US pair
570 // if (chrg == dphi) mretVal = 1;
571 // } else {
572 // int dpt = 1-2*signbit(mTrack1->Pt()- mTrack2->Pt()); // sign of DeltaPt
573 // if (chrg != dpt*dphi) mretVal = 1;
574 // }
575 // }
576  if (mretVal==0) return 0;
577  return ++(mCrossingCounter[mType]);
578 }
579 
580 inline int StEStructPairCuts::cutMerging2() {
581  // If tracks ever get close we cut them. Ignore TPC entrance.
582  // Doesn't seem to take care of crossing tracks.
583  // Turns out when one of the tracks leaves via the endcap the nominal
584  // Z separation is -1. We use only the XY separation in this case.
585  //
586  // Ooops. Looks like that final case of only using the XY separation when
587  // the nominal Z separation is negative allows cutting pairs where one track leaves
588  // the TPC endcap and the other has arbitrary eta. Not what I wanted.
589  // Use difference in radius instead of Z when track leaves via endcap.
590  // Also, if track does not make it to a given radius ignore the cut there.
591  if (!mMergingCut2) {
592  return 0;
593  }
594  // For normal eta and Z-Vertex cuts all tracks will exceed the MidTPC radius before
595  // passing the Z position of the endcap.
596  if ( (MidTpcZSeparation()<mMerging2[0]) && (MidTpcXYSeparation()<mMerging2[1]) ) {
597  return ++(mMergingCounter2[mType]);
598  }
599  // Possible for track to pass through endcap before getting to radius of OuterMidTPC.
600  // If tracks go through different endcaps do not cut.
601  // If only one track goes through an endcap or they go through same endcap then
602  // replace z difference cut with radial difference cut.
603  // Not sure XY separation cut is optimal in endcap case.
604  if ( OuterMidTpcZSeparation()>0 ) {
605  if ( OuterMidTpcZSeparation()<mMerging2[0] && OuterMidTpcXYSeparation()<mMerging2[1] ) {
606  return ++(mMergingCounter2[mType]);
607  }
608  } else if ( OuterMidTpcZSeparation()>-3 ) {
609  if ( (TpcRadialEndCapSeparation()<mMerging2[0]) && (OuterMidTpcXYSeparation()<mMerging2[1]) ) {
610  return ++(mMergingCounter2[mType]);
611  }
612  }
613  if (NominalTpcZExitSeparation() > 0) {
614  if ( (NominalTpcZExitSeparation()<mMerging2[0]) && (NominalTpcXYExitSeparation()<mMerging2[1]) ) {
615  return ++(mMergingCounter2[mType]);
616  }
617  } else if ( NominalTpcZExitSeparation() > -3) {
618  if ( (TpcRadialEndCapSeparation()<mMerging2[0]) && (NominalTpcXYExitSeparation()<mMerging2[1]) ) {
619  return ++(mMergingCounter2[mType]);
620  }
621  }
622  return 0;
623 }
624 inline int StEStructPairCuts::cutCrossing2() {
625  // For tracks that cross.
626  // This removes a slot narrow in Z separation but perhaps wide in XY separation.
627  // I find that low-y_t, high-y_t pairs can be far apart at the middle of the
628  // tpc but close together near the outer radius, so require large z separation
629  // at both radii.
630 
631  if (!mCrossingCut2) {
632  return 0;
633  }
634  // If tracks are far enough apart in Z don't bother to look for crossing.
635  // If one (or both) tracks pass through endcap before OuterMid point
636  // we ignore that check.
637  if ( OuterMidTpcZSeparation()>0 ) {
638  if ( MidTpcZSeparation()>mCrossing2[0] && OuterMidTpcZSeparation()>mCrossing2[1] ) {
639  return 0;
640  }
641  } else if ( OuterMidTpcZSeparation()==-3 ) {
642  return 0;
643  } else {
644  if ( MidTpcZSeparation()>mCrossing2[0] ) {
645  return 0;
646  }
647  }
648 // Following is the old cut. Problem for low pt tracks (which probably makes almost
649 // no difference overall.
650 // if ( (MidTpcZSeparation()>mCrossing2[0]) && (OuterMidTpcZSeparation()>mCrossing2[1]) ) {
651 // return 0;
652 // }
653  if (!TracksCrossInPhi()) {
654  return 0;
655  }
656  return ++(mCrossingCounter2[mType]);
657 }
658 inline int StEStructPairCuts::cutLUT() {
659  // Use a Look Up Table to determine which pairs get close enough
660  // to affect each other.
661 
662  if (!mLUTCut) {
663  return 0;
664  }
665  double delPhi = mTrack1->Phi() - mTrack2->Phi();
666  double delEta = mTrack1->Eta() - mTrack2->Eta();
667  if (mLUT->cut(mTrack1->Curvature(),mTrack2->Curvature(),delPhi,delEta)) {
668  return ++(mLUTCounter[mType]);
669  }
670  return 0;
671 }
672 
673 inline int StEStructPairCuts::cutMass() {
674  if (!(mpionOtherMassCut | mpionpionMassCut | mpionKaonMassCut | mpionprotonMassCut |
675  mKaonOtherMassCut | mKaonKaonMassCut | mKaonprotonMassCut | mprotonOtherMassCut |
676  mprotonprotonMassCut | mOtherOtherMassCut)) {
677  return 0;
678  }
679 
680  int mode[4][4] = {{9, 0, 4, 7}, {0, 1, 2, 3}, {4, 2, 5, 6}, {7, 3, 6, 8}};
681  int it1 = Track1()->PID();
682  int it2 = Track2()->PID();
683  if (it1 < 0 || 3 < it1) {
684  return 0;
685  }
686  if (it2 < 0 || 3 < it2) {
687  return 0;
688  }
689  int iBin = mode[it1][it2];
690 
691  double e, e1, e2, p1, p2, p[3], m, m1, m2;
692  p1 = Track1()->Ptot();
693  p2 = Track2()->Ptot();
694  p[0] = Track1()->Px() + Track2()->Px();
695  p[1] = Track1()->Py() + Track2()->Py();
696  p[2] = Track1()->Pz() + Track2()->Pz();
697  float Mass[] = {0.1396, 0.1396, 0.497, 0.9383};
698  float Mass2[] = {0.01949, 0.01949, 0.247, 0.880};
699  if (9 == iBin) {
700  // For o-o try using m1 = m2 = 0.
701  m1 = 0;
702  m2 = 0;
703  e1 = p1;
704  e2 = p2;
705  } else {
706  m1 = Mass[it1];
707  m2 = Mass[it2];
708  e1 = sqrt(p1*p1 + Mass2[it1]);
709  e2 = sqrt(p2*p2 + Mass2[it2]);
710  }
711  e = e1 + e2;
712  m = sqrt(e*e - p[0]*p[0] - p[1]*p[1] - p[2]*p[2]);
713  // Cut on invariant mass to keep or exclude resonances.
714  switch(iBin) {
715  case 0: { // pion-Other
716  if (mpionOtherMassCut) {
717  if (mpionOtherMass[0] < mpionOtherMass[1]) { // Exclude mass within window
718  if ((mpionOtherMass[0] < m) && (m < mpionOtherMass[1])) {
719  return ++(mpionOtherMassCounter[mType]);
720  }
721  } else { // Keep only pairs within mass windo
722  if ((m < mpionOtherMass[1]) || (mpionOtherMass[0] < m)) {
723  return ++(mpionOtherMassCounter[mType]);
724  }
725  }
726  }
727  break;
728  }
729  case 1: { // pion-pion
730  if (mpionpionMassCut) {
731  if (mpionpionMass[0] < mpionpionMass[1]) { // Exclude mass within window
732  if ((mpionpionMass[0] < m) && (m < mpionpionMass[1])) {
733  return ++(mpionpionMassCounter[mType]);
734  }
735  } else { // Keep only pairs within mass windo
736  if ((m < mpionpionMass[1]) || (mpionpionMass[0] < m)) {
737  return ++(mpionpionMassCounter[mType]);
738  }
739  }
740  }
741  break;
742  }
743  case 2: { // pion-Kaon
744  if (mpionKaonMassCut) {
745  if (mpionKaonMass[0] < mpionKaonMass[1]) { // Exclude mass within window
746  if ((mpionKaonMass[0] < m) && (m < mpionKaonMass[1])) {
747  return ++(mpionKaonMassCounter[mType]);
748  }
749  } else { // Keep only pairs within mass windo
750  if ((m < mpionKaonMass[1]) || (mpionKaonMass[0] < m)) {
751  return ++(mpionKaonMassCounter[mType]);
752  }
753  }
754  }
755  break;
756  }
757  case 3: { // pion-proton
758  if (mpionprotonMassCut) {
759  if (mpionprotonMass[0] < mpionprotonMass[1]) { // Exclude mass within window
760  if ((mpionprotonMass[0] < m) && (m < mpionprotonMass[1])) {
761  return ++(mpionprotonMassCounter[mType]);
762  }
763  } else { // Keep only pairs within mass windo
764  if ((m < mpionprotonMass[1]) || (mpionprotonMass[0] < m)) {
765  return ++(mpionprotonMassCounter[mType]);
766  }
767  }
768  }
769  break;
770  }
771  case 4: { // Kaon-Other
772  if (mKaonOtherMassCut) {
773  if (mKaonOtherMass[0] < mKaonOtherMass[1]) { // Exclude mass within window
774  if ((mKaonOtherMass[0] < m) && (m < mKaonOtherMass[1])) {
775  return ++(mKaonOtherMassCounter[mType]);
776  }
777  } else { // Keep only pairs within mass windo
778  if ((m < mKaonOtherMass[1]) || (mKaonOtherMass[0] < m)) {
779  return ++(mKaonOtherMassCounter[mType]);
780  }
781  }
782  }
783  break;
784  }
785  case 5: { // Kaon-Kaon
786  if (mKaonKaonMassCut) {
787  if (mKaonKaonMass[0] < mKaonKaonMass[1]) { // Exclude mass within window
788  if ((mKaonKaonMass[0] < m) && (m < mKaonKaonMass[1])) {
789  return ++(mKaonKaonMassCounter[mType]);
790  }
791  } else { // Keep only pairs within mass windo
792  if ((m < mKaonKaonMass[1]) || (mKaonKaonMass[0] < m)) {
793  return ++(mKaonKaonMassCounter[mType]);
794  }
795  }
796  }
797  break;
798  }
799  case 6: { // Kaon-proton
800  if (mKaonprotonMassCut) {
801  if (mKaonprotonMass[0] < mKaonprotonMass[1]) { // Exclude mass within window
802  if ((mKaonprotonMass[0] < m) && (m < mKaonprotonMass[1])) {
803  return ++(mKaonprotonMassCounter[mType]);
804  }
805  } else { // Keep only pairs within mass windo
806  if ((m < mKaonprotonMass[1]) || (mKaonprotonMass[0] < m)) {
807  return ++(mKaonprotonMassCounter[mType]);
808  }
809  }
810  }
811  break;
812  }
813  case 7: { // proton-Other
814  if (mprotonOtherMassCut) {
815  if (mprotonOtherMass[0] < mprotonOtherMass[1]) { // Exclude mass within window
816  if ((mprotonOtherMass[0] < m) && (m < mprotonOtherMass[1])) {
817  return ++(mprotonOtherMassCounter[mType]);
818  }
819  } else { // Keep only pairs within mass windo
820  if ((m < mprotonOtherMass[1]) || (mprotonOtherMass[0] < m)) {
821  return ++(mprotonOtherMassCounter[mType]);
822  }
823  }
824  }
825  break;
826  }
827  case 8: { // proton-proton
828  if (mprotonprotonMassCut) {
829  if (mprotonprotonMass[0] < mprotonprotonMass[1]) { // Exclude mass within window
830  if ((mprotonprotonMass[0] < m) && (m < mprotonprotonMass[1])) {
831  return ++(mprotonprotonMassCounter[mType]);
832  }
833  } else { // Keep only pairs within mass windo
834  if ((m < mprotonprotonMass[1]) || (mprotonprotonMass[0] < m)) {
835  return ++(mprotonprotonMassCounter[mType]);
836  }
837  }
838  }
839  break;
840  }
841  case 9: { // Other-Other
842  if (mOtherOtherMassCut) {
843  if (mOtherOtherMass[0] < mOtherOtherMass[1]) { // Exclude mass within window
844  if ((mOtherOtherMass[0] < m) && (m < mOtherOtherMass[1])) {
845  return ++(mOtherOtherMassCounter[mType]);
846  }
847  } else { // Keep only pairs within mass windo
848  if ((m < mOtherOtherMass[1]) || (mOtherOtherMass[0] < m)) {
849  return ++(mOtherOtherMassCounter[mType]);
850  }
851  }
852  }
853  break;
854  }
855  }
856  return 0;
857 }
858 
859 inline int StEStructPairCuts::cutDeltaPhiH(){
860  if(!mdeltaPhiCut) return 0;
861  mretVal=cutDeltaPhi();
862  mvalues[mdphiName.idx]=mdeltaPhi;
863  return mretVal;
864 }
865 
866 inline int StEStructPairCuts::cutDeltaEtaH(){
867  if(!mdeltaEtaCut){
868  mdeltaEta=DeltaEta(); // still needed for other cuts
869  return 0;
870  }
871  mretVal=cutDeltaEta();
872  mvalues[mdetaName.idx]=mdeltaEta;
873  return mretVal;
874 }
875 
876 inline int StEStructPairCuts::cutDeltaMtH(){
877  if(!mdeltaMtCut) return 0;
878  mretVal=cutDeltaMt();
879  mvalues[mdmtName.idx]=mdeltaMt;
880  return mretVal;
881 }
882 
883 inline int StEStructPairCuts::cutqInvH(){
884  if(!mqInvCut) return 0;
885  mretVal=cutqInv();
886  mvalues[mqInvName.idx]=mqInvariant;
887  return mretVal;
888 }
889 
890 inline int StEStructPairCuts::cutEntranceSepH(){
891  if(!mEntSepCut) return 0;
892  mretVal=cutEntranceSep();
893  mvalues[mEntSepName.idx]=mEntranceSeparation;
894  return mretVal;
895 }
896 
897 inline int StEStructPairCuts::cutMidTpcSepLSH(){
898  if(!mMidTpcSepLSCut) return 0;
899  mretVal=cutMidTpcSepLS();
900  mvalues[mMidTpcSepLSName.idx]=mMidTpcSeparationLS;
901  return mretVal;
902 }
903 
904 inline int StEStructPairCuts::cutMidTpcSepUSH(){
905  if(!mMidTpcSepUSCut) return 0;
906  mretVal=cutMidTpcSepUS();
907  mvalues[mMidTpcSepUSName.idx]=mMidTpcSeparationUS;
908  return mretVal;
909 }
910 
911 
912 inline int StEStructPairCuts::cutExitSepH(){
913  if(!mExitSepCut) return 0;
914  mretVal=cutExitSep();
915  mvalues[mExitSepName.idx]=mExitSeparation;
916  return mretVal;
917 }
918 
919 inline int StEStructPairCuts::cutQualityH(){
920  if(!mQualityCut) return 0;
921  mretVal=cutQuality();
922  mvalues[mQualityName.idx]=mQualityVal;
923  return mretVal;
924 }
925 
926 
927 inline int StEStructPairCuts::correlationDepth(){
928 
929  unsigned int am1=mTrack1->TopologyMapData(0)&65535;
930  unsigned int am2=mTrack2->TopologyMapData(0)&65535;
931  if(am1==0 || am2==0) return 0;
932  if(am1==am2) return 1;
933 
934  unsigned int msk=(65535<<16);
935 
936  unsigned int bm1=mTrack1->TopologyMapData(0)&msk;
937  unsigned int bm2=mTrack2->TopologyMapData(0)&msk;
938  bm1=(16>>bm1);
939  bm2=(16>>bm2);
940 
941  if(am1==bm2 || bm1==am2) return 2;
942  if(bm1==bm2){
943  if(bm1==0) return 0;
944  return 3;
945  }
946 
947  unsigned int cm1=mTrack1->TopologyMapData(1)&65535;
948  unsigned int cm2=mTrack2->TopologyMapData(1)&65535;
949  if(cm1==cm2 && cm1==0) return 0;
950 
951  if(am1==cm2 || am2==cm1) return 4;
952  if( ((bm1==cm2) && bm1!=0) || ((bm2==cm1)&&(bm2!=0)) ) return 4;
953  if( cm1==cm2) return 5;
954 
955  unsigned int dm1=mTrack1->TopologyMapData(1)&msk;
956  unsigned int dm2=mTrack2->TopologyMapData(1)&msk;
957  dm1=(16>>dm1);
958  dm2=(16>>dm2);
959  if(dm1==dm2 && dm1==0) return 0;
960 
961  if(am1==dm2 || am2==dm1) return 6;
962  if( ((bm1==dm2) && bm1!=0) || ((bm2==dm1)&&(bm2!=0)) ) return 6;
963  if( ((cm1==dm2) && cm1!=0) || ((cm2==dm1)&&(cm2!=0)) ) return 6;
964  if( dm1==dm2 && dm1!=0) return 7;
965 
966  return 0;
967 }
968 
969 
970 
971 #endif
972 
973 /***********************************************************************
974  *
975  * $Log: StEStructPairCuts.h,v $
976  * Revision 1.21 2012/11/16 21:22:27 prindle
977  * 2ptCorrelations: SS, AS histograms. Get eta limits from cuts. Fit PtAll histogram. Add histograms to keep track of eta, phi limits. A few more histograms
978  * Binning: Add quality cut.
979  * CutBin: modify mode9
980  * PairCuts: modify goodDeltaZ for case of one track leaving via endcap.
981  *
982  * Revision 1.20 2010/09/02 21:24:08 prindle
983  * 2ptCorrelations: Fill histograms for event mixing information
984  * Option for common mixing buffer
985  * Switch to selectively fill QInv histograms (which take a long time)
986  * CutBin: Moved PID code to Track class from Pair class. Needed to update this code.
987  * PairCuts: Moved PID code from here to Track class.
988  * Removed unnecessary creation of StThreeVector which seem to take a long time
989  * Add ToF momentum cuts, modify dEdx momentum cuts. (Now allow dEdx to be
990  * be identified up to 15GeV/c, ToF up to 10GeV/c.)
991  *
992  * Revision 1.19 2010/03/02 21:45:28 prindle
993  * Had a problem with pair cuts when one track exited via endplate
994  * Calculate maxDEta properly
995  * Warning if you try turning histograms for pair cuts on
996  *
997  * Revision 1.18 2009/11/09 21:32:41 prindle
998  * Fix warnings about casting char * to a const char * by redeclaring as const char *.
999  *
1000  * Revision 1.17 2009/05/08 00:09:55 prindle
1001  * In 2ptCorrelations we added switches to select blocks of histograms to fill.
1002  * (See constructor in StEStruct2ptCorrelations.cxx)
1003  * Use a brute force method for checking crossing cuts. I had too many corner
1004  * cases with my clever check.
1005  * In Binning, change Yt limit and add methods for accessing number of histogram bins
1006  * to use (used in Support)
1007  *
1008  * Revision 1.16 2008/12/02 23:45:07 prindle
1009  * Changed switchYt to switchXX (etc.) to better reflect function.
1010  * Change minYt to 1.0 in Binning so YtYt histogram doesn't have empty lower bin (pt = 0.164 for yt = 1.0)
1011  * In CutBin: remove initPtBin
1012  * add mode 8
1013  * add notSymmetrized (used in Support)
1014  * Added LUT (Look Up Table) for pair cuts. Experimental for now.
1015  * Modified cutMerging2 (to look at track separation at a few radii)
1016  * and cutCrossing2 so it doesn't accidentally reject almost back to back tracks.
1017  *
1018  * Revision 1.15 2008/05/01 23:39:14 prindle
1019  * I was using a triangular region for the merging cut. Decided a rectangular
1020  * region was safer.
1021  *
1022  * Revision 1.14 2008/03/19 22:06:01 prindle
1023  * Added doInvariantMass flag.
1024  * Added some plots in pairDensityHistograms.
1025  * SetZOffset used to only be done when doPairDensity was true.
1026  * Moved creating/copying pairDensity histograms to same place as other histograms.
1027  * Added cutBinHistMode
1028  * mode3 neck was defined as yt1<2.2 && yt2<2.2 (and not soft)
1029  * now is 1.8<yt1<2.2 && 1.8<yt2<2.2
1030  * Added gooddzdxy, Merging2 and Crossing2 to pair cuts.
1031  *
1032  * Revision 1.13 2007/11/26 19:55:25 prindle
1033  * In 2ptCorrelations: Support for keeping all z-bins of selected centralities
1034  * Change way \hat{p_t} is calculated for parent distributions in pid case.
1035  * Binning Added parent binning (for \hat{p_t}
1036  * CutBin: Mode 5 extensively modified.
1037  * Added invariant mass cuts (probably a bad idea in general.)
1038  *
1039  * Revision 1.12 2007/01/26 17:17:11 msd
1040  * Implemented new binning scheme: dEta stored in array with bin centered at zero, dPhi array has bins centered at zero and pi. Final DEtaDPhi has 25x25 bins with dPhi bin width of pi/12 so all major angles are centered in bins.
1041  *
1042  * Revision 1.11 2006/10/02 22:21:04 prindle
1043  * Store only quadrant of eta_Delta - phi_Delta array/histogram.
1044  * Store half of eta_Sigma - phi_Delta array/histogram.
1045  * This required modifications in Binning.
1046  * I had a bug in the pair loop (which left +- not fully symmetrized)
1047  * and had to make changes in cut bins for mode 5 (and 3 I think)
1048  * when I fixed this.
1049  * Also change crossing cut to use only two parameters, the sign of
1050  * the magnetic field being taken from the MuDst.
1051  *
1052  * Revision 1.10 2006/04/10 23:42:32 porter
1053  * Added sameSide() & awaySide() methods to PairCut (so only defined in 1 place)
1054  * and added the eta_delta weighting as a binned correctin defined by the eta-limits in
1055  * the StEStructBinning object
1056  *
1057  * Revision 1.9 2006/04/06 01:01:22 prindle
1058  *
1059  * New mode in CutBin, 5, to do pid correlations. There is still an issue
1060  * of how to set the pt ranges allowed for the different particle types.
1061  * For data we probably want to restrict p to below 1GeV for pi and K, but
1062  * for Hijing and Pythia we can have perfect pid. Currently cuts are type
1063  * into the code (so you have to re-compile to change them.)
1064  *
1065  * In the Correlations code I split -+ from +- and am keeping track of
1066  * pt for each cut bin. These required changes in the Support code.
1067  *
1068  * Revision 1.8 2006/04/04 22:10:13 porter
1069  * a handful of changes (specific to correlations)
1070  * - added StEStructQAHists so that if NOT input frm Maker, each analysis has its own
1071  * - used ability to get any max,min val from the cut class - or z-vertex binning
1072  * - put z-vertex binning into 1 place
1073  * - switched back 1st line of pair cut method to keep pair if good, not to reject if bad.
1074  * - Pair cut object is now pointer in correlations
1075  * - some diagnostic printouts available from macro
1076  * - Duncan's delta-phi binning change
1077  *
1078  * Revision 1.7 2006/02/22 22:05:19 prindle
1079  * Removed all references to multRef (?)
1080  * Added cut mode 5 for particle identified correlations.
1081  * Other cut modes should be same as before
1082  *
1083  * Revision 1.6 2005/09/14 17:14:25 msd
1084  * Large update, added new pair-cut system, added pair density plots for new analysis mode (4), added event mixing cuts (rewrote buffer for this)
1085  *
1086  * Revision 1.5 2005/03/03 01:30:44 porter
1087  * updated StEStruct2ptCorrelations to include pt-correlations and removed
1088  * old version of pt-correlations from chunhuih (StEStruct2ptPtNbar)
1089  *
1090  * Revision 1.4 2004/06/25 03:11:50 porter
1091  * New cut-binning implementation and modified pair-cuts for chunhui to review
1092  *
1093  * Revision 1.3 2004/06/16 20:00:43 chunhuih
1094  *
1095  * changed one more post-increment operator to pre-increment operator.
1096  *
1097  * Revision 1.2 2004/03/19 19:07:43 chunhuih
1098  *
1099  * Use pre-increment instead of post-increment operators for the return values
1100  * of a set of cut methods. This returns the correct cut value for the first
1101  * time the counting array is incremented.
1102  *
1103  * Revision 1.1 2003/10/15 18:20:46 porter
1104  * initial check in of Estruct Analysis maker codes.
1105  *
1106  *
1107  *********************************************************************/
1108 
1109 
1110 
1111 
1112 
1113 
1114 
StEStructPairCuts::mcutMode
int mcutMode
just a dummy holder used a lot
Definition: StEStructPairCuts.h:121
Generated by   doxygen 1.8.5