StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructPairCuts.cxx
1 /**********************************************************************
2  *
3  * $Id: StEStructPairCuts.cxx,v 1.15 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 #include "StEStructPairCuts.h"
14 #include <math.h>
15 #include <stdlib.h>
16 
17 ClassImp(StEStructPairCuts)
18 
19 StEStructPairCuts::StEStructPairCuts(): StEStructCuts(){ init(); };
20 StEStructPairCuts::StEStructPairCuts(const char* cutfileName): StEStructCuts(cutfileName) { init(); };
21 
22 StEStructPairCuts::~StEStructPairCuts() {
23  delete mLUT;
24 };
25 
26 void StEStructPairCuts::init(){
27 
28  strcpy(mcutTypeName,"Pair");
29  initCuts();
30  initNames();
31  mLUT = new StEStructPairLUT();
32  if(isLoaded())loadCuts();
33 
34 }
35 
36 void StEStructPairCuts::initCuts(){
37 
38  mdphi[0]=mdphi[1]=0;
39  mdeta[0]=mdeta[1]=0;
40  mgooddzdxy[0]=mgooddzdxy[1]=0;
41  mdmt[0]=mdmt[1]=0;
42  mqInv[0]=mqInv[1]=0;
43  mEntSep[0]=mEntSep[1]=0;
44  mExitSep[0]=mExitSep[1]=0;
45  mQuality[0]=mQuality[1]=0;
46  mMidTpcSepLS[0]=mMidTpcSepLS[1]=0;
47  mMidTpcSepUS[0]=mMidTpcSepUS[1]=0;
48  mHBT[0]=mHBT[1]=mHBT[2]=mHBT[3]=0;
49  mCoulomb[0]=mCoulomb[1]=mCoulomb[2]=mCoulomb[3]=0;
50  mMerging[0]=mMerging[1]=0;
51  mMerging2[0]=mMerging2[1]=0;
52  mCrossing[0]=mCrossing[1]=0;
53  mCrossing2[0]=mCrossing2[1]=0;
54  mLUTParams[0]=mLUTParams[1]=0;
55 
56  // Would be nice to have some way to change these without recompiling.
57  // Normally upper pt cuts are 999, 1.0, 1.0, 1.5 for all, pi, K, p
58  // Want to run Pythia without upper cut (I think we actually use 20 in cuts file.)
59  // I don't know right off what reasonable upper limits are on pid.
60  // dEdx can be done in the relativistic rise region, but probably Kaons can't be distinguished.
61  // With a requirement track is 2\sigma from any other particle we are probably ok without upper limits on data.
62  // For comparison with Pythia, Hijing and Therminator we need to think about it.
63  mdEdxMomentumCut[0][0] = 0.0;
64  mdEdxMomentumCut[0][1] = 999.0;
65  mdEdxMomentumCut[1][0] = 0.1; // pi
66  mdEdxMomentumCut[1][1] = 15.0;
67  mdEdxMomentumCut[2][0] = 0.1; // K
68  mdEdxMomentumCut[2][1] = 15.0;
69  mdEdxMomentumCut[3][0] = 0.2; // p
70  mdEdxMomentumCut[3][1] = 15.0;
71 
72  mToFMomentumCut[0][0] = 0.0;
73  mToFMomentumCut[0][1] = 999.0;
74  mToFMomentumCut[1][0] = 0.1; // pi
75  mToFMomentumCut[1][1] = 10.0;
76  mToFMomentumCut[2][0] = 0.1; // K
77  mToFMomentumCut[2][1] = 10.0;
78  mToFMomentumCut[3][0] = 0.2; // p
79  mToFMomentumCut[3][1] = 10.0;
80 
81  mdeltaPhiCut=mdeltaEtaCut=mGooddeltaZdeltaXYCut=mdeltaMtCut=mqInvCut=mEntSepCut=mExitSepCut=mQualityCut=mMidTpcSepLSCut=mMidTpcSepUSCut=false;
82  mHBTCut=mCoulombCut=mMergingCut=mCrossingCut=mMergingCut2=mCrossingCut2=mLUTCut = false;
83  mpionMomentumCut=mKaonMomentumCut=mprotonMomentumCut = false;
84  mpionOtherMassCut=mpionpionMassCut=mpionKaonMassCut=mpionprotonMassCut = false;
85  mKaonOtherMassCut=mKaonKaonMassCut=mKaonprotonMassCut=mprotonOtherMassCut = false;
86  mprotonprotonMassCut=mOtherOtherMassCut = false;
87 
88 
89  for(int i=0;i<4;i++) {
90  mdphiCounter[i]=mdetaCounter[i]=mgooddzdxyCounter[i]=mdmtCounter[i]=mqInvCounter[i]=mEntSepCounter[i]=0;
91  mExitSepCounter[i]=mQualityCounter[i]=msplitLSCounter[i]=msplitUSCounter[i]=0;
92  mHBTCounter[i]=mCoulombCounter[i]=mMergingCounter[i]=mCrossingCounter[i]=mMergingCounter2[i]=mCrossingCounter2[i]=mLUTCounter[i]=0;
93  mpionMomentumCounter[i]=mKaonMomentumCounter[i]=mprotonMomentumCounter[i]=0;
94  mpionOtherMassCounter[i]=mpionpionMassCounter[i]=mpionKaonMassCounter[i]=mpionprotonMassCounter[i]=0;
95  mKaonOtherMassCounter[i]=mKaonKaonMassCounter[i]=mKaonprotonMassCounter[i]=mprotonOtherMassCounter[i]=0;
96  mprotonprotonMassCounter[i]=mOtherOtherMassCounter[i]=0;
97  }
98  mdeltaPhi=mdeltaEta=mdeltaMt=mqInvariant= mEntranceSeparation=mExitSeparation=mQualityVal=mMidTpcSeparationLS=mMidTpcSeparationUS=0;
99 
100  mapMask0 = 0xFFFFFF00;
101  mapMask1 = 0x1FFFFF;
102  for(int i=0;i<32;i++)bitI[i]=1UL<<i;
103 
104  mZoffset = 0;
105 
106 };
107 
108 void StEStructPairCuts::initNames(){
109 
110  strcpy(mdphiName.name,"DeltaPhi");
111  strcpy(mdetaName.name,"DeltaEta");
112  strcpy(mgooddzdxyName.name,"GoodDeltaZDeltaXY");
113  strcpy(mdmtName.name,"DeltaMt");
114  strcpy(mqInvName.name,"qInv");
115  strcpy(mEntSepName.name,"EntranceSep");
116  strcpy(mExitSepName.name,"ExitSep");
117  strcpy(mQualityName.name,"Quality");
118  strcpy(mMidTpcSepLSName.name,"MidTpcSepLikeSign");
119  strcpy(mMidTpcSepUSName.name,"MidTpcSepUnlikeSign");
120  strcpy(mHBTName.name,"HBT");
121  strcpy(mCoulombName.name,"Coulomb");
122  strcpy(mMergingName.name,"Merging");
123  strcpy(mMergingName2.name,"Merging2");
124  strcpy(mCrossingName.name,"Crossing");
125  strcpy(mCrossingName2.name,"Crossing2");
126  strcpy(mLUTName.name,"LUT");
127  strcpy(mpionMomentumName.name,"pionMomentumRange");
128  strcpy(mKaonMomentumName.name,"KaonMomentumRange");
129  strcpy(mprotonMomentumName.name,"protonMomentumRange");
130 
131  strcpy(mpionOtherMassName.name,"pionOtherMass");
132  strcpy(mpionpionMassName.name,"pionpionMass");
133  strcpy(mpionKaonMassName.name,"pionKaonMass");
134  strcpy(mpionprotonMassName.name,"pionprotonMass");
135  strcpy(mKaonOtherMassName.name,"KaonOtherMass");
136  strcpy(mKaonKaonMassName.name,"KaonKaonMass");
137  strcpy(mKaonprotonMassName.name,"KaonprotonMass");
138  strcpy(mprotonOtherMassName.name,"protonOtherMass");
139  strcpy(mprotonprotonMassName.name,"protonprotonMass");
140  strcpy(mOtherOtherMassName.name,"OtherOtherMass");
141 
142 };
143 
144 
145 bool StEStructPairCuts::loadBaseCuts(const char* name, const char** vals, int nvals){
146 
147  // cout<<" Cut Name="<<name<<endl;
148  if(!strcmp(name,mdphiName.name)){
149  mdphi[0]=(float)(M_PI*atof(vals[0])); mdphi[1]=(float)(M_PI*atof(vals[1]));
150  mdphiName.idx=createCutHists(name,mdphi);
151  mdeltaPhiCut=true;
152  return true;
153  }
154 
155  if(!strcmp(name,mdetaName.name)){
156  mdeta[0]=atof(vals[0]); mdeta[1]=atof(vals[1]);
157  mdetaName.idx=createCutHists(name,mdeta);
158  mdeltaEtaCut=true;
159  return true;
160  }
161 
162  if(!strcmp(name,mgooddzdxyName.name)){
163  mgooddzdxy[0]=atof(vals[0]); mgooddzdxy[1]=atof(vals[1]);
164  mgooddzdxyName.idx=createCutHists(name,mgooddzdxy);
165  mGooddeltaZdeltaXYCut=true;
166  return true;
167  }
168 
169  if(!strcmp(name,mdmtName.name)){
170  mdmt[0]=atof(vals[0]); mdmt[1]=atof(vals[1]);
171  mdmtName.idx=createCutHists(name,mdmt);
172  mdeltaMtCut=true;
173  return true;
174  }
175 
176  if(!strcmp(name,mqInvName.name)){
177  mqInv[0]=atof(vals[0]); mqInv[1]=atof(vals[1]);
178  mqInvName.idx=createCutHists(name,mqInv);
179  mqInvCut=true;
180  return true;
181  }
182 
183  if(!strcmp(name,mEntSepName.name)){
184  mEntSep[0]=atof(vals[0]); mEntSep[1]=atof(vals[1]);
185  mEntSepName.idx=createCutHists(name,mEntSep);
186  mEntSepCut=true;
187  return true;
188  }
189 
190  if(!strcmp(name,mExitSepName.name)){
191  mExitSep[0]=atof(vals[0]); mExitSep[1]=atof(vals[1]);
192  mExitSepName.idx=createCutHists(name,mExitSep);
193  mExitSepCut=true;
194  return true;
195  }
196 
197  if(!strcmp(name,mQualityName.name)){
198  mQuality[0]=atof(vals[0]); mQuality[1]=atof(vals[1]);
199  mQualityName.idx=createCutHists(name,mQuality);
200  mQualityCut=true;
201  return true;
202  }
203 
204  if(!strcmp(name,mMidTpcSepLSName.name)){
205  mMidTpcSepLS[0]=atof(vals[0]); mMidTpcSepLS[1]=atof(vals[1]);
206  mMidTpcSepLSName.idx=createCutHists(name,mMidTpcSepLS);
207  mMidTpcSepLSCut=true;
208  return true;
209  }
210 
211  if(!strcmp(name,mMidTpcSepUSName.name)){
212  mMidTpcSepUS[0]=atof(vals[0]); mMidTpcSepUS[1]=atof(vals[1]);
213  mMidTpcSepUSName.idx=createCutHists(name,mMidTpcSepUS);
214  mMidTpcSepUSCut=true;
215  return true;
216  }
217 
218  if(!strcmp(name,mHBTName.name)){
219  mHBT[0]=atof(vals[0]); mHBT[1]=atof(vals[1]);
220  mHBT[2]=atof(vals[2]); mHBT[3]=atof(vals[3]);
221  //mHBTName.idx=createCutHists(name,mHBT,4); // not making cut histograms
222  mHBTCut=true;
223  cout << " Loading HBT cut with range of cuts = "<<mHBT[0]<<","<<mHBT[1]<<","<<mHBT[2]<<","<<mHBT[3]<<endl;
224  return true;
225  }
226 
227  if(!strcmp(name,mCoulombName.name)){
228  mCoulomb[0]=atof(vals[0]); mCoulomb[1]=atof(vals[1]);
229  mCoulomb[2]=atof(vals[2]); mCoulomb[3]=atof(vals[3]);
230  //mCoulombName.idx=createCutHists(name,mCoulomb,4); // not making cut histograms
231  mCoulombCut=true;
232  cout << " Loading Coulomb cut with range of cuts = "<<mCoulomb[0]<<","<<mCoulomb[1]<<","<<mCoulomb[2]<<","<<mCoulomb[3]<<endl;
233  return true;
234  }
235 
236  if(!strcmp(name,mMergingName.name)){
237  mMerging[0]=atof(vals[0]); mMerging[1]=atof(vals[1]);
238  //mMergingName.idx=createCutHists(name,mMerging,2); // not making cut histograms
239  mMergingCut=true;
240  cout << " Loading Merging cut with range of cuts = "<<mMerging[0]<<","<<mMerging[1]<<endl;
241  return true;
242  }
243 
244  if(!strcmp(name,mMergingName2.name)){
245  mMerging2[0]=atof(vals[0]); mMerging2[1]=atof(vals[1]);
246  //mMergingName2.idx=createCutHists(name,mMerging2,2); // not making cut histograms
247  mMergingCut2=true;
248  cout << " Loading Merging2 cut with range of cuts = "<<mMerging2[0]<<","<<mMerging2[1]<<endl;
249  return true;
250  }
251 
252  if(!strcmp(name,mCrossingName.name)){
253  mCrossing[0]=atof(vals[0]); mCrossing[1]=atof(vals[1]);
254  //mCrossingName.idx=createCutHists(name,mCrossing,2); // not making cut histograms
255  mCrossingCut=true;
256  cout << " Loading Crossing cut with range of cuts = "<<mCrossing[0]<<","<<mCrossing[1]<<endl;
257  return true;
258  }
259 
260  if(!strcmp(name,mCrossingName2.name)){
261  mCrossing2[0]=atof(vals[0]); mCrossing2[1]=atof(vals[1]);
262  //mCrossingName2.idx=createCutHists(name,mCrossing2,2); // not making cut histograms
263  mCrossingCut2=true;
264  cout << " Loading Crossing2 cut with range of cuts = "<<mCrossing2[0]<<","<<mCrossing2[1]<<endl;
265  return true;
266  }
267  if(!strcmp(name,mLUTName.name)){
268  mLUTParams[0]=atof(vals[0]); mLUTParams[1]=atof(vals[1]);
269  mLUTCut=true;
270  cout << " Calculating Look Up Table with dXY cut = "<<mLUTParams[0]<<" and dZ cut = "<<mLUTParams[1]<< " (takes a little while)" <<endl;
271  mLUT->mDelXYCut = mLUTParams[0];
272  mLUT->mDelZCut = mLUTParams[1];
273  mLUT->initHists();
274  mLUT->fillLUTs();
275  cout << " LUT has been calculated " <<endl;
276  return true;
277  }
278 
279  if(!strcmp(name,mpionMomentumName.name)){
280  mdEdxMomentumCut[1][0]=atof(vals[0]); mdEdxMomentumCut[1][1]=atof(vals[1]);
281  mToFMomentumCut[1][0]=atof(vals[0]); mToFMomentumCut[1][1]=atof(vals[1]);
282  mpionMomentumCut=true;
283  cout << " Loading pion momentum cut with range of values = "<<mdEdxMomentumCut[1][0]<<","<<mdEdxMomentumCut[1][1]<<endl;
284  return true;
285  }
286 
287  if(!strcmp(name,mKaonMomentumName.name)){
288  mdEdxMomentumCut[2][0]=atof(vals[0]); mdEdxMomentumCut[2][1]=atof(vals[1]);
289  mToFMomentumCut[2][0]=atof(vals[0]); mToFMomentumCut[2][1]=atof(vals[1]);
290  mKaonMomentumCut=true;
291  cout << " Loading Kaon momentum cut with range of values = "<<mdEdxMomentumCut[2][0]<<","<<mdEdxMomentumCut[2][1]<<endl;
292  return true;
293  }
294 
295  if(!strcmp(name,mprotonMomentumName.name)){
296  mdEdxMomentumCut[3][0]=atof(vals[0]); mdEdxMomentumCut[3][1]=atof(vals[1]);
297  mToFMomentumCut[3][0]=atof(vals[0]); mToFMomentumCut[3][1]=atof(vals[1]);
298  mprotonMomentumCut=true;
299  cout << " Loading proton momentum cut with range of values = "<<mdEdxMomentumCut[3][0]<<","<<mdEdxMomentumCut[3][1]<<endl; return true;
300  }
301 
302  if(!strcmp(name,mpionOtherMassName.name)){
303  mpionOtherMass[0]=atof(vals[0]); mpionOtherMass[1]=atof(vals[1]);
304  mpionOtherMassCut=true;
305  cout << " Loading pion-Other mass cut with range of values = "<<mpionOtherMass[0]<<","<<mpionOtherMass[1]<<endl; return true;
306  }
307  if(!strcmp(name,mpionpionMassName.name)){
308  mpionpionMass[0]=atof(vals[0]); mpionpionMass[1]=atof(vals[1]);
309  mpionpionMassCut=true;
310  cout << " Loading pion-pion mass cut with range of values = "<<mpionpionMass[0]<<","<<mpionpionMass[1]<<endl; return true;
311  }
312  if(!strcmp(name,mpionKaonMassName.name)){
313  mpionKaonMass[0]=atof(vals[0]); mpionKaonMass[1]=atof(vals[1]);
314  mpionKaonMassCut=true;
315  cout << " Loading pion-Kaon mass cut with range of values = "<<mpionKaonMass[0]<<","<<mpionKaonMass[1]<<endl; return true;
316  }
317  if(!strcmp(name,mpionprotonMassName.name)){
318  mpionprotonMass[0]=atof(vals[0]); mpionprotonMass[1]=atof(vals[1]);
319  mpionprotonMassCut=true;
320  cout << " Loading pion-proton mass cut with range of values = "<<mpionprotonMass[0]<<","<<mpionprotonMass[1]<<endl; return true;
321  }
322  if(!strcmp(name,mKaonOtherMassName.name)){
323  mKaonOtherMass[0]=atof(vals[0]); mKaonOtherMass[1]=atof(vals[1]);
324  mKaonOtherMassCut=true;
325  cout << " Loading Kaon-Other mass cut with range of values = "<<mKaonOtherMass[0]<<","<<mKaonOtherMass[1]<<endl; return true;
326  }
327  if(!strcmp(name,mKaonKaonMassName.name)){
328  mKaonKaonMass[0]=atof(vals[0]); mKaonKaonMass[1]=atof(vals[1]);
329  mKaonKaonMassCut=true;
330  cout << " Loading Kaon-Kaon mass cut with range of values = "<<mKaonKaonMass[0]<<","<<mKaonKaonMass[1]<<endl; return true;
331  }
332  if(!strcmp(name,mKaonprotonMassName.name)){
333  mKaonprotonMass[0]=atof(vals[0]); mKaonprotonMass[1]=atof(vals[1]);
334  mKaonprotonMassCut=true;
335  cout << " Loading Kaon-proton mass cut with range of values = "<<mKaonprotonMass[0]<<","<<mKaonprotonMass[1]<<endl; return true;
336  }
337  if(!strcmp(name,mprotonOtherMassName.name)){
338  mprotonOtherMass[0]=atof(vals[0]); mprotonOtherMass[1]=atof(vals[1]);
339  mprotonOtherMassCut=true;
340  cout << " Loading proton-Other mass cut with range of values = "<<mprotonOtherMass[0]<<","<<mprotonOtherMass[1]<<endl; return true;
341  }
342  if(!strcmp(name,mprotonprotonMassName.name)){
343  mprotonprotonMass[0]=atof(vals[0]); mprotonprotonMass[1]=atof(vals[1]);
344  mprotonprotonMassCut=true;
345  cout << " Loading proton-proton mass cut with range of values = "<<mprotonprotonMass[0]<<","<<mprotonprotonMass[1]<<endl; return true;
346  }
347  if(!strcmp(name,mOtherOtherMassName.name)){
348  mOtherOtherMass[0]=atof(vals[0]); mOtherOtherMass[1]=atof(vals[1]);
349  mOtherOtherMassCut=true;
350  cout << " Loading Other-Other mass cut with range of values = "<<mOtherOtherMass[0]<<","<<mOtherOtherMass[1]<<endl; return true;
351  }
352 
353  // cout<<" didn't find any cut with this name "<<endl;
354  return false;
355 };
356 
357 void StEStructPairCuts::printCutCounts(ostream& ofs, const char* cutName,int c1,int c2){
358  ofs<<cutName<<c1<<" + "<<c2<<" = "<<c1+c2<<endl;
359 }
360 
361 
362 void StEStructPairCuts::printCutStats(ostream& ofs){
363 
364  // ofs<<"# ******************************************** "<<endl;
365  // ofs<<"# *************** Pair Cuts ****************** "<<endl;
366  // ofs<<"# *** format = Cut, minvalue, maxvalue *** "<<endl;
367  // ofs<<"# *** Sib LS + US = Total *** "<<endl;
368  // ofs<<"# *** Mix LS + US = Total *** "<<endl;
369  // ofs<<"# ******************************************** "<<endl;
370  ofs<<endl;
371  const char* cutTypes[]={"#--- Sibling Pairs : LS + US = ",
372  "#--- Mixed Pairs : LS + US = "};
373  if(mdeltaPhiCut){
374  ofs<<mdphiName.name<<","<<mdphi[0]/M_PI<<","<<mdphi[1]/M_PI<<"\t\t\t"<<" # pair dphi cut"<<endl;
375  printCutCounts(ofs,cutTypes[0],mdphiCounter[0],mdphiCounter[1]);
376  printCutCounts(ofs,cutTypes[1],mdphiCounter[2],mdphiCounter[3]);
377  }
378  if(mdeltaEtaCut){
379  ofs<<mdetaName.name<<","<<mdeta[0]<<","<<mdeta[1]<<"\t\t\t"<<" # pair deta cut"<<endl;
380  printCutCounts(ofs,cutTypes[0],mdetaCounter[0],mdetaCounter[1]);
381  printCutCounts(ofs,cutTypes[1],mdetaCounter[2],mdetaCounter[3]);
382  }
383  if(mGooddeltaZdeltaXYCut){
384  ofs<<mgooddzdxyName.name<<","<<mgooddzdxy[0]<<","<<mgooddzdxy[1]<<"\t\t\t"<<" # pair deta cut"<<endl;
385  printCutCounts(ofs,cutTypes[0],mgooddzdxyCounter[0],mgooddzdxyCounter[1]);
386  printCutCounts(ofs,cutTypes[1],mgooddzdxyCounter[2],mgooddzdxyCounter[3]);
387  }
388 
389  if(mdeltaMtCut){
390  ofs<<mdmtName.name<<","<<mdmt[0]<<","<<mdmt[1]<<"\t\t\t"<<" # pair dmt cut"<<endl;
391  printCutCounts(ofs,cutTypes[0],mdmtCounter[0],mdmtCounter[1]);
392  printCutCounts(ofs,cutTypes[1],mdmtCounter[2],mdmtCounter[3]);
393  }
394 
395  if(mqInvCut){
396  ofs<<mqInvName.name<<","<<mqInv[0]<<","<<mqInv[1]<<"\t\t\t"<<" # pair qInv cut"<<endl;
397  printCutCounts(ofs,cutTypes[0],mqInvCounter[0],mqInvCounter[1]);
398  printCutCounts(ofs,cutTypes[1],mqInvCounter[2],mqInvCounter[3]);
399  }
400 
401  if(mEntSepCut){
402  ofs<<mEntSepName.name<<","<<mEntSep[0]<<","<<mEntSep[1]<<"\t\t\t"<<" # pair EntSep cut"<<endl;
403  printCutCounts(ofs,cutTypes[0],mEntSepCounter[0],mEntSepCounter[1]);
404  printCutCounts(ofs,cutTypes[1],mEntSepCounter[2],mEntSepCounter[3]);
405  }
406 
407  if(mExitSepCut){
408  ofs<<mExitSepName.name<<","<<mExitSep[0]<<","<<mExitSep[1]<<"\t\t\t"<<" # pair ExitSep cut"<<endl;
409  printCutCounts(ofs,cutTypes[0],mExitSepCounter[0],mExitSepCounter[1]);
410  printCutCounts(ofs,cutTypes[1],mExitSepCounter[2],mExitSepCounter[3]);
411  }
412 
413  if(mQualityCut){
414  ofs<<mQualityName.name<<","<<mQuality[0]<<","<<mQuality[1]<<"\t\t\t"<<" # pair Quality cut"<<endl;
415  printCutCounts(ofs,cutTypes[0],mQualityCounter[0],mQualityCounter[1]);
416  printCutCounts(ofs,cutTypes[1],mQualityCounter[2],mQualityCounter[3]);
417  }
418 
419  if(mMidTpcSepLSCut){
420  ofs<<mMidTpcSepLSName.name<<","<<mMidTpcSepLS[0]<<","<<mMidTpcSepLS[1]<<"\t\t\t"<<" # pair MidTpcSepLS cut"<<endl;
421  printCutCounts(ofs,cutTypes[0],msplitLSCounter[0],msplitLSCounter[1]);
422  printCutCounts(ofs,cutTypes[1],msplitLSCounter[2],msplitLSCounter[3]);
423  }
424 
425  if(mMidTpcSepUSCut){
426  ofs<<mMidTpcSepUSName.name<<","<<mMidTpcSepUS[0]<<","<<mMidTpcSepUS[1]<<"\t\t\t"<<" # pair MidTpcSepUS cut"<<endl;
427  printCutCounts(ofs,cutTypes[0],msplitUSCounter[0],msplitUSCounter[1]);
428  printCutCounts(ofs,cutTypes[1],msplitUSCounter[2],msplitUSCounter[3]);
429  }
430 
431  if(mHBTCut){
432  ofs<<mHBTName.name<<","<<mHBT[0]<<","<<mHBT[1]<<","<<mHBT[2]<<","<<mHBT[3]<<"\t\t"<<" # pair HBT cut"<<endl;
433  printCutCounts(ofs,cutTypes[0],mHBTCounter[0],mHBTCounter[1]);
434  printCutCounts(ofs,cutTypes[1],mHBTCounter[2],mHBTCounter[3]);
435  }
436  if(mCoulombCut){
437  ofs<<mCoulombName.name<<","<<mCoulomb[0]<<","<<mCoulomb[1]<<","<<mCoulomb[2]<<","<<mCoulomb[3]<<"\t"<<" # pair Coulomb cut"<<endl;
438  printCutCounts(ofs,cutTypes[0],mCoulombCounter[0],mCoulombCounter[1]);
439  printCutCounts(ofs,cutTypes[1],mCoulombCounter[2],mCoulombCounter[3]);
440  }
441  if(mMergingCut){
442  ofs<<mMergingName.name<<","<<mMerging[0]<<","<<mMerging[1]<<"\t\t\t"<<" # pair Merging cut"<<endl;
443  printCutCounts(ofs,cutTypes[0],mMergingCounter[0],mMergingCounter[1]);
444  printCutCounts(ofs,cutTypes[1],mMergingCounter[2],mMergingCounter[3]);
445  }
446  if(mMergingCut2){
447  ofs<<mMergingName2.name<<","<<mMerging2[0]<<","<<mMerging2[1]<<"\t\t\t"<<" # pair Merging2 cut"<<endl;
448  printCutCounts(ofs,cutTypes[0],mMergingCounter2[0],mMergingCounter2[1]);
449  printCutCounts(ofs,cutTypes[1],mMergingCounter2[2],mMergingCounter2[3]);
450  }
451  if(mCrossingCut){
452  ofs<<mCrossingName.name<<","<<mCrossing[0]<<","<<mCrossing[1]<<"\t\t"<<" # pair Crossing cut"<<endl;
453  printCutCounts(ofs,cutTypes[0],mCrossingCounter[0],mCrossingCounter[1]);
454  printCutCounts(ofs,cutTypes[1],mCrossingCounter[2],mCrossingCounter[3]);
455  }
456  if(mCrossingCut2){
457  ofs<<mCrossingName2.name<<","<<mCrossing2[0]<<","<<mCrossing2[1]<<"\t\t"<<" # pair Crossing2 cut"<<endl;
458  printCutCounts(ofs,cutTypes[0],mCrossingCounter2[0],mCrossingCounter2[1]);
459  printCutCounts(ofs,cutTypes[1],mCrossingCounter2[2],mCrossingCounter2[3]);
460  }
461  if(mLUTCut){
462  ofs<<mLUTName.name<<","<<mLUTParams[0]<<","<<mLUTParams[1]<<"\t\t"<<" # pair LUT cut"<<endl;
463  printCutCounts(ofs,cutTypes[0],mLUTCounter[0],mLUTCounter[1]);
464  printCutCounts(ofs,cutTypes[1],mLUTCounter[2],mLUTCounter[3]);
465  }
466  if(mpionMomentumCut){
467  ofs<<mpionMomentumName.name<<","<<mdEdxMomentumCut[1][0]<<","<<mdEdxMomentumCut[1][1]<<"\t\t"<<" # pion momentum range cut"<<endl;
468  printCutCounts(ofs,cutTypes[0],mpionMomentumCounter[0],mpionMomentumCounter[1]);
469  printCutCounts(ofs,cutTypes[1],mpionMomentumCounter[2],mpionMomentumCounter[3]);
470  }
471  if(mKaonMomentumCut){
472  ofs<<mKaonMomentumName.name<<","<<mdEdxMomentumCut[2][0]<<","<<mdEdxMomentumCut[2][1]<<"\t\t"<<" # Kaon momentum range cut"<<endl;
473  printCutCounts(ofs,cutTypes[0],mKaonMomentumCounter[0],mKaonMomentumCounter[1]);
474  printCutCounts(ofs,cutTypes[1],mKaonMomentumCounter[2],mKaonMomentumCounter[3]);
475  }
476  if(mprotonMomentumCut){
477  ofs<<mprotonMomentumName.name<<","<<mdEdxMomentumCut[3][0]<<","<<mdEdxMomentumCut[3][1]<<"\t\t"<<" # proton momentum range cut"<<endl;
478  printCutCounts(ofs,cutTypes[0],mprotonMomentumCounter[0],mprotonMomentumCounter[1]);
479  printCutCounts(ofs,cutTypes[1],mprotonMomentumCounter[2],mprotonMomentumCounter[3]);
480  }
481 
482  if(mpionOtherMassCut){
483  ofs<<mpionOtherMassName.name<<","<<mpionOtherMass[0]<<","<<mpionOtherMass[1]<<"\t\t"<<" # pion-Other mass cut range cut"<<endl;
484  printCutCounts(ofs,cutTypes[0],mpionOtherMassCounter[0],mpionOtherMassCounter[1]);
485  printCutCounts(ofs,cutTypes[1],mpionOtherMassCounter[2],mpionOtherMassCounter[3]);
486  }
487  if(mpionpionMassCut){
488  ofs<<mpionpionMassName.name<<","<<mpionpionMass[0]<<","<<mpionpionMass[1]<<"\t\t"<<" # pion-pion mass cut range cut"<<endl;
489  printCutCounts(ofs,cutTypes[0],mpionpionMassCounter[0],mpionpionMassCounter[1]);
490  printCutCounts(ofs,cutTypes[1],mpionpionMassCounter[2],mpionpionMassCounter[3]);
491  }
492  if(mpionKaonMassCut){
493  ofs<<mpionKaonMassName.name<<","<<mpionKaonMass[0]<<","<<mpionKaonMass[1]<<"\t\t"<<" # pion-Kaon mass cut range cut"<<endl;
494  printCutCounts(ofs,cutTypes[0],mpionKaonMassCounter[0],mpionKaonMassCounter[1]);
495  printCutCounts(ofs,cutTypes[1],mpionKaonMassCounter[2],mpionKaonMassCounter[3]);
496  }
497  if(mpionprotonMassCut){
498  ofs<<mpionprotonMassName.name<<","<<mpionprotonMass[0]<<","<<mpionprotonMass[1]<<"\t\t"<<" # pion-proton mass cut range cut"<<endl;
499  printCutCounts(ofs,cutTypes[0],mpionprotonMassCounter[0],mpionprotonMassCounter[1]);
500  printCutCounts(ofs,cutTypes[1],mpionprotonMassCounter[2],mpionprotonMassCounter[3]);
501  }
502  if(mKaonOtherMassCut){
503  ofs<<mKaonOtherMassName.name<<","<<mKaonOtherMass[0]<<","<<mKaonOtherMass[1]<<"\t\t"<<" # pion-Other mass cut range cut"<<endl;
504  printCutCounts(ofs,cutTypes[0],mKaonOtherMassCounter[0],mKaonOtherMassCounter[1]);
505  printCutCounts(ofs,cutTypes[1],mKaonOtherMassCounter[2],mKaonOtherMassCounter[3]);
506  }
507  if(mKaonKaonMassCut){
508  ofs<<mKaonKaonMassName.name<<","<<mKaonKaonMass[0]<<","<<mKaonKaonMass[1]<<"\t\t"<<" # Kaon-Kaon mass cut range cut"<<endl;
509  printCutCounts(ofs,cutTypes[0],mKaonKaonMassCounter[0],mKaonKaonMassCounter[1]);
510  printCutCounts(ofs,cutTypes[1],mKaonKaonMassCounter[2],mKaonKaonMassCounter[3]);
511  }
512  if(mKaonprotonMassCut){
513  ofs<<mKaonprotonMassName.name<<","<<mKaonprotonMass[0]<<","<<mKaonprotonMass[1]<<"\t\t"<<" # Kaon-proton mass cut range cut"<<endl;
514  printCutCounts(ofs,cutTypes[0],mKaonprotonMassCounter[0],mKaonprotonMassCounter[1]);
515  printCutCounts(ofs,cutTypes[1],mKaonprotonMassCounter[2],mKaonprotonMassCounter[3]);
516  }
517  if(mprotonOtherMassCut){
518  ofs<<mprotonOtherMassName.name<<","<<mprotonOtherMass[0]<<","<<mprotonOtherMass[1]<<"\t\t"<<" # proton-Other mass cut range cut"<<endl;
519  printCutCounts(ofs,cutTypes[0],mprotonOtherMassCounter[0],mprotonOtherMassCounter[1]);
520  printCutCounts(ofs,cutTypes[1],mprotonOtherMassCounter[2],mprotonOtherMassCounter[3]);
521  }
522  if(mprotonprotonMassCut){
523  ofs<<mprotonprotonMassName.name<<","<<mprotonprotonMass[0]<<","<<mprotonprotonMass[1]<<"\t\t"<<" # proton-proton mass cut range cut"<<endl;
524  printCutCounts(ofs,cutTypes[0],mprotonprotonMassCounter[0],mprotonprotonMassCounter[1]);
525  printCutCounts(ofs,cutTypes[1],mprotonprotonMassCounter[2],mprotonprotonMassCounter[3]);
526  }
527  if(mOtherOtherMassCut){
528  ofs<<mOtherOtherMassName.name<<","<<mOtherOtherMass[0]<<","<<mOtherOtherMass[1]<<"\t\t"<<" # Other-Other mass cut range cut"<<endl;
529  printCutCounts(ofs,cutTypes[0],mOtherOtherMassCounter[0],mOtherOtherMassCounter[1]);
530  printCutCounts(ofs,cutTypes[1],mOtherOtherMassCounter[2],mOtherOtherMassCounter[3]);
531  }
532 
533 
534 
535  // ofs<<"# ******************************************** "<<endl<<endl;
536 
537 }
538 
539 //------------------------------------------------------------
540 int StEStructPairCuts::cutPair(){
541  // return 0 to use pair, return 1 to cut pair
542 
543  // if(cutDeltaPhi() || cutDeltaEta() || cutDeltaMt()) return 1;
544  //if(goodDeltaPhi() && goodDeltaEta() && goodDeltaMt()) return 0;
545 
546  // *** new test ***
547  /* cout.precision(3);
548  cout <<MidTpcXYSeparation()<<"\t"<<MidTpcZSeparation()<<"\t"<<MidTpcSeparation();
549  cout << "::\t"<<NominalTpcXYExitSeparation()<<"\t"<<NominalTpcZExitSeparation()<<"\t"<<NominalTpcExitSeparation()<<endl;
550  //cout << endl;
551  if (cutCoulomb()) cout << "COUL\t"<<DeltaEta()<<"\t"<<DeltaPhi()<<"\t"<<DeltaPt()<<"\t"<<mTrack1->Pt()<<"\t"<<mTrack2->Pt()<<endl;
552  if (cutHBT()) cout << "HBT\t"<<mType<<endl;
553  if (cutMerging()) cout << "MERG\t"<<NominalTpcAvgXYSeparation()<<"\t"<<NominalTpcAvgZSeparation()<<endl;
554  if (cutCrossing()) {
555  cout << "CROS\t";
556  if (mTrack1->Charge()>0) cout << "+ ";
557  else cout << "- ";
558  if (mTrack2->Charge()>0) cout << "+\t";
559  else cout << "-\t";
560  cout <<MidTpcXYSeparation()<<"\t"<<MidTpcZSeparation();
561  float dphi = mTrack1->Phi()-mTrack2->Phi(); // signed DeltaPhi
562  float dpt = mTrack1->Pt()- mTrack2->Pt(); // signed DeltaPt
563  cout << dphi << "\t" << dpt << endl;
564  }*/
565  // ***
566 
567  if (goodDeltaZ() || goodDeltaXY()) {
568  return 0;
569  }
570  if(cutMerging() || cutCrossing()) {
571  return 1;
572  }
573  if(cutMerging2() || cutCrossing2()) {
574  return 1;
575  }
576  if(cutLUT()) {
577  return 1;
578  }
579  if(cutCoulomb() || cutHBT()) {
580  return 1;
581  }
582 
583  if (!mdeltaEta) {
584  mdeltaEta = fabs(DeltaEta()); // may have been set above
585  }
586 
587  if(mdeltaEta<0.03){
588 
589  //--> qInv and EntSep are combined for speed & small delta eta
590  if(cutqInvORNominalEntranceSep()) return 1;
591 
592  if(mType==1 || mType==3) {
593  if(cutMidTpcSepUS()) return 1;
594  } else {
595  if(cutMidTpcSepLS()) return 1;
596  }
597  }
598 
599  if (cutQuality()) {
600  return 1;
601  }
602 
603  if(cutMass()) return 1;
604 
605  // if(cutExitSep() || cutQuality()) return 1;
606  return 0;
607 }
608 
609 //------------------------------------------------------------
610 int StEStructPairCuts::cutPairHistograms(){
611 
612  //>>>>> djp I haven't been updating this code with to include my attempts at pair cuts.
613  //>>>>> Should do this sometime.
614 
615  // much much slower cut code - should be run on a small sub-sample to view
616  // results of histograms. Tricky piece here is the connection between pair
617  // pair types and cut and event inter-cut relations.
618  // where I need to set the value used for histogramming to be
619  // at a specific overflow: xhi+2.0*(xhi-xlo) when this particular cut is
620  // not in play for this specific Pair.
621 
622  int cutCount=0;
623 
624  cutCount+=cutDeltaPhiH();
625  cutCount+=cutDeltaEtaH();
626  cutCount+=cutDeltaMtH();
627 
628  cutCount+=cutHBT();
629  cutCount+=cutCoulomb();
630  cutCount+=cutMerging();
631  cutCount+=cutCrossing();
632 
633  if(mdeltaEta<0.03){
634  cutCount+=cutqInvH();
635  cutCount+=cutEntranceSepH();
636  if(mType==1 || mType==3){
637  cutCount+=cutMidTpcSepUSH();
638  if(mMidTpcSepLSCut)mvalues[mMidTpcSepLSName.idx]=mMidTpcSepLS[1]+2.0*(mMidTpcSepLS[1]-mMidTpcSepLS[0]);
639  } else {
640  cutCount+=cutMidTpcSepLSH();
641  if(mMidTpcSepUSCut)mvalues[mMidTpcSepUSName.idx]=mMidTpcSepUS[1]+2.0*(mMidTpcSepUS[1]-mMidTpcSepUS[0]);
642  }
643  } else {
644  if(mqInvCut)mvalues[mqInvName.idx]=mqInv[1]+2.0*(mqInv[1]-mqInv[0]);
645  if(mEntSepCut)mvalues[mEntSepName.idx]=mEntSep[1]+2.0*(mEntSep[1]-mEntSep[0]);
646  if(mMidTpcSepUSCut)mvalues[mMidTpcSepUSName.idx]=mMidTpcSepUS[1]+2.0*(mMidTpcSepUS[1]-mMidTpcSepUS[0]);
647  if(mMidTpcSepLSCut)mvalues[mMidTpcSepLSName.idx]=mMidTpcSepLS[1]+2.0*(mMidTpcSepLS[1]-mMidTpcSepLS[0]);
648  }
649 
650  cutCount+=cutExitSepH();
651  cutCount+=cutQualityH();
652 
653  fillHistograms((cutCount==0));
654 
655  return cutCount;
656 }
657 
658 StLorentzVectorF
659 StEStructPairCuts::fourMomentumSum() const
660 {
661  StLorentzVectorF temp = mTrack1->FourMomentum()+mTrack2->FourMomentum();
662  return temp;
663 }
664 //-----------------------------------------------------------------
665 StLorentzVectorF
666 StEStructPairCuts::fourMomentumDiff() const
667 {
668  StLorentzVectorF temp = mTrack1->FourMomentum()-mTrack2->FourMomentum();
669  return temp;
670 }
671 
672 double
673 StEStructPairCuts::quality() const {
674 
675 
676  unsigned long padRow1To24Track1 = mTrack1->TopologyMapData(0) & mapMask0;
677  unsigned long padRow25To45Track1 = mTrack1->TopologyMapData(1) & mapMask1;
678  unsigned long padRow1To24Track2 = mTrack2->TopologyMapData(0) & mapMask0;
679  unsigned long padRow25To45Track2 = mTrack2->TopologyMapData(1) & mapMask1;
680 
681  // AND logic
682  unsigned long bothPads1To24 = padRow1To24Track1 & padRow1To24Track2;
683  unsigned long bothPads25To45 = padRow25To45Track1 & padRow25To45Track2;
684  // XOR logic
685  unsigned long onePad1To24 = padRow1To24Track1 ^ padRow1To24Track2;
686  unsigned long onePad25To45 = padRow25To45Track1 ^ padRow25To45Track2;
687  int ibits;
688  int Quality = 0;
689 
690  double normQual = 0.0;
691 
692 
693  int MaxQuality = mTrack1->TopologyMapTPCNHits() + mTrack2->TopologyMapTPCNHits();
694 
695  for (ibits=8;ibits<=31;ibits++) {
696  // bitI = 0;
697  // bitI |= 1UL<<(ibits);
698  if ( onePad1To24 & bitI[ibits] ) {
699  Quality++;
700  // continue;
701  } else if ( bothPads1To24 & bitI[ibits] ) {
702  Quality--;
703  }
704  }
705 
706  for (ibits=0;ibits<=20;ibits++) {
707  // bitI = 0;
708  //bitI = 1UL<<(ibits);
709  if ( onePad25To45 & bitI[ibits] ) {
710  Quality++;
711  //continue;
712  } else if ( bothPads25To45 & bitI[ibits] ) {
713  Quality--;
714  }
715  }
716 
717  normQual = (double)Quality/((double)MaxQuality );
718  return ( normQual );
719 
720 }
721 
722 //--------------------------------------------------------
723 double
724 StEStructPairCuts::OpeningAngle() const {
725  StThreeVectorF p1 = mTrack1->FourMomentum().vect();
726  StThreeVectorF p2 = mTrack2->FourMomentum().vect();
727  double dAngInv = 57.296*acos((p1.dot(p2))/(p1.mag()*p2.mag()));
728  return (dAngInv);
729 }
730 
731 //--------------------------------------------------------
732 // Zoffset correction is needed for mixed events.
733 double
734 StEStructPairCuts::NominalTpcExitSeparation() const {
735  StThreeVectorF off(0,0,mZoffset);
736  StThreeVectorF diff = mTrack1->NominalTpcExitPoint() - mTrack2->NominalTpcExitPoint() + off;
737  return (double)(diff.mag());
738 }
739 
740 double
741 StEStructPairCuts::NominalTpcXYExitSeparation() const {
742  StThreeVectorF diff = mTrack1->NominalTpcExitPoint() - mTrack2->NominalTpcExitPoint();
743  return (double)(diff.perp());
744 }
745 
746 double
747 StEStructPairCuts::NominalTpcZExitSeparation() const {
748  // There is a problem when both tracks exit the TPC through the same endcap, then both exit
749  // points have z = (+/-) 200, and the z separation is exactly zero.
750  // Distinguish casese
751  // -1 One track crosses endcap.
752  // -2 Both tracks cross same endcap.
753  // -3 Both tracks cross endcap but different endcaps.
754  // Not sure what I was thinking, comparing floating point numbers with ==.
755  // if ( fabs(mTrack1->NominalTpcExitPoint().z())==200 || fabs(mTrack2->NominalTpcExitPoint().z())==200 ) return -1;
756  if ( mTrack1->EndCapOuter()!=0 || mTrack2->EndCapOuter()!=0 ) {
757  if ( mTrack1->EndCapOuter()==mTrack2->EndCapOuter() ) {
758  return -2;
759  } else if ( mTrack1->EndCapOuter()==0 || mTrack2->EndCapOuter()==0 ) {
760  return -1;
761  } else {
762  return -3;
763  }
764  }
765  double diff = mTrack1->NominalTpcExitPoint().z() - mTrack2->NominalTpcExitPoint().z();
766  return (double)(fabs(diff+mZoffset));
767 
768 }
769 double
770 StEStructPairCuts::TpcRadialEndCapSeparation() const {
771  // Use Radial dseparation at endcap in some cases when one or both tracks cross
772  // endcap before getting to some radius.
773  float rdiff = mTrack1->EndCapRadius() - mTrack2->EndCapRadius();
774  return (double)(fabs(rdiff));
775 
776 }
777 
778 //--------------------------------------------------------
779 double
780 StEStructPairCuts::NominalTpcEntranceSeparation() const {
781  StThreeVectorF off(0,0,mZoffset);
782  StThreeVectorF diff = mTrack1->NominalTpcEntrancePoint() - mTrack2->NominalTpcEntrancePoint() + off;
783  return (double)(diff.mag());
784 }
785 
786 double
787 StEStructPairCuts::NominalTpcXYEntranceSeparation() const {
788  StThreeVectorF diff = mTrack1->NominalTpcEntrancePoint() - mTrack2->NominalTpcEntrancePoint();
789  return (double)(diff.perp());
790 }
791 double
792 StEStructPairCuts::NominalTpcZEntranceSeparation() const {
793  double diff = mTrack1->NominalTpcEntrancePoint().z() - mTrack2->NominalTpcEntrancePoint().z();
794  return (double)(fabs(diff + mZoffset));
795 }
796 
797 //--------------------------------------------------------
798 double
799 StEStructPairCuts::NominalTpcAvgXYSeparation() const {
800  double x1=NominalTpcXYEntranceSeparation();
801  double x2=MidTpcXYSeparation();
802  double x3=NominalTpcXYExitSeparation();
803  return (x1+x2+x3)/3.;
804 }
805 double
806 StEStructPairCuts::NominalTpcAvgZSeparation() const {
807  double x1=NominalTpcZEntranceSeparation();
808  double x2=MidTpcZSeparation();
809  double x3=NominalTpcZExitSeparation();
810  if (x3==-1) return (x1+x2)/2.; // if particle exited the endcap, exlude from average
811  else return (x1+x2+x3)/3.;
812 }
813 bool StEStructPairCuts::TracksCrossInPhi() const {
814  // Try a direct (but must be quite a bit slower than the cross product)
815  // calculation of phi differences. I think the cross product method had too many corner cases.
816  // Really want dPhi between -pi and pi. I was making sure they were between -2pi and 2pi.
817  double phiEnt1 = mTrack1->NominalTpcEntrancePoint().phi();
818  double phiEnt2 = mTrack2->NominalTpcEntrancePoint().phi();
819  double dphiEnt = phiEnt1 - phiEnt2;
820  while (dphiEnt > M_PI) {
821  dphiEnt -= 2*M_PI;
822  }
823  while (dphiEnt < -M_PI) {
824  dphiEnt += 2*M_PI;
825  }
826  double phiExit1 = mTrack1->NominalTpcExitPoint().phi();
827  double phiExit2 = mTrack2->NominalTpcExitPoint().phi();
828  double dphiExit = phiExit1 - phiExit2;
829  while (dphiExit > M_PI) {
830  dphiExit -= 2*M_PI;
831  }
832  while (dphiExit < -M_PI) {
833  dphiExit += 2*M_PI;
834  }
835  if ((dphiEnt*dphiExit < 0) && (fabs(dphiEnt) < M_PI/4)) {
836  return true;
837  }
838  return false;
839 }
840 
841 //--------------------------------------------------------
842 double
843 StEStructPairCuts::MidTpcSeparation() const {
844  StThreeVectorF off(0,0,mZoffset);
845  StThreeVectorF diff = mTrack1->MidTpcPoint() - mTrack2->MidTpcPoint() + off;
846  //cout << "Mid XY Sep" << diff.perp() << endl;
847  return (double)(diff.mag());
848 }
849 
850 double
851 StEStructPairCuts::MidTpcXYSeparation() const {
852  StThreeVectorF diff = mTrack1->MidTpcPoint() - mTrack2->MidTpcPoint();
853  //cout << "Mid XY Sep" << diff.perp() << endl;
854  return (double)(diff.perp());
855 }
856 
857 double
858 StEStructPairCuts::MidTpcZSeparation() const {
859  double diff = mTrack1->MidTpcPoint().z() - mTrack2->MidTpcPoint().z();
860  //cout << "Mid Z Sep" << fabs(diff.z()) << endl;
861  return (double)(fabs(diff+mZoffset));
862 }
863 //--------------------------------------------------------
864 double StEStructPairCuts::OuterMidTpcSeparation() const {
865  StThreeVectorF off(0,0,mZoffset);
866  StThreeVectorF diff = mTrack1->OuterMidTpcPoint() - mTrack2->OuterMidTpcPoint() + off;
867  return (double)(diff.mag());
868 }
869 
870 double StEStructPairCuts::OuterMidTpcXYSeparation() const {
871  StThreeVectorF diff = mTrack1->OuterMidTpcPoint() - mTrack2->OuterMidTpcPoint();
872  return (double)(diff.perp());
873 }
874 
875 double StEStructPairCuts::OuterMidTpcZSeparation() const {
876  // See comments for NominalTpcZExitSeparation()
877  if ( mTrack1->EndCapOuterMid()!=0 || mTrack2->EndCapOuterMid()!=0 ) {
878  if ( mTrack1->EndCapOuterMid()==mTrack2->EndCapOuterMid() ) {
879  return -2;
880  } else if ( mTrack1->EndCapOuterMid()==0 || mTrack2->EndCapOuterMid()==0 ) {
881  return -1;
882  } else {
883  return -3;
884  }
885  }
886  double diff = mTrack1->OuterMidTpcPoint().z() - mTrack2->OuterMidTpcPoint().z();
887  return (double)(fabs(diff+mZoffset));
888 }
889 
890 
891 /***********************************************************************
892  *
893  * $Log: StEStructPairCuts.cxx,v $
894  * Revision 1.15 2012/11/16 21:22:27 prindle
895  * 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
896  * Binning: Add quality cut.
897  * CutBin: modify mode9
898  * PairCuts: modify goodDeltaZ for case of one track leaving via endcap.
899  *
900  * Revision 1.14 2010/09/02 21:24:08 prindle
901  * 2ptCorrelations: Fill histograms for event mixing information
902  * Option for common mixing buffer
903  * Switch to selectively fill QInv histograms (which take a long time)
904  * CutBin: Moved PID code to Track class from Pair class. Needed to update this code.
905  * PairCuts: Moved PID code from here to Track class.
906  * Removed unnecessary creation of StThreeVector which seem to take a long time
907  * Add ToF momentum cuts, modify dEdx momentum cuts. (Now allow dEdx to be
908  * be identified up to 15GeV/c, ToF up to 10GeV/c.)
909  *
910  * Revision 1.13 2010/03/02 21:45:28 prindle
911  * Had a problem with pair cuts when one track exited via endplate
912  * Calculate maxDEta properly
913  * Warning if you try turning histograms for pair cuts on
914  *
915  * Revision 1.12 2009/11/09 21:32:41 prindle
916  * Fix warnings about casting char * to a const char * by redeclaring as const char *.
917  *
918  * Revision 1.11 2009/05/08 00:09:55 prindle
919  * In 2ptCorrelations we added switches to select blocks of histograms to fill.
920  * (See constructor in StEStruct2ptCorrelations.cxx)
921  * Use a brute force method for checking crossing cuts. I had too many corner
922  * cases with my clever check.
923  * In Binning, change Yt limit and add methods for accessing number of histogram bins
924  * to use (used in Support)
925  *
926  * Revision 1.10 2008/12/02 23:45:06 prindle
927  * Changed switchYt to switchXX (etc.) to better reflect function.
928  * Change minYt to 1.0 in Binning so YtYt histogram doesn't have empty lower bin (pt = 0.164 for yt = 1.0)
929  * In CutBin: remove initPtBin
930  * add mode 8
931  * add notSymmetrized (used in Support)
932  * Added LUT (Look Up Table) for pair cuts. Experimental for now.
933  * Modified cutMerging2 (to look at track separation at a few radii)
934  * and cutCrossing2 so it doesn't accidentally reject almost back to back tracks.
935  *
936  * Revision 1.9 2008/03/19 22:06:01 prindle
937  * Added doInvariantMass flag.
938  * Added some plots in pairDensityHistograms.
939  * SetZOffset used to only be done when doPairDensity was true.
940  * Moved creating/copying pairDensity histograms to same place as other histograms.
941  * Added cutBinHistMode
942  * mode3 neck was defined as yt1<2.2 && yt2<2.2 (and not soft)
943  * now is 1.8<yt1<2.2 && 1.8<yt2<2.2
944  * Added gooddzdxy, Merging2 and Crossing2 to pair cuts.
945  *
946  * Revision 1.8 2007/11/26 19:55:25 prindle
947  * In 2ptCorrelations: Support for keeping all z-bins of selected centralities
948  * Change way \hat{p_t} is calculated for parent distributions in pid case.
949  * Binning Added parent binning (for \hat{p_t}
950  * CutBin: Mode 5 extensively modified.
951  * Added invariant mass cuts (probably a bad idea in general.)
952  *
953  * Revision 1.7 2007/05/27 22:45:03 msd
954  * Added new cut bin modes 2 (soft/hard SS/AS), 6 (z-vertex binning), and 7 (modes 2*6).
955  * Fixed bug in merging cut.
956  * Added a few histograms to 2pt corr.
957  *
958  * Revision 1.6 2007/01/26 17:17:10 msd
959  * 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.
960  *
961  * Revision 1.5 2006/10/02 22:21:02 prindle
962  * Store only quadrant of eta_Delta - phi_Delta array/histogram.
963  * Store half of eta_Sigma - phi_Delta array/histogram.
964  * This required modifications in Binning.
965  * I had a bug in the pair loop (which left +- not fully symmetrized)
966  * and had to make changes in cut bins for mode 5 (and 3 I think)
967  * when I fixed this.
968  * Also change crossing cut to use only two parameters, the sign of
969  * the magnetic field being taken from the MuDst.
970  *
971  * Revision 1.4 2006/04/04 22:10:13 porter
972  * a handful of changes (specific to correlations)
973  * - added StEStructQAHists so that if NOT input frm Maker, each analysis has its own
974  * - used ability to get any max,min val from the cut class - or z-vertex binning
975  * - put z-vertex binning into 1 place
976  * - switched back 1st line of pair cut method to keep pair if good, not to reject if bad.
977  * - Pair cut object is now pointer in correlations
978  * - some diagnostic printouts available from macro
979  * - Duncan's delta-phi binning change
980  *
981  * Revision 1.3 2005/09/14 17:14:25 msd
982  * Large update, added new pair-cut system, added pair density plots for new analysis mode (4), added event mixing cuts (rewrote buffer for this)
983  *
984  * Revision 1.2 2004/06/25 03:11:49 porter
985  * New cut-binning implementation and modified pair-cuts for chunhui to review
986  *
987  * Revision 1.1 2003/10/15 18:20:46 porter
988  * initial check in of Estruct Analysis maker codes.
989  *
990  *
991  *********************************************************************/
992 
993 
994 
995 
Generated by   doxygen 1.8.5