00001
00002
00003
00004
00005
00006
00007
00008
00010
00011 #include "NCUtils/Cuts/NCAnalysisCutsNC.h"
00012 #include "NCUtils/NCType.h"
00013 #include "NCUtils/NCUtility.h"
00014 #include "AnalysisNtuples/ANtpEventInfoNC.h"
00015 #include "AnalysisNtuples/ANtpShowerInfoNC.h"
00016 #include "AnalysisNtuples/ANtpBeamInfo.h"
00017 #include "AnalysisNtuples/ANtpHeaderInfo.h"
00018 #include "AnalysisNtuples/ANtpTrackInfoNC.h"
00019 #include "AnalysisNtuples/ANtpTruthInfoBeam.h"
00020 #include "AnalysisNtuples/ANtpRecoInfo.h"
00021 #include "MessageService/MsgService.h"
00022 #include "Conventions/Detector.h"
00023 #include "Conventions/Munits.h"
00024 #include "Conventions/SimFlag.h"
00025 #include "Conventions/PlaneView.h"
00026 #include "Conventions/PlaneCoverage.h"
00027 #include "DataUtil/PlaneOutline.h"
00028
00029 #include "TMath.h"
00030
00031 using namespace BeamType;
00032 using std::vector;
00033 using NC::Utility::SQR;
00034
00035 CVSID("$Id: NCAnalysisCutsNC.cxx,v 1.14 2010/02/15 14:50:13 pittam Exp $");
00036
00037
00038 NCAnalysisCutsNC::NCAnalysisCutsNC() :
00039 NCAnalysisCuts()
00040 {
00041 }
00042
00043
00044 NCAnalysisCutsNC::~NCAnalysisCutsNC()
00045 {
00046 }
00047
00048
00049 bool NCAnalysisCutsNC::IsGoodBeamEvent()
00050 {
00051
00052
00053
00054
00055 ReleaseType::Release_t softType( NCAnalysisCuts::GetReleaseType() );
00056 Bool_t isMC( ReleaseType::IsMC(softType) );
00057
00058
00059
00060 if ( fEventInfo.header->detector == Detector::kNear ) return true;
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085 {
00086
00087 if (fEventInfo.header->events == 0) return false;
00088
00089 if (fEventInfo.header->snarlPulseHeight==0.) return false;
00090 Float_t evtPHfrac(fEventInfo.event->pulseHeight
00091 / fEventInfo.header->snarlPulseHeight);
00092 if ( ( fEventInfo.header->events == 2 && evtPHfrac < 0.75)
00093 || fEventInfo.header->events > 2 ) return false;
00094 }
00095
00096
00097 if(IsCosmicInSpillOx()) return false;
00098
00099 if(IsFibreNoiseInSpillOx()) return false;
00100
00101 if(IsLIInSpillOx()) return false;
00102
00103
00104 if(!isMC){
00105 double dt_spill = fEventInfo.event->stripTime1st * Munits::s;
00106 dt_spill -= fEventInfo.beam->nearestNSToSpill * Munits::ns;
00107
00108 if ( dt_spill < -2. * Munits::microsecond
00109 || dt_spill > 12. * Munits::microsecond ){
00110 return false;
00111 }
00112 }
00113
00114
00115 return true;
00116 }
00117
00118
00119 bool NCAnalysisCutsNC::IsStoppingBeamMuon()
00120 {
00121 bool stoppingMuon = false;
00122
00123 bool endZOK = false;
00124 bool xyOK = false;
00125
00126 if(fEventInfo.header->detector==(int)Detector::kFar){
00127 if((fEventInfo.track->endZ > 0.5
00128 && fEventInfo.track->endZ < 12.)
00129 || (fEventInfo.track->endZ > 16.5
00130 && fEventInfo.track->endZ < 28.))
00131 endZOK = true;
00132
00133 if(TMath::Sqrt(fEventInfo.track->endX*fEventInfo.track->endX
00134 +fEventInfo.track->endY*fEventInfo.track->endY) < 3.5)
00135 xyOK = true;
00136
00137 if(TMath::Abs(fEventInfo.track->traceEndZ) > 0.5 && xyOK && endZOK)
00138 stoppingMuon = true;
00139 }
00140 else if(fEventInfo.header->detector==(int)Detector::kNear){
00141
00142
00143
00144
00145 double u = TMath::Sqrt(0.5)*(fEventInfo.track->endX+fEventInfo.track->endY);
00146 double v = TMath::Sqrt(0.5)*(fEventInfo.track->endY-fEventInfo.track->endX);
00147
00148 if(0.3 < u && u < 1.8
00149 && -1.8 < v && v < -0.3
00150 && fEventInfo.track->endX < 2.4
00151 && TMath::Sqrt(u*u + v*v) > 0.8
00152 && fEventInfo.track->endZ < 15. ) stoppingMuon = true;
00153
00154 }
00155
00156 return stoppingMuon;
00157 }
00158
00159
00160 bool NCAnalysisCutsNC::InBeamFiducialVolume()
00161 {
00162
00163 if (fEventInfo.header->detector == (int)Detector::kFar){
00164
00165
00166
00167
00168
00169
00170 Float_t dist_edge(fEventInfo.event->vtxMetersToCloseEdge);
00171 Float_t dist_coil(fEventInfo.event->vtxMetersToCoil);
00172 Float_t vtxZ(fEventInfo.event->vtxZ);
00173 if ( (fEventInfo.event->tracks > 0)
00174 && ( (fEventInfo.track->planes - fEventInfo.shower->planes)>0)) {
00175 dist_edge = fEventInfo.track->vtxMetersToCloseEdge;
00176 dist_coil = fEventInfo.track->vtxMetersToCoil;
00177 vtxZ = fEventInfo.track->vtxZ - NCType::kTrackVtxAdjustment;
00178 }
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191 if ( vtxZ < 0.21 * Munits::m || vtxZ > 28.96 * Munits::m ) return false;
00192 if ( vtxZ > 13.72 * Munits::m && vtxZ < 16.12 * Munits::m ) return false;
00193
00194
00195
00196 if ( dist_edge < 0.4 * Munits::m ) return false;
00197
00198
00199
00200 if ( dist_coil < 0.6 * Munits::m ) return false;
00201 }
00202
00203
00204 else if(fEventInfo.header->detector == (int)Detector::kNear){
00205
00206
00207
00208
00209
00210 Float_t dist_edge(fEventInfo.event->vtxMetersToCloseEdge);
00211 Float_t vtxZ(fEventInfo.event->vtxZ);
00212 if ( (fEventInfo.event->tracks > 0)
00213 && ( (fEventInfo.track->planes - fEventInfo.shower->planes)>0) ) {
00214 dist_edge = fEventInfo.track->vtxMetersToCloseEdge;
00215 vtxZ = fEventInfo.track->vtxZ - NCType::kTrackVtxAdjustment;
00216 }
00217
00218
00219
00220
00221
00222 if ( vtxZ > 4.7368 * Munits::m || vtxZ < 1.7 * Munits::m ) return false;
00223
00224
00225
00226
00227 if ( dist_edge < 0.5 * Munits::m ) return false;
00228
00229
00230 }
00231
00232
00233 return true;
00234 }
00235
00236
00237 bool NCAnalysisCutsNC::InFiducialVolumeTrue()
00238 {
00239 Detector::Detector_t det = Detector::kFar;
00240 PlaneCoverage::PlaneCoverage_t outline = PlaneCoverage::kComplete;
00241 if(fEventInfo.header->detector == (int)Detector::kNear){
00242 det = Detector::kNear;
00243 outline = PlaneCoverage::kNearPartial;
00244 }
00245
00246 PlaneOutline po;
00247 float distU = 0.;
00248 float distV = 0.;
00249 float xedge = 0.;
00250 float yedge = 0.;
00251 float x = fEventInfo.truth->nuVtxX;
00252 float y = fEventInfo.truth->nuVtxY;
00253 float z = fEventInfo.truth->nuVtxZ;
00254
00255 po.DistanceToOuterEdge(x, y, PlaneView::kU,
00256 outline, distU, xedge, yedge);
00257 po.DistanceToOuterEdge(x, y, PlaneView::kV,
00258 outline, distV, xedge, yedge);
00259
00260 bool insideU = false;
00261
00262 if ( outline != PlaneCoverage::kNearPartial ) {
00263
00264
00265 insideU = ( (x*x + y*y) < 0.5 * Munits::m2 );
00266
00267
00268 }
00269
00270 bool insideV = insideU;
00271 if ( !insideU ){
00272 insideU = po.IsInside(x, y, PlaneView::kU, outline);
00273 insideV = po.IsInside(x, y, PlaneView::kV, outline);
00274 }
00275 distU *= ( insideU ? 1. : -1. );
00276 distV *= ( insideV ? 1. : -1. );
00277
00278 float vtxMetersToCloseEdge = 0.5*(distU + distV);
00279 float vtxMetersToCoil = TMath::Sqrt(x*x + y*y);
00280
00281
00282
00283 if(det == Detector::kFar){
00284
00285
00286
00287
00288
00289
00290
00291 if ( z < 0.2328*Munits::m || z > (28.8608+0.119) * Munits::m ) return false;
00292 if ( z > (13.6208+0.119)*Munits::m && z < 16.1408 * Munits::m ) return false;
00293
00294
00295
00296 if ( vtxMetersToCloseEdge < 0.4 * Munits::m ) return false;
00297
00298
00299
00300 if ( vtxMetersToCoil < 0.6 * Munits::m ) return false;
00301 }
00302
00303
00304 else if(det == Detector::kNear){
00305
00306
00307
00308
00309
00310 if ( z > 4.7368 * Munits::m || z < 1.728 * Munits::m ) return false;
00311
00312
00313
00314
00315 if ( vtxMetersToCloseEdge < 0.5 * Munits::m ) return false;
00316
00317
00318 }
00319
00320 return true;
00321 }
00322
00323
00324 bool NCAnalysisCutsNC::PassesFinalSelection()
00325 {
00326 if(!InBeamFiducialVolume()) return false;
00327 return true;
00328 }
00329
00330
00331 bool NCAnalysisCutsNC::PassesFinalSelection(ANtpRecoInfo *recoInfo)
00332 {
00333 if(recoInfo->inFiducialVolume < 1) return false;
00334 return true;
00335 }
00336
00337
00338 bool NCAnalysisCutsNC::IsFibreNoiseInSpillOx()
00339 {
00340 ReleaseType::Release_t softType(GetReleaseType());
00341
00342
00343 if(ReleaseType::IsCedar(softType))
00344 {
00345
00346 if(fEventInfo.event->pulseHeight < 2500 &&
00347 fEventInfo.event->totalStrips <= 8) return true;
00348
00349
00350
00351 if(!ReleaseType::IsBirch(softType))
00352 {
00353 if(fEventInfo.event->pulseHeight < 5000 &&
00354 fEventInfo.event->totalStrips <= 4) return true;
00355 }
00356 }
00357 else
00358 {
00359 if(fEventInfo.event->pulseHeight < 3750 &&
00360 fEventInfo.event->totalStrips <= 8) return true;
00361 if(fEventInfo.event->pulseHeight < 2000 &&
00362 fEventInfo.event->totalStrips > 8) return true;
00363
00364 }
00365 return false;
00366 }
00367
00368
00369
00370 bool NCAnalysisCutsNC::IsCosmicInSpillOx()
00371 {
00372
00373 Float_t strips_f(fEventInfo.event->totalStrips);
00374 Float_t planes_f(fEventInfo.event->planes);
00375 if ( planes_f == 0.0 || strips_f/(planes_f*planes_f) >= 1 ) return true;
00376
00377
00378
00379 if ( fEventInfo.event->showers > 0 ){
00380 Float_t tRMS2 = SQR(fEventInfo.shower->transverseRMSU);
00381 tRMS2 += SQR(fEventInfo.shower->transverseRMSV);
00382 Float_t tRMS = TMath::Sqrt( tRMS2 );
00383 if ( ( 0.3 + 0.1901*TMath::Log10(planes_f) ) < tRMS ) return true;
00384 }
00385
00386
00387 if ( fEventInfo.event->tracks > 0 ){
00388
00389 if ( fEventInfo.track->dcosZVtx < 0.4 ) return true;
00390
00391
00392 Bool_t upExitTrack(false);
00393 if ( ( fEventInfo.track->endZ > 28.78
00394 || fEventInfo.track->endMetersToCloseEdge < 0.5 )
00395 && fEventInfo.track->endY > -1.657 ) upExitTrack = true;
00396
00397 if ( upExitTrack
00398 && fEventInfo.track->dtdz < -0.5
00399 && -fEventInfo.track->dcosYVtx < -0.4 ) return true;
00400
00401 }
00402
00403
00404 return false;
00405 }
00406
00407
00408 bool NCAnalysisCutsNC::IsLIInSpillOx()
00409 {
00410
00411
00412 ReleaseType::Release_t softType(GetReleaseType());
00413
00414
00415 if(ReleaseType::IsCedar(softType))
00416 {
00417 return fEventInfo.header->isLIold || fEventInfo.header->triggerPMTTime > 0;
00418 }
00419 else
00420 {
00421 return fEventInfo.header->isLI || fEventInfo.header->triggerPMTTime > 0;
00422 }
00423 }