StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
AnaCuts.cxx
1 #include <Riostream.h>
2 #include <TString.h>
3 #include <TMath.h>
4 #include <assert.h>
5 
6 #include <StEmcPool/StPhotonCommon/MyEvent.h>
7 #include <StEmcPool/StPhotonCommon/MyPoint.h>
8 //#include <StEmcPool/StPhotonCommon/MyMcTrack.h>
9 
10 #include "AnaCuts.h"
11 using namespace std;
12 ClassImp(AnaCuts)
13 
14 AnaCuts::AnaCuts(const char* flag)
15 {
16  cout<<"constructing CUTS for "<<flag<<endl;
17 
18  if(strcmp(flag,"dAu")==0){
19 
20  vertexxCUT=5.0;
21  vertexyCUT=5.0;
22  vertexzCUT=60.;
23  neutralCUT=5.;//10, note this is cm's, 1 strip = 1.7cm or so.
24  photonCUT=15.;
25  etaHitsCUT=2;//2;
26  phiHitsCUT=1;//2;
27  energyRatioCUT=0.;//0
28 
29  etaMinCUT=.1;//.1
30  etaMaxCUT=.9;//.9
31  rapidityMinCUT=.1;//0.
32  rapidityMaxCUT=.9;//1.
33  rapPionMinCUT=.1;//.1
34  rapPionMaxCUT=.9;//.9
35  asymmetryCUT=0.7;//.7
36  ratioCUT=0.8;//.8
37 
38  softTrigHT1=0.;
39  softTrigHT2=0.;
40 
41  nPtBinsMB=8;
42  ptBinsMB.Set(nPtBinsMB+1);
43  for(Int_t i=0;i<=nPtBinsMB;i++) ptBinsMB[i]=(Double_t)0.5*i;
44  ptBinsMB[7]=4.;
45  ptBinsMB[8]=5.;
46 
47  nPtBinsHT1=14;
48  ptBinsHT1.Set(nPtBinsHT1+1);
49  for(Int_t i=0;i<=nPtBinsHT1;i++) ptBinsHT1[i]=(Double_t)1.*i;
50  ptBinsHT1[13]=(Double_t)13.5;
51  ptBinsHT1[14]=(Double_t)15.;
52 
53  nPtBinsHT2=14;
54  ptBinsHT2.Set(nPtBinsHT2+1);
55  for(Int_t i=0;i<=nPtBinsHT2;i++) ptBinsHT2[i]=(Double_t)1.*i;
56  ptBinsHT2[13]=(Double_t)13.5;
57  ptBinsHT2[14]=(Double_t)15.;
58 
59 
60  nPtBinsEffMB=nPtBinsMB;//MB
61  ptBinsEffMB=ptBinsMB;//MB
62  nPtBinsEffHT1=nPtBinsHT1;
63  ptBinsEffHT1=ptBinsHT1;
64  nPtBinsEffHT2=nPtBinsHT2;
65  ptBinsEffHT2=ptBinsHT2;
66 
67  nMinvBinsMB=500;
68  mInvLowMB=0.0;
69  mInvHighMB=5.0;
70  nMinvBinsHT1=250;
71  mInvLowHT1=0.0;
72  mInvHighHT1=5.0;
73  nMinvBinsHT2=250;
74  mInvLowHT2=0.0;
75  mInvHighHT2=5.0;
76 
77  nPtBinsEff=40;
78  ptBinsEff.Set(nPtBinsEff+1);
79  for(Int_t i=0;i<=nPtBinsEff;i++){
80  ptBinsEff[i]=i*0.5;
81  }
82 
83  ht1AdcCUT=287;
84  ht2AdcCUT=447;
85 
86  isHot.Set(4800);
87  isHot.Reset();
88  Int_t badIdsDAU[]={82, 156, 240, 264, 719, 733, 818, 853, 931, 1289, 1355, 1358, 1359, 1472,
89  1527, 1759, 1943, 1965, 2005, 2011, 2025, 2090, 2130, 2131, 2152, 2170, 2369, -1};
90  Int_t entryDAU=0;
91  while(badIdsDAU[entryDAU]>-1) {isHot[badIdsDAU[entryDAU]-1]=1; entryDAU++;}
92 
93  timesSigma=2.0;
94  }
95  else if(strcmp(flag,"pp05")==0 || strcmp(flag,"pythia")==0){
96 
97  vertexxCUT=5.0;
98  vertexyCUT=5.0;
99  vertexzCUT=60.;//60.
100  neutralCUT=5.;//10, note this is cm's, 1 strip = 1.7cm or so.
101  photonCUT=15.;
102  etaHitsCUT=2;//2;
103  phiHitsCUT=1;//1;
104  energyRatioCUT=0.;//0
105 
106  etaMinCUT=.2;//.1
107  etaMaxCUT=.8;//.9
108  rapidityMinCUT=.1;//0.
109  rapidityMaxCUT=.9;//1.
110  rapPionMinCUT=.1;//.1
111  rapPionMaxCUT=.9;//.9
112  asymmetryCUT=.7;//.7
113  ratioCUT=1.;
114 
115  softTrigHT1=0.;
116  softTrigHT2=0.;
117 
118  nPtBinsMB=8;
119  ptBinsMB.Set(nPtBinsMB+1);
120  for(Int_t i=0;i<=nPtBinsMB;i++) ptBinsMB[i]=(Double_t)0.5*i;
121  ptBinsMB[7]=4.0;
122  ptBinsMB[8]=5.0;
123 
124 
125  nPtBinsHT1=14;
126  ptBinsHT1.Set(nPtBinsHT1+1);
127  for(Int_t i=0;i<=nPtBinsHT1;i++) ptBinsHT1[i]=(Double_t)1.*i;
128  ptBinsHT1[13]=(Double_t)13.5;
129  ptBinsHT1[14]=(Double_t)15.0;
130 
131  nPtBinsHT2=14;
132  ptBinsHT2.Set(nPtBinsHT2+1);
133  for(Int_t i=0;i<=nPtBinsHT2;i++) ptBinsHT2[i]=(Double_t)1.*i;
134  ptBinsHT2[13]=(Double_t)13.5;
135  ptBinsHT2[14]=(Double_t)15.0;
136 
137  nPtBinsEffMB=nPtBinsMB;//MB
138  ptBinsEffMB=ptBinsMB;//MB
139  nPtBinsEffHT1=nPtBinsHT1;
140  ptBinsEffHT1=ptBinsHT1;
141  nPtBinsEffHT2=nPtBinsHT2;
142  ptBinsEffHT2=ptBinsHT2;
143 
144  nMinvBinsMB=500;
145  mInvLowMB=0.0;
146  mInvHighMB=5.0;
147  nMinvBinsHT1=250;
148  mInvLowHT1=0.0;
149  mInvHighHT1=5.0;
150  nMinvBinsHT2=250;
151  mInvLowHT2=0.0;
152  mInvHighHT2=5.0;
153 
154  nPtBinsEff=40;
155  ptBinsEff.Set(nPtBinsEff+1);
156  for(Int_t i=0;i<=nPtBinsEff;i++){
157  ptBinsEff[i]=i*0.5;
158  }
159 
160  ht1AdcCUT=0;//420;//prelim. rejection of obvious bad triggers
161  ht2AdcCUT=0;//550;//but for pp05 trigger must be simulated.
162 
163  isHot.Set(4800);
164  isHot.Reset();
165  Int_t badIdsPP[]={750, 757, 762, 768, 812, -1};
166  Int_t entryPP=0;
167  while(badIdsPP[entryPP]>-1) {isHot[badIdsPP[entryPP]-1]=1; entryPP++;}
168 
169  timesSigma=2.0;
170  }
171  else{
172  cout<<"error constructing cuts"<<endl;
173  assert(0);
174  }
175 
176 
177  if(0){
178  nPtBinsEffMB=40;
179  ptBinsEffMB.Set(nPtBinsEffMB+1);
180  for(Int_t i=0;i<=nPtBinsEffMB;i++) ptBinsEffMB[i]=(Double_t).5*i;
181  nPtBinsEffHT1=40;
182  ptBinsEffHT1.Set(nPtBinsEffHT1+1);
183  for(Int_t i=0;i<=nPtBinsEffHT1;i++) ptBinsEffHT1[i]=(Double_t).5*i;
184  nPtBinsEffHT2=40;
185  ptBinsEffHT2.Set(nPtBinsEffHT2+1);
186  for(Int_t i=0;i<=nPtBinsEffHT2;i++) ptBinsEffHT2[i]=(Double_t).5*i;
187  }
188 
189 }
190 AnaCuts::~AnaCuts()
191 {
192  cout<<"cuts destructed!"<<endl;
193 }
194 Bool_t AnaCuts::isPointOK(MyPoint *p,TVector3 vpos)
195 {
196  TVector3 pos=p->position();
197  TVector3 mom=pos-vpos;
198  mom.SetMag(p->energy());
199  if(p->distanceToTrack()<neutralCUT) return kFALSE;
200  if(!(p->nHitsEta()>=etaHitsCUT && p->nHitsPhi()>=phiHitsCUT))
201  {
202  if(!(p->nHitsEta()>=phiHitsCUT && p->nHitsPhi()>=etaHitsCUT))
203  {
204  return kFALSE;
205  }
206  }
207  if(pos.PseudoRapidity()<etaMinCUT) return kFALSE;
208  if(pos.PseudoRapidity()>etaMaxCUT) return kFALSE;
209 
210  //new
211  if((p->energyEta()+p->energyPhi())/p->energy()<energyRatioCUT) return kFALSE;
212 
213  return kTRUE;
214 }
215 
216 Bool_t AnaCuts::isEventOK(MyEvent *ev, const char *flag)
217 {
218  TVector3 vPos=ev->vertex();
219  Float_t ratioE=TMath::Abs(ev->energyInBarrel()+ev->momentumInTpc())>0. ?
220  ev->energyInBarrel()/(ev->energyInBarrel()+ev->momentumInTpc()) : 0.;
221  //event cuts
222  if(strcmp(flag,"dAu")==0){
223  if(TMath::Abs(vPos.Z())>vertexzCUT) return kFALSE;
224  }
225  if(strcmp(flag,"pp05")==0){
226  if(TMath::Abs(ev->bbcVertexZ())>vertexzCUT) return kFALSE;
227  }
228  if(TMath::Abs(vPos.X())>vertexxCUT) return kFALSE;
229  if(TMath::Abs(vPos.Y())>vertexyCUT) return kFALSE;
230 
231  if(ev->trigger()<1) return kFALSE;
232  //hot towers and beam bg:
233  if(ev->trigger()&2 || ev->trigger()&4){
234  if(ev->highTowerId()>0 && isHot[ev->highTowerId() - 1]){
235  if(ev->trigger()&1) ev->setTrigger(1);
236  else return false;
237  }
238  if(ev->highTowerEnergy()<0.001){
239  if(ev->trigger()&1) ev->setTrigger(1);
240  else return false;
241  }
242  if(ratioE>ratioCUT){
243  if(ev->trigger()&1) ev->setTrigger(1);
244  else return false;
245  }
246  }
247  if(ev->trigger()&4 && ev->highTowerAdc()<=ht2AdcCUT){
248  ev->setTrigger(ev->trigger()-4);
249  }
250  if(ev->trigger()&2 && ev->highTowerAdc()<=ht1AdcCUT){
251  ev->setTrigger(1);
252  }
253 
254  return true;
255 }
256 void AnaCuts::printCuts()
257 {
258  cout<<"----- event cuts ------"<<endl;
259  cout<<"vertexxCUT="<<vertexxCUT<<endl;
260  cout<<"vertexyCUT="<<vertexyCUT<<endl;
261  cout<<"vertexzCUT="<<vertexzCUT<<endl;
262  cout<<"ratioCUT="<<ratioCUT<<endl<<endl;
263 
264  cout<<"----- point cuts ------"<<endl;
265  cout<<"neutralCUT="<<neutralCUT<<endl;
266  cout<<"chargedCUT="<<photonCUT<<endl;
267  cout<<"etaHitsCUT="<<etaHitsCUT<<endl;
268  cout<<"phiHitsCUT="<<phiHitsCUT<<endl;
269  cout<<"etaMinCUT="<<etaMinCUT<<endl;
270  cout<<"etaMaxCUT="<<etaMaxCUT<<endl;
271  cout<<"rapidityMinCUT="<<rapidityMinCUT<<endl;
272  cout<<"rapidityMaxCUT="<<rapidityMaxCUT<<endl;
273  cout<<"rapPionMinCUT="<<rapPionMinCUT<<endl;
274  cout<<"rapPionMaxCUT="<<rapPionMaxCUT<<endl;
275  cout<<"asymmetryCUT="<<asymmetryCUT<<endl;
276  cout<<"ht1AdcCUT="<<ht1AdcCUT<<endl;
277  cout<<"ht2AdcCUT="<<ht2AdcCUT<<endl<<endl;
278 
279 
280  cout<<endl<<"------ binning ------"<<endl;
281 
282  //bins
283  cout<<"mb bins={";
284  for(Int_t i=0;i<=nPtBinsMB;i++) cout<<ptBinsMB[i]<<",";
285  cout<<"};"<<endl;
286  cout<<"ht1 bins={";
287  for(Int_t i=0;i<=nPtBinsHT1;i++) cout<<ptBinsHT1[i]<<",";
288  cout<<"};"<<endl;
289  cout<<"ht2 bins={";
290  for(Int_t i=0;i<=nPtBinsHT2;i++) cout<<ptBinsHT2[i]<<",";
291  cout<<"};"<<endl;
292  cout<<"Eff bins={";
293  for(Int_t i=0;i<=nPtBinsEff;i++) cout<<ptBinsEff[i]<<",";
294  cout<<"};"<<endl;
295 
296  cout<<"inv mass bins mb: "<<nMinvBinsMB<<" "<<mInvLowMB<<" "<<mInvHighMB<<endl;
297  cout<<"inv mass bins ht1: "<<nMinvBinsHT1<<" "<<mInvLowHT1<<" "<<mInvHighHT1<<endl;
298  cout<<"inv mass bins ht2: "<<nMinvBinsHT2<<" "<<mInvLowHT2<<" "<<mInvHighHT2<<endl;
299 
300  //hot towers:
301  cout<<endl<<endl<<"------- hot towers -------"<<endl;
302  for(Int_t k=1;k<=4800;k++)
303  {
304  if(isHot[k-1]) cout<<k<<", ";
305  }
306  cout<<endl<<"-----------------------"<<endl<<endl<<endl;
307 }
MyPoint
Definition: MyPoint.h:7
Generated by   doxygen 1.8.5