00001 #include "Mad/MadAbID.h"
00002
00003 #include "CandNtupleSR/NtpSRTrack.h"
00004 #include "CandNtupleSR/NtpSRShower.h"
00005 #include "CandNtupleSR/NtpSRStrip.h"
00006
00007 #include "MessageService/MsgService.h"
00008
00009 #include <cmath>
00010
00011 ClassImp(MadAbID)
00012
00013 CVSID("$Id: MadAbID.cxx,v 1.4 2007/06/19 21:38:04 petyt Exp $");
00014
00015 MadAbID::MadAbID()
00016 {
00017 fCosTheta_1D = 0;
00018 fX_1D = 0;
00019 fY_1D = 0;
00020 fQ2_1D = 0;
00021 fW2_1D = 0;
00022 fCosTheta_2D = 0;
00023 fX_2D = 0;
00024 fY_2D = 1;
00025 fQ2_2D = 0;
00026 fW2_2D = 0;
00027 fTrackPHfrac_1D = 0;
00028 fTrackPHmean_1D = 0;
00029 fTrackQPsigmaQP_1D = 0;
00030 fTrackLikePlanes_1D = 0;
00031 fTrackPHfrac_2D = 0;
00032 fTrackPHmean_2D = 1;
00033 fTrackQPsigmaQP_2D = 1;
00034 fTrackLikePlanes_2D = 1;
00035 fTrackCharge_1D = 1;
00036 fTrackEnergy_1D = 1;
00037 fEventEnergy_1D = 0;
00038 fNormalization = 1;
00039
00040 hCosTheta_1D_CC = 0;
00041 hX_1D_CC = 0;
00042 hY_1D_CC = 0;
00043 hQ2_1D_CC = 0;
00044 hW2_1D_CC = 0;
00045 hCosTheta_2D_CC = 0;
00046 hX_2D_CC = 0;
00047 hY_2D_CC = 0;
00048 hQ2_2D_CC = 0;
00049 hW2_2D_CC = 0;
00050 hE_2D_CC = 0;
00051 hTrackPHfrac_1D_CC = 0;
00052 hTrackPHmean_1D_CC = 0;
00053 hTrackQPsigmaQP_1D_CC = 0;
00054 hTrackLikePlanes_1D_CC = 0;
00055 hTrackPHfrac_2D_CC = 0;
00056 hTrackPHmean_2D_CC = 0;
00057 hTrackQPsigmaQP_2D_CC = 0;
00058 hTrackLikePlanes_2D_CC = 0;
00059 hTrackPlanes_2D_CC_1 = 0;
00060 hTrackPlanes_2D_CC_2 = 0;
00061 hTrackCharge_1D_CC = 0;
00062 hTrackEnergy_1D_CC = 0;
00063 hEventEnergy_1D_CC = 0;
00064
00065 hCosTheta_1D_NC = 0;
00066 hX_1D_NC = 0;
00067 hY_1D_NC = 0;
00068 hQ2_1D_NC = 0;
00069 hW2_1D_NC = 0;
00070 hCosTheta_2D_NC = 0;
00071 hX_2D_NC = 0;
00072 hY_2D_NC = 0;
00073 hQ2_2D_NC = 0;
00074 hW2_2D_NC = 0;
00075 hE_2D_NC = 0;
00076 hTrackPHfrac_1D_NC = 0;
00077 hTrackPHmean_1D_NC = 0;
00078 hTrackQPsigmaQP_1D_NC = 0;
00079 hTrackLikePlanes_1D_NC = 0;
00080 hTrackPHfrac_2D_NC = 0;
00081 hTrackPHmean_2D_NC = 0;
00082 hTrackQPsigmaQP_2D_NC = 0;
00083 hTrackLikePlanes_2D_NC = 0;
00084 hTrackPlanes_2D_NC_1 = 0;
00085 hTrackPlanes_2D_NC_2 = 0;
00086 hTrackCharge_1D_NC = 0;
00087 hTrackEnergy_1D_NC = 0;
00088 hEventEnergy_1D_NC = 0;
00089
00090 hNormalization = 0;
00091
00092 RecoRun = -1;
00093 RecoSnarl = -1;
00094 RecoEvent = -1;
00095
00096 PidRun = -1;
00097 PidSnarl = -1;
00098 PidEvent = -1;
00099
00100 fReadPdfs = 0;
00101 fWritePdfs = 0;
00102 fMakePdfs = 0;
00103 fFoundPdfs = 0;
00104
00105 fPdfsFileIn = "pdfs.in.root";
00106 fPdfsFileOut = "pdfs.out.root";
00107
00108 fStripList = new TObjArray[500];
00109 fTrackStripList = new TObjArray[500];
00110 fTrkShowerStripList = new TObjArray[500];
00111 fShwShowerStripList = new TObjArray[500];
00112 }
00113
00114 MadAbID::~MadAbID()
00115 {
00116 if(hCosTheta_1D_CC) delete hCosTheta_1D_CC;
00117 if(hX_1D_CC) delete hX_1D_CC;
00118 if(hY_1D_CC) delete hY_1D_CC;
00119 if(hQ2_1D_CC) delete hQ2_1D_CC;
00120 if(hW2_1D_CC) delete hW2_1D_CC;
00121 if(hCosTheta_2D_CC) delete hCosTheta_2D_CC;
00122 if(hX_2D_CC) delete hX_2D_CC;
00123 if(hY_2D_CC) delete hY_2D_CC;
00124 if(hQ2_2D_CC) delete hQ2_2D_CC;
00125 if(hW2_2D_CC) delete hW2_2D_CC;
00126 if(hE_2D_CC) delete hE_2D_CC;
00127 if(hTrackPHfrac_1D_CC) delete hTrackPHfrac_1D_CC;
00128 if(hTrackPHmean_1D_CC) delete hTrackPHmean_1D_CC;
00129 if(hTrackQPsigmaQP_1D_CC) delete hTrackQPsigmaQP_1D_CC;
00130 if(hTrackLikePlanes_1D_CC) delete hTrackLikePlanes_1D_CC;
00131 if(hTrackPHfrac_2D_CC) delete hTrackPHfrac_2D_CC;
00132 if(hTrackPHmean_2D_CC) delete hTrackPHmean_2D_CC;
00133 if(hTrackQPsigmaQP_2D_CC) delete hTrackQPsigmaQP_2D_CC;
00134 if(hTrackLikePlanes_2D_CC) delete hTrackLikePlanes_2D_CC;
00135 if(hTrackPlanes_2D_CC_1) delete hTrackPlanes_2D_CC_1;
00136 if(hTrackPlanes_2D_CC_2) delete hTrackPlanes_2D_CC_2;
00137 if(hTrackCharge_1D_CC) delete hTrackCharge_1D_CC;
00138 if(hTrackEnergy_1D_CC) delete hTrackEnergy_1D_CC;
00139 if(hEventEnergy_1D_CC) delete hEventEnergy_1D_CC;
00140
00141 if(hCosTheta_1D_NC) delete hCosTheta_1D_NC;
00142 if(hX_1D_NC) delete hX_1D_NC;
00143 if(hY_1D_NC) delete hY_1D_NC;
00144 if(hQ2_1D_NC) delete hQ2_1D_NC;
00145 if(hW2_1D_NC) delete hW2_1D_NC;
00146 if(hCosTheta_2D_NC) delete hCosTheta_2D_NC;
00147 if(hX_2D_NC) delete hX_2D_NC;
00148 if(hY_2D_NC) delete hY_2D_NC;
00149 if(hQ2_2D_NC) delete hQ2_2D_NC;
00150 if(hW2_2D_NC) delete hW2_2D_NC;
00151 if(hE_2D_NC) delete hE_2D_NC;
00152 if(hTrackPHfrac_1D_NC) delete hTrackPHfrac_1D_NC;
00153 if(hTrackPHmean_1D_NC) delete hTrackPHmean_1D_NC;
00154 if(hTrackQPsigmaQP_1D_NC) delete hTrackQPsigmaQP_1D_NC;
00155 if(hTrackLikePlanes_1D_NC) delete hTrackLikePlanes_1D_NC;
00156 if(hTrackPHfrac_2D_NC) delete hTrackPHfrac_2D_NC;
00157 if(hTrackPHmean_2D_NC) delete hTrackPHmean_2D_NC;
00158 if(hTrackQPsigmaQP_2D_NC) delete hTrackQPsigmaQP_2D_NC;
00159 if(hTrackLikePlanes_2D_NC) delete hTrackLikePlanes_2D_NC;
00160 if(hTrackPlanes_2D_NC_1) delete hTrackPlanes_2D_NC_1;
00161 if(hTrackPlanes_2D_NC_2) delete hTrackPlanes_2D_NC_2;
00162 if(hTrackCharge_1D_NC) delete hTrackCharge_1D_NC;
00163 if(hTrackEnergy_1D_NC) delete hTrackEnergy_1D_NC;
00164 if(hEventEnergy_1D_NC) delete hEventEnergy_1D_NC;
00165
00166 if(hNormalization) delete hNormalization;
00167
00168 delete [] fStripList;
00169 delete [] fTrackStripList;
00170 delete [] fTrkShowerStripList;
00171 delete [] fShwShowerStripList;
00172 }
00173
00174 Int_t MadAbID::PidCosTheta()
00175 {
00176 if( fCosTheta_1D ) return 1;
00177 else if( fCosTheta_2D ) return 2;
00178 else return 0;
00179 }
00180
00181 Int_t MadAbID::PidX()
00182 {
00183 if( fX_1D ) return 1;
00184 else if( fX_2D ) return 2;
00185 else return 0;
00186 }
00187
00188 Int_t MadAbID::PidY()
00189 {
00190 if( fY_1D ) return 1;
00191 else if( fY_2D ) return 2;
00192 else return 0;
00193 }
00194
00195 Int_t MadAbID::PidW2()
00196 {
00197 if( fW2_1D ) return 1;
00198 else if( fW2_2D ) return 2;
00199 else return 0;
00200 }
00201
00202 Int_t MadAbID::PidQ2()
00203 {
00204 if( fQ2_1D ) return 1;
00205 else if( fQ2_2D ) return 2;
00206 else return 0;
00207 }
00208
00209 Int_t MadAbID::PidTrackPHfrac()
00210 {
00211 if( fTrackPHfrac_1D ) return 1;
00212 else if( fTrackPHfrac_2D ) return 2;
00213 else return 0;
00214 }
00215
00216 Int_t MadAbID::PidTrackPHmean()
00217 {
00218 if( fTrackPHmean_1D ) return 1;
00219 else if( fTrackPHmean_2D ) return 2;
00220 else return 0;
00221 }
00222
00223 Int_t MadAbID::PidTrackQPsigmaQP()
00224 {
00225 if( fTrackQPsigmaQP_1D ) return 1;
00226 else if( fTrackQPsigmaQP_2D ) return 2;
00227 else return 0;
00228 }
00229
00230 Int_t MadAbID::PidTrackLikePlanes()
00231 {
00232 if( fTrackLikePlanes_1D ) return 1;
00233 else if( fTrackLikePlanes_2D ) return 2;
00234 else return 0;
00235 }
00236
00237 Int_t MadAbID::PidTrackCharge()
00238 {
00239 if( fTrackCharge_1D ) return 1;
00240 else return 0;
00241 }
00242
00243 Int_t MadAbID::PidTrackEnergy()
00244 {
00245 if( fTrackEnergy_1D ) return 1;
00246 else return 0;
00247 }
00248
00249 Int_t MadAbID::PidEventEnergy()
00250 {
00251 if( fEventEnergy_1D ) return 1;
00252 else return 0;
00253 }
00254
00255 Int_t MadAbID::PidNormalization()
00256 {
00257 if( fNormalization ) return 1;
00258 else return 0;
00259 }
00260
00261 void MadAbID::UseCosTheta(Int_t n)
00262 {
00263 switch( n ){
00264 case 0: fCosTheta_1D = 0; fCosTheta_2D = 0; break;
00265 case 1: fCosTheta_1D = 1; fCosTheta_2D = 0; break;
00266 case 2: fCosTheta_1D = 0; fCosTheta_2D = 1; break;
00267 default: fCosTheta_1D = 0; fCosTheta_2D = 1; break;
00268 }
00269 }
00270
00271 void MadAbID::UseX(Int_t n)
00272 {
00273 switch( n ){
00274 case 0: fX_1D = 0; fX_2D = 0; break;
00275 case 1: fX_1D = 1; fX_2D = 0; break;
00276 case 2: fX_1D = 0; fX_2D = 1; break;
00277 default: fX_1D = 0; fX_2D = 1; break;
00278 }
00279 }
00280
00281 void MadAbID::UseY(Int_t n)
00282 {
00283 switch( n ){
00284 case 0: fY_1D = 0; fY_2D = 0; break;
00285 case 1: fY_1D = 1; fY_2D = 0; break;
00286 case 2: fY_1D = 0; fY_2D = 1; break;
00287 default: fY_1D = 0; fY_2D = 1; break;
00288 }
00289 }
00290
00291 void MadAbID::UseW2(Int_t n)
00292 {
00293 switch( n ){
00294 case 0: fQ2_1D = 0; fQ2_2D = 0; break;
00295 case 1: fQ2_1D = 1; fQ2_2D = 0; break;
00296 case 2: fQ2_1D = 0; fQ2_2D = 1; break;
00297 default: fQ2_1D = 0; fQ2_2D = 1; break;
00298 }
00299 }
00300
00301 void MadAbID::UseQ2(Int_t n)
00302 {
00303 switch( n ){
00304 case 0: fW2_1D = 0; fW2_2D = 0; break;
00305 case 1: fW2_1D = 1; fW2_2D = 0; break;
00306 case 2: fW2_1D = 0; fW2_2D = 1; break;
00307 default: fW2_1D = 0; fW2_2D = 1; break;
00308 }
00309 }
00310
00311 void MadAbID::UseTrackPHfrac(Int_t n)
00312 {
00313 switch( n ){
00314 case 0: fTrackPHfrac_1D = 0; fTrackPHfrac_2D = 0; break;
00315 case 1: fTrackPHfrac_1D = 1; fTrackPHfrac_2D = 0; break;
00316 case 2: fTrackPHfrac_1D = 0; fTrackPHfrac_2D = 1; break;
00317 default: fTrackPHfrac_1D = 0; fTrackPHfrac_2D = 1; break;
00318 }
00319 }
00320
00321 void MadAbID::UseTrackPHmean(Int_t n)
00322 {
00323 switch( n ){
00324 case 0: fTrackPHmean_1D = 0; fTrackPHmean_2D = 0; break;
00325 case 1: fTrackPHmean_1D = 1; fTrackPHmean_2D = 0; break;
00326 case 2: fTrackPHmean_1D = 0; fTrackPHmean_2D = 1; break;
00327 default: fTrackPHmean_1D = 0; fTrackPHmean_2D = 1; break;
00328 }
00329 }
00330
00331 void MadAbID::UseTrackQPsigmaQP(Int_t n)
00332 {
00333 switch( n ){
00334 case 0: fTrackQPsigmaQP_1D = 0; fTrackQPsigmaQP_2D = 0; break;
00335 case 1: fTrackQPsigmaQP_1D = 1; fTrackQPsigmaQP_2D = 0; break;
00336 case 2: fTrackQPsigmaQP_1D = 0; fTrackQPsigmaQP_2D = 1; break;
00337 default: fTrackQPsigmaQP_1D = 0; fTrackQPsigmaQP_2D = 1; break;
00338 }
00339 }
00340
00341 void MadAbID::UseTrackLikePlanes(Int_t n)
00342 {
00343 switch( n ){
00344 case 0: fTrackLikePlanes_1D = 0; fTrackLikePlanes_2D = 0; break;
00345 case 1: fTrackLikePlanes_1D = 1; fTrackLikePlanes_2D = 0; break;
00346 case 2: fTrackLikePlanes_1D = 0; fTrackLikePlanes_2D = 1; break;
00347 default: fTrackLikePlanes_1D = 0; fTrackLikePlanes_2D = 1; break;
00348 }
00349 }
00350
00351 void MadAbID::UseTrackCharge(Int_t n)
00352 {
00353 switch( n ){
00354 case 0: fTrackCharge_1D = 0; break;
00355 case 1: fTrackCharge_1D = 1; break;
00356 default: fTrackCharge_1D = 1; break;
00357 }
00358 }
00359
00360 void MadAbID::UseTrackEnergy(Int_t n)
00361 {
00362 switch( n ){
00363 case 0: fTrackEnergy_1D = 0; break;
00364 case 1: fTrackEnergy_1D = 1; break;
00365 default: fTrackEnergy_1D = 1; break;
00366 }
00367 }
00368
00369 void MadAbID::UseEventEnergy(Int_t n)
00370 {
00371 switch( n ){
00372 case 0: fEventEnergy_1D = 0; break;
00373 case 1: fEventEnergy_1D = 1; break;
00374 default: fEventEnergy_1D = 1; break;
00375 }
00376 }
00377
00378 void MadAbID::UseNormalization(Int_t n)
00379 {
00380 switch( n ){
00381 case 0: fNormalization = 0; break;
00382 case 1: fNormalization = 1; break;
00383 default: fNormalization = 1; break;
00384 }
00385 }
00386
00387 void MadAbID::Reset()
00388 {
00389 fCosTheta_1D = 0;
00390 fX_1D = 0;
00391 fY_1D = 0;
00392 fQ2_1D = 0;
00393 fW2_1D = 0;
00394 fCosTheta_2D = 0;
00395 fX_2D = 0;
00396 fY_2D = 0;
00397 fQ2_2D = 0;
00398 fW2_2D = 0;
00399 fTrackPHfrac_1D = 0;
00400 fTrackPHmean_1D = 0;
00401 fTrackQPsigmaQP_1D = 0;
00402 fTrackLikePlanes_1D = 0;
00403 fTrackPHfrac_2D = 0;
00404 fTrackPHmean_2D = 0;
00405 fTrackQPsigmaQP_2D = 0;
00406 fTrackLikePlanes_2D = 0;
00407 fTrackCharge_1D = 0;
00408 fTrackEnergy_1D = 0;
00409 fEventEnergy_1D = 0;
00410 fNormalization = 0;
00411 }
00412
00413 Bool_t MadAbID::PassCuts(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00414 {
00415 MSG("MadAb",Msg::kDebug) << " *** MadAbID::PassCuts(...) *** " << endl;
00416
00417
00418
00419 this->MakeRecoVariables(ntpEvent,ntpRecord);
00420
00421
00422
00423
00424
00425
00426
00427 Bool_t passfail = 0;
00428
00429 MSG("MadAb",Msg::kVerbose) << " Pass/Fail Conditions: " << endl
00430 << " Track Planes = " << TRKplanes << endl
00431 << " Track Fit Pass/Fail = " << TRKFITpass << endl
00432 << " Track Vertex Containment = " << RECOvtxcontained << endl;
00433
00434 if( TRKplanes>0
00435
00436 && RECOvtxcontained ){
00437 passfail = 1;
00438 }
00439
00440 MSG("MadAb",Msg::kVerbose) << " Pass/Fail = " << passfail << endl;
00441
00442 return passfail;
00443 }
00444
00445 Double_t MadAbID::GetRecoCosTheta(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00446 {
00447
00448
00449 this->MakeRecoVariables(ntpEvent,ntpRecord);
00450
00451 return RECOdircosneu;
00452 }
00453
00454 Double_t MadAbID::GetRecoX(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00455 {
00456
00457
00458 this->MakeRecoVariables(ntpEvent,ntpRecord);
00459
00460 return RECOx;
00461 }
00462
00463 Double_t MadAbID::GetRecoY(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00464 {
00465
00466
00467 this->MakeRecoVariables(ntpEvent,ntpRecord);
00468
00469 return RECOy;
00470 }
00471
00472 Double_t MadAbID::GetRecoQ2(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00473 {
00474
00475
00476 this->MakeRecoVariables(ntpEvent,ntpRecord);
00477
00478 return RECOq2;
00479 }
00480
00481 Double_t MadAbID::GetRecoW2(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00482 {
00483
00484
00485 this->MakeRecoVariables(ntpEvent,ntpRecord);
00486
00487 return RECOw2;
00488 }
00489
00490 Double_t MadAbID::GetRecoE(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00491 {
00492
00493
00494 this->MakeRecoVariables(ntpEvent,ntpRecord);
00495
00496 return RECOenu;
00497 }
00498
00499 Int_t MadAbID::GetTrackPlanes(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00500 {
00501
00502
00503 this->MakeRecoVariables(ntpEvent,ntpRecord);
00504
00505 return TRKspanplanes;
00506 }
00507
00508 Int_t MadAbID::GetTrackLikePlanes(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00509 {
00510
00511
00512 this->MakeRecoVariables(ntpEvent,ntpRecord);
00513
00514 return TRKtrackplanes;
00515 }
00516
00517 Int_t MadAbID::GetTrackEMCharge(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00518 {
00519
00520
00521 this->MakeRecoVariables(ntpEvent,ntpRecord);
00522
00523 return TRKFITemcharge;
00524 }
00525
00526 Double_t MadAbID::GetTrackQPsigmaQP(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00527 {
00528
00529
00530 this->MakeRecoVariables(ntpEvent,ntpRecord);
00531
00532 return fabs(TRKFITqpsigmaqp);
00533 }
00534
00535 Double_t MadAbID::GetTrackPHfrac(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00536 {
00537
00538
00539 this->MakeRecoVariables(ntpEvent,ntpRecord);
00540
00541 Double_t reco_phfrac = -999.9;
00542 if( RECOtotph>0.0 ) reco_phfrac = TRKtotph/RECOtotph;
00543
00544 return reco_phfrac;
00545 }
00546
00547 Double_t MadAbID::GetTrackPHmean(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00548 {
00549
00550
00551 this->MakeRecoVariables(ntpEvent,ntpRecord);
00552
00553 Double_t reco_phmean = -999.9;
00554 if( TRKtrackplanes>2 ) reco_phmean = TRKtrkph/TRKtrackplanes;
00555 else if(TRKplanes>0 ) reco_phmean = TRKtotph/TRKplanes;
00556
00557 return reco_phmean;
00558 }
00559
00560 Double_t MadAbID::CalcPID(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00561 {
00562 return this->GetPid(ntpEvent,ntpRecord);
00563 }
00564
00565 Bool_t MadAbID::GetPass(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00566 {
00567
00568
00569 this->MakePidVariables(ntpEvent,ntpRecord);
00570
00571 return PIDPASSFAIL;
00572 }
00573
00574 Double_t MadAbID::GetPid(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00575 {
00576
00577
00578 this->MakePidVariables(ntpEvent,ntpRecord);
00579
00580 return PID;
00581 }
00582
00583 Double_t MadAbID::GetProbCC(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00584 {
00585
00586
00587 this->MakePidVariables(ntpEvent,ntpRecord);
00588
00589 return PROBCC;
00590 }
00591
00592 Double_t MadAbID::GetProbNC(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
00593 {
00594
00595
00596 this->MakePidVariables(ntpEvent,ntpRecord);
00597
00598 return PROBNC;
00599 }
00600
00601 void MadAbID::ReadPDFs(const char* filename)
00602 {
00603 MSG("MadAb",Msg::kDebug) << " *** MadAbID::ReadPDFs(...) *** " << endl;
00604
00605 fPdfsFileIn = filename;
00606 fReadPdfs = 1;
00607
00608 if( fReadPdfs ){
00609
00610
00611
00612 if(hCosTheta_1D_CC) delete hCosTheta_1D_CC;
00613 if(hX_1D_CC) delete hX_1D_CC;
00614 if(hY_1D_CC) delete hY_1D_CC;
00615 if(hQ2_1D_CC) delete hQ2_1D_CC;
00616 if(hW2_1D_CC) delete hW2_1D_CC;
00617 if(hCosTheta_2D_CC) delete hCosTheta_2D_CC;
00618 if(hX_2D_CC) delete hX_2D_CC;
00619 if(hY_2D_CC) delete hY_2D_CC;
00620 if(hQ2_2D_CC) delete hQ2_2D_CC;
00621 if(hW2_2D_CC) delete hW2_2D_CC;
00622 if(hE_2D_CC) delete hE_2D_CC;
00623 if(hTrackPHfrac_1D_CC) delete hTrackPHfrac_1D_CC;
00624 if(hTrackPHmean_1D_CC) delete hTrackPHmean_1D_CC;
00625 if(hTrackQPsigmaQP_1D_CC) delete hTrackQPsigmaQP_1D_CC;
00626 if(hTrackLikePlanes_1D_CC) delete hTrackLikePlanes_1D_CC;
00627 if(hTrackPHfrac_2D_CC) delete hTrackPHfrac_2D_CC;
00628 if(hTrackPHmean_2D_CC) delete hTrackPHmean_2D_CC;
00629 if(hTrackQPsigmaQP_2D_CC) delete hTrackQPsigmaQP_2D_CC;
00630 if(hTrackLikePlanes_2D_CC) delete hTrackLikePlanes_2D_CC;
00631 if(hTrackPlanes_2D_CC_1) delete hTrackPlanes_2D_CC_1;
00632 if(hTrackPlanes_2D_CC_2) delete hTrackPlanes_2D_CC_2;
00633 if(hTrackCharge_1D_CC) delete hTrackCharge_1D_CC;
00634 if(hTrackEnergy_1D_CC) delete hTrackEnergy_1D_CC;
00635 if(hEventEnergy_1D_CC) delete hEventEnergy_1D_CC;
00636
00637 if(hCosTheta_1D_NC) delete hCosTheta_1D_NC;
00638 if(hX_1D_NC) delete hX_1D_NC;
00639 if(hY_1D_NC) delete hY_1D_NC;
00640 if(hQ2_1D_NC) delete hQ2_1D_NC;
00641 if(hW2_1D_NC) delete hW2_1D_NC;
00642 if(hCosTheta_2D_NC) delete hCosTheta_2D_NC;
00643 if(hX_2D_NC) delete hX_2D_NC;
00644 if(hY_2D_NC) delete hY_2D_NC;
00645 if(hQ2_2D_NC) delete hQ2_2D_NC;
00646 if(hW2_2D_NC) delete hW2_2D_NC;
00647 if(hE_2D_NC) delete hE_2D_NC;
00648 if(hTrackPHfrac_1D_NC) delete hTrackPHfrac_1D_NC;
00649 if(hTrackPHmean_1D_NC) delete hTrackPHmean_1D_NC;
00650 if(hTrackQPsigmaQP_1D_NC) delete hTrackQPsigmaQP_1D_NC;
00651 if(hTrackLikePlanes_1D_NC) delete hTrackLikePlanes_1D_NC;
00652 if(hTrackPHfrac_2D_NC) delete hTrackPHfrac_2D_NC;
00653 if(hTrackPHmean_2D_NC) delete hTrackPHmean_2D_NC;
00654 if(hTrackQPsigmaQP_2D_NC) delete hTrackQPsigmaQP_2D_NC;
00655 if(hTrackLikePlanes_2D_NC) delete hTrackLikePlanes_2D_NC;
00656 if(hTrackPlanes_2D_NC_1) delete hTrackPlanes_2D_NC_1;
00657 if(hTrackPlanes_2D_NC_2) delete hTrackPlanes_2D_NC_2;
00658 if(hTrackCharge_1D_NC) delete hTrackCharge_1D_NC;
00659 if(hTrackEnergy_1D_NC) delete hTrackEnergy_1D_NC;
00660 if(hEventEnergy_1D_NC) delete hEventEnergy_1D_NC;
00661
00662 if(hNormalization) delete hNormalization;
00663
00664 fFoundPdfs = 0;
00665
00666
00667
00668 MSG("MadAb",Msg::kInfo) << " reading PDFs from file " << endl;
00669
00670 TDirectory* tmpd = 0;
00671 tmpd = gDirectory;
00672 TFile* file = new TFile(fPdfsFileIn.Data(),"read");
00673
00674 if( file->IsOpen() ){
00675 hCosTheta_1D_CC = (TH1D*)(file->Get("hCosTheta_1D_CC"));
00676 hX_1D_CC = (TH1D*)(file->Get("hX_1D_CC"));
00677 hY_1D_CC = (TH1D*)(file->Get("hY_1D_CC"));
00678 hQ2_1D_CC = (TH1D*)(file->Get("hQ2_1D_CC"));
00679 hW2_1D_CC = (TH1D*)(file->Get("hW2_1D_CC"));
00680 hCosTheta_2D_CC = (TH2D*)(file->Get("hCosTheta_2D_CC"));
00681 hX_2D_CC = (TH2D*)(file->Get("hX_2D_CC"));
00682 hY_2D_CC = (TH2D*)(file->Get("hY_2D_CC"));
00683 hQ2_2D_CC = (TH2D*)(file->Get("hQ2_2D_CC"));
00684 hW2_2D_CC = (TH2D*)(file->Get("hW2_2D_CC"));
00685 hE_2D_CC = (TH1D*)(file->Get("hE_2D_CC"));
00686 hTrackPHfrac_1D_CC = (TH1D*)(file->Get("hTrackPHfrac_1D_CC"));
00687 hTrackPHmean_1D_CC = (TH1D*)(file->Get("hTrackPHmean_1D_CC"));
00688 hTrackQPsigmaQP_1D_CC = (TH1D*)(file->Get("hTrackQPsigmaQP_1D_CC"));
00689 hTrackLikePlanes_1D_CC = (TH1D*)(file->Get("hTrackLikePlanes_1D_CC"));
00690 hTrackPHfrac_2D_CC = (TH2D*)(file->Get("hTrackPHfrac_2D_CC"));
00691 hTrackPHmean_2D_CC = (TH2D*)(file->Get("hTrackPHmean_2D_CC"));
00692 hTrackQPsigmaQP_2D_CC = (TH2D*)(file->Get("hTrackQPsigmaQP_2D_CC"));
00693 hTrackLikePlanes_2D_CC = (TH2D*)(file->Get("hTrackLikePlanes_2D_CC"));
00694 hTrackPlanes_2D_CC_1 = (TH1D*)(file->Get("hTrackPlanes_2D_CC_1"));
00695 hTrackPlanes_2D_CC_2 = (TH1D*)(file->Get("hTrackPlanes_2D_CC_2"));
00696 hTrackCharge_1D_CC = (TH1D*)(file->Get("hTrackCharge_1D_CC"));
00697 hTrackEnergy_1D_CC = (TH1D*)(file->Get("hTrackEnergy_1D_CC"));
00698 hEventEnergy_1D_CC = (TH1D*)(file->Get("hEventEnergy_1D_CC"));
00699 hCosTheta_1D_NC = (TH1D*)(file->Get("hCosTheta_1D_NC"));
00700 hX_1D_NC = (TH1D*)(file->Get("hX_1D_NC"));
00701 hY_1D_NC = (TH1D*)(file->Get("hY_1D_NC"));
00702 hQ2_1D_NC = (TH1D*)(file->Get("hQ2_1D_NC"));
00703 hW2_1D_NC = (TH1D*)(file->Get("hW2_1D_NC"));
00704 hCosTheta_2D_NC = (TH2D*)(file->Get("hCosTheta_2D_NC"));
00705 hX_2D_NC = (TH2D*)(file->Get("hX_2D_NC"));
00706 hY_2D_NC = (TH2D*)(file->Get("hY_2D_NC"));
00707 hQ2_2D_NC = (TH2D*)(file->Get("hQ2_2D_NC"));
00708 hW2_2D_NC = (TH2D*)(file->Get("hW2_2D_NC"));
00709 hE_2D_NC = (TH1D*)(file->Get("hE_2D_NC"));
00710 hTrackPHfrac_1D_NC = (TH1D*)(file->Get("hTrackPHfrac_1D_NC"));
00711 hTrackPHmean_1D_NC = (TH1D*)(file->Get("hTrackPHmean_1D_NC"));
00712 hTrackQPsigmaQP_1D_NC = (TH1D*)(file->Get("hTrackQPsigmaQP_1D_NC"));
00713 hTrackLikePlanes_1D_NC = (TH1D*)(file->Get("hTrackLikePlanes_1D_NC"));
00714 hTrackPHfrac_2D_NC = (TH2D*)(file->Get("hTrackPHfrac_2D_NC"));
00715 hTrackPHmean_2D_NC = (TH2D*)(file->Get("hTrackPHmean_2D_NC"));
00716 hTrackQPsigmaQP_2D_NC = (TH2D*)(file->Get("hTrackQPsigmaQP_2D_NC"));
00717 hTrackLikePlanes_2D_NC = (TH2D*)(file->Get("hTrackLikePlanes_2D_NC"));
00718 hTrackPlanes_2D_NC_1 = (TH1D*)(file->Get("hTrackPlanes_2D_NC_1"));
00719 hTrackPlanes_2D_NC_2 = (TH1D*)(file->Get("hTrackPlanes_2D_NC_2"));
00720 hTrackCharge_1D_NC = (TH1D*)(file->Get("hTrackCharge_1D_NC"));
00721 hTrackEnergy_1D_NC = (TH1D*)(file->Get("hTrackEnergy_1D_NC"));
00722 hEventEnergy_1D_NC = (TH1D*)(file->Get("hEventEnergy_1D_NC"));
00723 hNormalization = (TH1D*)(file->Get("hNormalization"));
00724 }
00725 else{
00726 MSG("MadAb",Msg::kWarning) << " ... unable to open file: " << fPdfsFileIn.Data() << endl;
00727 }
00728
00729 gDirectory = tmpd;
00730
00731 if( hCosTheta_1D_CC
00732 && hX_1D_CC
00733 && hY_1D_CC
00734 && hQ2_1D_CC
00735 && hW2_1D_CC
00736 && hCosTheta_2D_CC
00737 && hX_2D_CC
00738 && hY_2D_CC
00739 && hQ2_2D_CC
00740 && hW2_2D_CC
00741 && hE_2D_CC
00742 && hTrackPHfrac_1D_CC
00743 && hTrackPHmean_1D_CC
00744 && hTrackQPsigmaQP_1D_CC
00745 && hTrackLikePlanes_1D_CC
00746 && hTrackPHfrac_2D_CC
00747 && hTrackPHmean_2D_CC
00748 && hTrackQPsigmaQP_2D_CC
00749 && hTrackLikePlanes_2D_CC
00750 && hTrackPlanes_2D_CC_1
00751 && hTrackPlanes_2D_CC_2
00752 && hTrackCharge_1D_CC
00753 && hTrackEnergy_1D_CC
00754 && hEventEnergy_1D_CC
00755 && hCosTheta_1D_NC
00756 && hX_1D_NC
00757 && hY_1D_NC
00758 && hQ2_1D_NC
00759 && hW2_1D_NC
00760 && hCosTheta_2D_NC
00761 && hX_2D_NC
00762 && hY_2D_NC
00763 && hQ2_2D_NC
00764 && hW2_2D_NC
00765 && hE_2D_NC
00766 && hTrackPHfrac_1D_NC
00767 && hTrackPHmean_1D_NC
00768 && hTrackQPsigmaQP_1D_NC
00769 && hTrackLikePlanes_1D_NC
00770 && hTrackPHfrac_2D_NC
00771 && hTrackPHmean_2D_NC
00772 && hTrackQPsigmaQP_2D_NC
00773 && hTrackLikePlanes_2D_NC
00774 && hTrackPlanes_2D_NC_1
00775 && hTrackPlanes_2D_NC_2
00776 && hTrackCharge_1D_NC
00777 && hTrackEnergy_1D_NC
00778 && hEventEnergy_1D_NC
00779 && hNormalization ){
00780 MSG("MadAb",Msg::kInfo) << " ... found PDFs from file " << endl;
00781 fFoundPdfs = 1;
00782 }
00783 }
00784 }
00785
00786 void MadAbID::WritePDFs(const char* filename)
00787 {
00788 MSG("MadAb",Msg::kDebug) << " *** MadAbID::WritePDFs(...) *** " << endl;
00789
00790 fPdfsFileOut = filename;
00791 fWritePdfs = 1;
00792
00793
00794
00795 if( fWritePdfs ){
00796
00797
00798
00799 TFile* file = new TFile(fPdfsFileOut.Data(),"recreate");
00800
00801
00802
00803
00804
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828
00829
00830
00831
00832
00833
00834
00835
00836
00837
00838
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858 file->cd();
00859 hCosTheta_1D_CC->Write();
00860 hX_1D_CC->Write();
00861 hY_1D_CC->Write();
00862 hQ2_1D_CC->Write();
00863 hW2_1D_CC->Write();
00864 hCosTheta_2D_CC->Write();
00865 hX_2D_CC->Write();
00866 hY_2D_CC->Write();
00867 hQ2_2D_CC->Write();
00868 hW2_2D_CC->Write();
00869 hE_2D_CC->Write();
00870 hTrackPHfrac_1D_CC->Write();
00871 hTrackPHmean_1D_CC->Write();
00872 hTrackQPsigmaQP_1D_CC->Write();
00873 hTrackLikePlanes_1D_CC->Write();
00874 hTrackPHfrac_2D_CC->Write();
00875 hTrackPHmean_2D_CC->Write();
00876 hTrackQPsigmaQP_2D_CC->Write();
00877 hTrackLikePlanes_2D_CC->Write();
00878 hTrackPlanes_2D_CC_1->Write();
00879 hTrackPlanes_2D_CC_2->Write();
00880 hTrackCharge_1D_CC->Write();
00881 hTrackEnergy_1D_CC->Write();
00882 hEventEnergy_1D_CC->Write();
00883
00884 hCosTheta_1D_NC->Write();
00885 hX_1D_NC->Write();
00886 hY_1D_NC->Write();
00887 hQ2_1D_NC->Write();
00888 hW2_1D_NC->Write();
00889 hCosTheta_2D_NC->Write();
00890 hX_2D_NC->Write();
00891 hY_2D_NC->Write();
00892 hQ2_2D_NC->Write();
00893 hW2_2D_NC->Write();
00894 hE_2D_NC->Write();
00895 hTrackPHfrac_1D_NC->Write();
00896 hTrackPHmean_1D_NC->Write();
00897 hTrackQPsigmaQP_1D_NC->Write();
00898 hTrackLikePlanes_1D_NC->Write();
00899 hTrackPHfrac_2D_NC->Write();
00900 hTrackPHmean_2D_NC->Write();
00901 hTrackQPsigmaQP_2D_NC->Write();
00902 hTrackLikePlanes_2D_NC->Write();
00903 hTrackPlanes_2D_NC_1->Write();
00904 hTrackPlanes_2D_NC_2->Write();
00905 hTrackCharge_1D_NC->Write();
00906 hTrackEnergy_1D_NC->Write();
00907 hEventEnergy_1D_NC->Write();
00908
00909 hNormalization->Write();
00910
00911
00912
00913 file->Close();
00914
00915
00916
00917 }
00918 }
00919
00920 void MadAbID::NormalizePDFs()
00921 {
00922 MSG("MadAb",Msg::kDebug) << " *** MadAbID::NormalizePDFs() *** " << endl;
00923
00924 if( fFoundPdfs ){
00925 MSG("MadAb",Msg::kDebug) << " OVERALL NORMALIZATION: CC=" << hNormalization->GetBinContent(hNormalization->FindBin(1)) << " NC=" << hNormalization->GetBinContent(hNormalization->FindBin(0)) << endl;
00926
00927 MSG("MadAb",Msg::kDebug) << endl
00928 << " normalizing histograms [before]: " << endl
00929 << " hCosTheta_1D_CC: " << hCosTheta_1D_CC->Integral("width") << endl
00930 << " hX_1D_CC: " << hX_1D_CC->Integral("width") << endl
00931 << " hY_1D_CC: " << hY_1D_CC->Integral("width") << endl
00932 << " hQ2_1D_CC: " << hQ2_1D_CC->Integral("width") << endl
00933 << " hW2_1D_CC: " << hW2_1D_CC->Integral("width") << endl
00934 << " hCosTheta_2D_CC: " << hCosTheta_2D_CC->Integral("width") << endl
00935 << " hX_2D_CC: " << hX_2D_CC->Integral("width") << endl
00936 << " hY_2D_CC: " << hY_2D_CC->Integral("width") << endl
00937 << " hQ2_2D_CC: " << hQ2_2D_CC->Integral("width") << endl
00938 << " hW2_2D_CC: " << hW2_2D_CC->Integral("width") << endl
00939 << " hE_2D_CC: " << hE_2D_CC->Integral("width") << endl
00940 << " hTrackPHfrac_1D_CC: " << hTrackPHfrac_1D_CC->Integral("width") << endl
00941 << " hTrackPHmean_1D_CC: " << hTrackPHmean_1D_CC->Integral("width") << endl
00942 << " hTrackQPsigmaQP_1D_CC: " << hTrackQPsigmaQP_1D_CC->Integral("width") << endl
00943 << " hTrackLikePlanes1D_CC: " << hTrackLikePlanes_1D_CC->Integral("width") << endl
00944 << " hTrackPHfrac_2D_CC: " << hTrackPHfrac_2D_CC->Integral("width") << endl
00945 << " hTrackPHmean_2D_CC: " << hTrackPHmean_2D_CC->Integral("width") << endl
00946 << " hTrackQPsigmaQP_2D_CC: " << hTrackQPsigmaQP_2D_CC->Integral("width") << endl
00947 << " hTrackLikePlanes_2D_CC: " << hTrackLikePlanes_2D_CC->Integral("width") << endl
00948 << " hTrackPlanes_2D_CC_1: " << hTrackPlanes_2D_CC_1->Integral("width") << endl
00949 << " hTrackPlanes_2D_CC_2: " << hTrackPlanes_2D_CC_2->Integral("width") << endl
00950 << " hTrackCharge_1D_CC: " << hTrackCharge_1D_CC->Integral("width") << endl
00951 << " hTrackEnergy_1D_CC: " << hTrackEnergy_1D_CC->Integral("width") << endl
00952 << " hEventEnergy_1D_CC: " << hEventEnergy_1D_CC->Integral("width") << endl
00953 << " hCosTheta_1D_NC: " << hCosTheta_1D_NC->Integral("width") << endl
00954 << " hX_1D_NC: " << hX_1D_NC->Integral("width") << endl
00955 << " hY_1D_NC: " << hY_1D_NC->Integral("width") << endl
00956 << " hQ2_1D_NC: " << hQ2_1D_NC->Integral("width") << endl
00957 << " hW2_1D_NC: " << hW2_1D_NC->Integral("width") << endl
00958 << " hCosTheta_2D_NC: " << hCosTheta_2D_NC->Integral("width") << endl
00959 << " hX_2D_NC: " << hX_2D_NC->Integral("width") << endl
00960 << " hY_2D_NC: " << hY_2D_NC->Integral("width") << endl
00961 << " hQ2_2D_NC: " << hQ2_2D_NC->Integral("width") << endl
00962 << " hW2_2D_NC: " << hW2_2D_NC->Integral("width") << endl
00963 << " hE_2D_NC: " << hE_2D_NC->Integral("width") << endl
00964 << " hTrackPHfrac_1D_NC: " << hTrackPHfrac_1D_NC->Integral("width") << endl
00965 << " hTrackPHmean_1D_NC: " << hTrackPHmean_1D_NC->Integral("width") << endl
00966 << " hTrackQPsigmaQP_1D_NC: " << hTrackQPsigmaQP_1D_NC->Integral("width") << endl
00967 << " hTrackLikePlanes_1D_NC: " << hTrackLikePlanes_1D_NC->Integral("width") << endl
00968 << " hTrackPHfrac_2D_NC: " << hTrackPHfrac_2D_NC->Integral("width") << endl
00969 << " hTrackPHmean_2D_NC: " << hTrackPHmean_2D_NC->Integral("width") << endl
00970 << " hTrackQPsigmaQP_2D_NC: " << hTrackQPsigmaQP_2D_NC->Integral("width") << endl
00971 << " hTrackLikePlanes_2D_NC: " << hTrackLikePlanes_2D_NC->Integral("width") << endl
00972 << " hTrackPlanes_2D_NC_1: " << hTrackPlanes_2D_NC_1->Integral("width") << endl
00973 << " hTrackPlanes_2D_NC_2: " << hTrackPlanes_2D_NC_2->Integral("width") << endl
00974 << " hTrackCharge_1D_NC: " << hTrackCharge_1D_NC->Integral("width") << endl
00975 << " hTrackEnergy_1D_NC: " << hTrackEnergy_1D_NC->Integral("width") << endl
00976 << " hEventEnergy_1D_NC: " << hEventEnergy_1D_NC->Integral("width") << endl
00977 << " hNormalization: " << hNormalization->Integral("width") << endl;
00978
00979 if(hCosTheta_1D_CC->Integral("width")>0.0) hCosTheta_1D_CC->Scale(1.0/hCosTheta_1D_CC->Integral("width"));
00980 if(hX_1D_CC->Integral("width")>0.0) hX_1D_CC->Scale(1.0/hX_1D_CC->Integral("width"));
00981 if(hY_1D_CC->Integral("width")>0.0) hY_1D_CC->Scale(1.0/hY_1D_CC->Integral("width"));
00982 if(hQ2_1D_CC->Integral("width")>0.0) hQ2_1D_CC->Scale(1.0/hQ2_1D_CC->Integral("width"));
00983 if(hW2_1D_CC->Integral("width")>0.0) hW2_1D_CC->Scale(1.0/hW2_1D_CC->Integral("width"));
00984 if(hCosTheta_2D_CC->Integral("width")>0.0) hCosTheta_2D_CC->Scale(1.0/hCosTheta_2D_CC->Integral("width"));
00985 if(hX_2D_CC->Integral("width")>0.0) hX_2D_CC->Scale(1.0/hX_2D_CC->Integral("width"));
00986 if(hY_2D_CC->Integral("width")>0.0) hY_2D_CC->Scale(1.0/hY_2D_CC->Integral("width"));
00987 if(hQ2_2D_CC->Integral("width")>0.0) hQ2_2D_CC->Scale(1.0/hQ2_2D_CC->Integral("width"));
00988 if(hW2_2D_CC->Integral("width")>0.0) hW2_2D_CC->Scale(1.0/hW2_2D_CC->Integral("width"));
00989 if(hE_2D_CC->Integral("width")>0.0) hE_2D_CC->Scale(1.0/hE_2D_CC->Integral("width"));
00990 if(hTrackPHfrac_1D_CC->Integral("width")>0.0) hTrackPHfrac_1D_CC->Scale(1.0/hTrackPHfrac_1D_CC->Integral("width"));
00991 if(hTrackPHmean_1D_CC->Integral("width")>0.0) hTrackPHmean_1D_CC->Scale(1.0/hTrackPHmean_1D_CC->Integral("width"));
00992 if(hTrackQPsigmaQP_1D_CC->Integral("width")>0.0) hTrackQPsigmaQP_1D_CC->Scale(1.0/hTrackQPsigmaQP_1D_CC->Integral("width"));
00993 if(hTrackLikePlanes_1D_CC->Integral("width")>0.0) hTrackLikePlanes_1D_CC->Scale(1.0/hTrackLikePlanes_1D_CC->Integral("width"));
00994 if(hTrackPHfrac_2D_CC->Integral("width")>0.0) hTrackPHfrac_2D_CC->Scale(1.0/hTrackPHfrac_2D_CC->Integral("width"));
00995 if(hTrackPHmean_2D_CC->Integral("width")>0.0) hTrackPHmean_2D_CC->Scale(1.0/hTrackPHmean_2D_CC->Integral("width"));
00996 if(hTrackQPsigmaQP_2D_CC->Integral("width")>0.0) hTrackQPsigmaQP_2D_CC->Scale(1.0/hTrackQPsigmaQP_2D_CC->Integral("width"));
00997 if(hTrackLikePlanes_2D_CC->Integral("width")>0.0) hTrackLikePlanes_2D_CC->Scale(1.0/hTrackLikePlanes_2D_CC->Integral("width"));
00998 if(hTrackPlanes_2D_CC_1->Integral("width")>0.0) hTrackPlanes_2D_CC_1->Scale(1.0/hTrackPlanes_2D_CC_1->Integral("width"));
00999 if(hTrackPlanes_2D_CC_2->Integral("width")>0.0) hTrackPlanes_2D_CC_2->Scale(1.0/hTrackPlanes_2D_CC_2->Integral("width"));
01000 if(hTrackCharge_1D_CC->Integral("width")>0.0) hTrackCharge_1D_CC->Scale(1.0/hTrackCharge_1D_CC->Integral("width"));
01001 if(hTrackEnergy_1D_CC->Integral("width")>0.0) hTrackEnergy_1D_CC->Scale(1.0/hTrackEnergy_1D_CC->Integral("width"));
01002 if(hEventEnergy_1D_CC->Integral("width")>0.0) hEventEnergy_1D_CC->Scale(1.0/hEventEnergy_1D_CC->Integral("width"));
01003 if(hCosTheta_1D_NC->Integral("width")>0.0) hCosTheta_1D_NC->Scale(1.0/hCosTheta_1D_NC->Integral("width"));
01004 if(hX_1D_NC->Integral("width")>0.0) hX_1D_NC->Scale(1.0/hX_1D_NC->Integral("width"));
01005 if(hY_1D_NC->Integral("width")>0.0) hY_1D_NC->Scale(1.0/hY_1D_NC->Integral("width"));
01006 if(hQ2_1D_NC->Integral("width")>0.0) hQ2_1D_NC->Scale(1.0/hQ2_1D_NC->Integral("width"));
01007 if(hW2_1D_NC->Integral("width")>0.0) hW2_1D_NC->Scale(1.0/hW2_1D_NC->Integral("width"));
01008 if(hCosTheta_2D_NC->Integral("width")>0.0) hCosTheta_2D_NC->Scale(1.0/hCosTheta_2D_NC->Integral("width"));
01009 if(hX_2D_NC->Integral("width")>0.0) hX_2D_NC->Scale(1.0/hX_2D_NC->Integral("width"));
01010 if(hY_2D_NC->Integral("width")>0.0) hY_2D_NC->Scale(1.0/hY_2D_NC->Integral("width"));
01011 if(hQ2_2D_NC->Integral("width")>0.0) hQ2_2D_NC->Scale(1.0/hQ2_2D_NC->Integral("width"));
01012 if(hW2_2D_NC->Integral("width")>0.0) hW2_2D_NC->Scale(1.0/hW2_2D_NC->Integral("width"));
01013 if(hE_2D_NC->Integral("width")>0.0) hE_2D_NC->Scale(1.0/hE_2D_NC->Integral("width"));
01014 if(hTrackPHfrac_1D_NC->Integral("width")>0.0) hTrackPHfrac_1D_NC->Scale(1.0/hTrackPHfrac_1D_NC->Integral("width"));
01015 if(hTrackPHmean_1D_NC->Integral("width")>0.0) hTrackPHmean_1D_NC->Scale(1.0/hTrackPHmean_1D_NC->Integral("width"));
01016 if(hTrackQPsigmaQP_1D_NC->Integral("width")>0.0) hTrackQPsigmaQP_1D_NC->Scale(1.0/hTrackQPsigmaQP_1D_NC->Integral("width"));
01017 if(hTrackLikePlanes_1D_NC->Integral("width")>0.0) hTrackLikePlanes_1D_NC->Scale(1.0/hTrackLikePlanes_1D_NC->Integral("width"));
01018 if(hTrackPHfrac_2D_NC->Integral("width")>0.0) hTrackPHfrac_2D_NC->Scale(1.0/hTrackPHfrac_2D_NC->Integral("width"));
01019 if(hTrackPHmean_2D_NC->Integral("width")>0.0) hTrackPHmean_2D_NC->Scale(1.0/hTrackPHmean_2D_NC->Integral("width"));
01020 if(hTrackQPsigmaQP_2D_NC->Integral("width")>0.0) hTrackQPsigmaQP_2D_NC->Scale(1.0/hTrackQPsigmaQP_2D_NC->Integral("width"));
01021 if(hTrackLikePlanes_2D_NC->Integral("width")>0.0) hTrackLikePlanes_2D_NC->Scale(1.0/hTrackLikePlanes_2D_NC->Integral("width"));
01022 if(hTrackPlanes_2D_NC_1->Integral("width")>0.0) hTrackPlanes_2D_NC_1->Scale(1.0/hTrackPlanes_2D_NC_1->Integral("width"));
01023 if(hTrackPlanes_2D_NC_2->Integral("width")>0.0) hTrackPlanes_2D_NC_2->Scale(1.0/hTrackPlanes_2D_NC_2->Integral("width"));
01024 if(hTrackCharge_1D_NC->Integral("width")>0.0) hTrackCharge_1D_NC->Scale(1.0/hTrackCharge_1D_NC->Integral("width"));
01025 if(hTrackEnergy_1D_NC->Integral("width")>0.0) hTrackEnergy_1D_NC->Scale(1.0/hTrackEnergy_1D_NC->Integral("width"));
01026 if(hEventEnergy_1D_NC->Integral("width")>0.0) hEventEnergy_1D_NC->Scale(1.0/hEventEnergy_1D_NC->Integral("width"));
01027 if(hNormalization->Integral("width")>0.0) hNormalization->Scale(1.0/hNormalization->Integral("width"));
01028
01029 MSG("MadAb",Msg::kDebug) << endl
01030 << " normalizing histograms [after]: " << endl
01031 << " hCosTheta_1D_CC: " << hCosTheta_1D_CC->Integral("width") << endl
01032 << " hX_1D_CC: " << hX_1D_CC->Integral("width") << endl
01033 << " hY_1D_CC: " << hY_1D_CC->Integral("width") << endl
01034 << " hQ2_1D_CC: " << hQ2_1D_CC->Integral("width") << endl
01035 << " hW2_1D_CC: " << hW2_1D_CC->Integral("width") << endl
01036 << " hCosTheta_2D_CC: " << hCosTheta_2D_CC->Integral("width") << endl
01037 << " hX_2D_CC: " << hX_2D_CC->Integral("width") << endl
01038 << " hY_2D_CC: " << hY_2D_CC->Integral("width") << endl
01039 << " hQ2_2D_CC: " << hQ2_2D_CC->Integral("width") << endl
01040 << " hW2_2D_CC: " << hW2_2D_CC->Integral("width") << endl
01041 << " hE_2D_CC: " << hE_2D_CC->Integral("width") << endl
01042 << " hTrackPHfrac_1D_CC: " << hTrackPHfrac_1D_CC->Integral("width") << endl
01043 << " hTrackPHmean_1D_CC: " << hTrackPHmean_1D_CC->Integral("width") << endl
01044 << " hTrackQPsigmaQP_1D_CC: " << hTrackQPsigmaQP_1D_CC->Integral("width") << endl
01045 << " hTrackLikePlanes_1D_CC: " << hTrackLikePlanes_1D_CC->Integral("width") << endl
01046 << " hTrackPHfrac_2D_CC: " << hTrackPHfrac_2D_CC->Integral("width") << endl
01047 << " hTrackPHmean_2D_CC: " << hTrackPHmean_2D_CC->Integral("width") << endl
01048 << " hTrackQPsigmaQP_2D_CC: " << hTrackQPsigmaQP_2D_CC->Integral("width") << endl
01049 << " hTrackLikePlanes_2D_CC: " << hTrackLikePlanes_2D_CC->Integral("width") << endl
01050 << " hTrackPlanes_2D_CC_1: " << hTrackPlanes_2D_CC_1->Integral("width") << endl
01051 << " hTrackPlanes_2D_CC_2: " << hTrackPlanes_2D_CC_2->Integral("width") << endl
01052 << " hTrackCharge_1D_CC: " << hTrackCharge_1D_CC->Integral("width") << endl
01053 << " hTrackEnergy_1D_CC: " << hTrackEnergy_1D_CC->Integral("width") << endl
01054 << " hEventEnergy_1D_CC: " << hEventEnergy_1D_CC->Integral("width") << endl
01055 << " hCosTheta_1D_NC: " << hCosTheta_1D_NC->Integral("width") << endl
01056 << " hX_1D_NC: " << hX_1D_NC->Integral("width") << endl
01057 << " hY_1D_NC: " << hY_1D_NC->Integral("width") << endl
01058 << " hQ2_1D_NC: " << hQ2_1D_NC->Integral("width") << endl
01059 << " hW2_1D_NC: " << hW2_1D_NC->Integral("width") << endl
01060 << " hCosTheta_2D_NC: " << hCosTheta_2D_NC->Integral("width") << endl
01061 << " hX_2D_NC: " << hX_2D_NC->Integral("width") << endl
01062 << " hY_2D_NC: " << hY_2D_NC->Integral("width") << endl
01063 << " hQ2_2D_NC: " << hQ2_2D_NC->Integral("width") << endl
01064 << " hW2_2D_NC: " << hW2_2D_NC->Integral("width") << endl
01065 << " hE_2D_NC: " << hE_2D_NC->Integral("width") << endl
01066 << " hTrackPHfrac_1D_NC: " << hTrackPHfrac_1D_NC->Integral("width") << endl
01067 << " hTrackPHmean_1D_NC: " << hTrackPHmean_1D_NC->Integral("width") << endl
01068 << " hTrackQPsigmaQP_1D_NC: " << hTrackQPsigmaQP_1D_NC->Integral("width") << endl
01069 << " hTrackLikePlanes_1D_NC: " << hTrackLikePlanes_1D_NC->Integral("width") << endl
01070 << " hTrackPHfrac_2D_NC: " << hTrackPHfrac_2D_NC->Integral("width") << endl
01071 << " hTrackPHmean_2D_NC: " << hTrackPHmean_2D_NC->Integral("width") << endl
01072 << " hTrackQPsigmaQP_2D_NC: " << hTrackQPsigmaQP_2D_NC->Integral("width") << endl
01073 << " hTrackLikePlanes_2D_NC: " << hTrackLikePlanes_2D_NC->Integral("width") << endl
01074 << " hTrackPlanes_2D_NC_1: " << hTrackPlanes_2D_NC_1->Integral("width") << endl
01075 << " hTrackPlanes_2D_NC_2: " << hTrackPlanes_2D_NC_2->Integral("width") << endl
01076 << " hTrackCharge_1D_NC: " << hTrackCharge_1D_NC->Integral("width") << endl
01077 << " hTrackEnergy_1D_NC: " << hTrackEnergy_1D_NC->Integral("width") << endl
01078 << " hEventEnergy_1D_NC: " << hEventEnergy_1D_NC->Integral("width") << endl
01079 << " hNormalization: " << hNormalization->Integral("width") << endl;
01080 }
01081 }
01082
01083 void MadAbID::FillPDFs(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord, Int_t ccnc)
01084 {
01085 this->FillPDFs(ntpEvent,ntpRecord,ccnc,1.0);
01086 }
01087
01088 void MadAbID::FillPDFs(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord, Int_t ccnc, Double_t osc)
01089 {
01090 MSG("MadAb",Msg::kDebug) << " *** MadAbID::FillPDFs(...) *** " << endl;
01091
01092
01093
01094 if( !fFoundPdfs ) this->MakePDFs();
01095
01096
01097
01098 if( fFoundPdfs ){
01099
01100
01101
01102 this->MakeRecoVariables(ntpEvent,ntpRecord);
01103
01104 Int_t bin,binx,biny;
01105 Double_t widthx,widthy;
01106 Double_t area,weight,var;
01107 Double_t p,ptot,mean,sigma;
01108 Int_t min,max,trkplanes;
01109
01110
01111
01112 if( TRKplanes>0 ){
01113 var = RECOdircosneu;
01114 bin = hCosTheta_1D_CC->FindBin(var);
01115 area = hCosTheta_1D_CC->GetBinWidth(bin);
01116
01117 if( area!=0.0 ){
01118 weight = osc/area;
01119 if( ccnc==1 ) hCosTheta_1D_CC->Fill(var,weight);
01120 if( ccnc==0 ) hCosTheta_1D_NC->Fill(var,weight);
01121 }
01122 }
01123
01124
01125
01126 if( TRKplanes>0 ){
01127 var = RECOx;
01128 bin = hX_1D_CC->FindBin(var);
01129 area = hX_1D_CC->GetBinWidth(bin);
01130
01131 if( area!=0.0 ){
01132 weight = osc/area;
01133 if( ccnc==1 ) hX_1D_CC->Fill(var,weight);
01134 if( ccnc==0 ) hX_1D_NC->Fill(var,weight);
01135 }
01136 }
01137
01138
01139
01140 if( TRKplanes>0 ){
01141 var = RECOy;
01142 bin = hY_1D_CC->FindBin(var);
01143 area = hY_1D_CC->GetBinWidth(bin);
01144
01145 if( area!=0.0 ){
01146 weight = osc/area;
01147 if( ccnc==1 ) hY_1D_CC->Fill(var,weight);
01148 if( ccnc==0 ) hY_1D_NC->Fill(var,weight);
01149 }
01150 }
01151
01152
01153
01154 if( TRKplanes>0 ){
01155 var = RECOq2;
01156 bin = hQ2_1D_CC->FindBin(var);
01157 area = hQ2_1D_CC->GetBinWidth(bin);
01158
01159 if( area!=0.0 ){
01160 weight = osc/area;
01161 if( ccnc==1 ) hQ2_1D_CC->Fill(var,weight);
01162 if( ccnc==0 ) hQ2_1D_NC->Fill(var,weight);
01163 }
01164 }
01165
01166
01167
01168 if( TRKplanes>0 ){
01169 var = RECOw2;
01170 bin = hW2_1D_CC->FindBin(var);
01171 area = hW2_1D_CC->GetBinWidth(bin);
01172
01173 if( area!=0.0 ){
01174 weight = osc/area;
01175 if( ccnc==1 ) hW2_1D_CC->Fill(var,weight);
01176 if( ccnc==0 ) hW2_1D_NC->Fill(var,weight);
01177 }
01178 }
01179
01180
01181
01182 if( TRKplanes>0 ){
01183 var = RECOdircosneu;
01184 binx = hCosTheta_2D_CC->GetXaxis()->FindBin(RECOenu);
01185 widthx = hCosTheta_2D_CC->GetXaxis()->GetBinWidth(binx);
01186 biny = hCosTheta_2D_CC->GetYaxis()->FindBin(var);
01187 widthy = hCosTheta_2D_CC->GetYaxis()->GetBinWidth(biny);
01188 area = widthx*widthy;
01189
01190 if( area!=0.0 ){
01191 weight = osc/area;
01192 if( ccnc==1 ) hCosTheta_2D_CC->Fill(RECOenu,var,weight);
01193 if( ccnc==0 ) hCosTheta_2D_NC->Fill(RECOenu,var,weight);
01194 }
01195 }
01196
01197
01198
01199 if( TRKplanes>0 ){
01200 var = RECOx;
01201 binx = hX_2D_CC->GetXaxis()->FindBin(RECOenu);
01202 widthx = hX_2D_CC->GetXaxis()->GetBinWidth(binx);
01203 biny = hX_2D_CC->GetYaxis()->FindBin(var);
01204 widthy = hX_2D_CC->GetYaxis()->GetBinWidth(biny);
01205 area = widthx*widthy;
01206
01207 if( area!=0.0 ){
01208 weight = osc/area;
01209 if( ccnc==1 ) hX_2D_CC->Fill(RECOenu,var,weight);
01210 if( ccnc==0 ) hX_2D_NC->Fill(RECOenu,var,weight);
01211 }
01212 }
01213
01214
01215
01216 if( TRKplanes>0 ){
01217 var = RECOy;
01218 binx = hY_2D_CC->GetXaxis()->FindBin(RECOenu);
01219 widthx = hY_2D_CC->GetXaxis()->GetBinWidth(binx);
01220 biny = hY_2D_CC->GetYaxis()->FindBin(var);
01221 widthy = hY_2D_CC->GetYaxis()->GetBinWidth(biny);
01222 area = widthx*widthy;
01223
01224 if( area!=0.0 ){
01225 weight = osc/area;
01226 if( ccnc==1 ) hY_2D_CC->Fill(RECOenu,var,weight);
01227 if( ccnc==0 ) hY_2D_NC->Fill(RECOenu,var,weight);
01228 }
01229 }
01230
01231
01232
01233 if( TRKplanes>0 ){
01234 var = RECOq2;
01235 binx = hQ2_2D_CC->GetXaxis()->FindBin(RECOenu);
01236 widthx = hQ2_2D_CC->GetXaxis()->GetBinWidth(binx);
01237 biny = hQ2_2D_CC->GetYaxis()->FindBin(var);
01238 widthy = hQ2_2D_CC->GetYaxis()->GetBinWidth(biny);
01239 area = widthx*widthy;
01240
01241 if( area!=0.0 ){
01242 weight = osc/area;
01243 if( ccnc==1 ) hQ2_2D_CC->Fill(RECOenu,var,weight);
01244 if( ccnc==0 ) hQ2_2D_NC->Fill(RECOenu,var,weight);
01245 }
01246 }
01247
01248
01249
01250 if( TRKplanes>0 ){
01251 var = RECOw2;
01252 binx = hW2_2D_CC->GetXaxis()->FindBin(RECOenu);
01253 widthx = hW2_2D_CC->GetXaxis()->GetBinWidth(binx);
01254 biny = hW2_2D_CC->GetYaxis()->FindBin(var);
01255 widthy = hW2_2D_CC->GetYaxis()->GetBinWidth(biny);
01256 area = widthx*widthy;
01257
01258 if( area!=0.0 ){
01259 weight = osc/area;
01260 if( ccnc==1 ) hW2_2D_CC->Fill(RECOenu,var,weight);
01261 if( ccnc==0 ) hW2_2D_NC->Fill(RECOenu,var,weight);
01262 }
01263 }
01264
01265
01266
01267 if( TRKplanes>0 ){
01268 bin = hE_2D_CC->FindBin(RECOenu);
01269 area = hE_2D_CC->GetBinWidth(bin);
01270
01271 if( area!=0.0 ){
01272 weight = osc/area;
01273 if( ccnc==1 ) hE_2D_CC->Fill(RECOenu,weight);
01274 if( ccnc==0 ) hE_2D_NC->Fill(RECOenu,weight);
01275 }
01276 }
01277
01278
01279
01280 if( TRKplanes>0 ){
01281 var = TRKtotph/RECOtotph;
01282 bin = hTrackPHfrac_1D_CC->FindBin(var);
01283 area = hTrackPHfrac_1D_CC->GetBinWidth(bin);
01284
01285 if( area!=0.0 ){
01286 weight = osc/area;
01287 if( ccnc==1 ) hTrackPHfrac_1D_CC->Fill(var,weight);
01288 if( ccnc==0 ) hTrackPHfrac_1D_NC->Fill(var,weight);
01289 }
01290 }
01291
01292
01293
01294 if( TRKplanes>0 ){
01295 if(TRKtrackplanes>2) var = TRKtrkph/TRKtrackplanes;
01296 else var = TRKtotph/TRKplanes;
01297 bin = hTrackPHmean_1D_CC->FindBin(var);
01298 area = hTrackPHmean_1D_CC->GetBinWidth(bin);
01299
01300 if( area!=0.0 ){
01301 weight = osc/area;
01302 if( ccnc==1 ) hTrackPHmean_1D_CC->Fill(var,weight);
01303 if( ccnc==0 ) hTrackPHmean_1D_NC->Fill(var,weight);
01304 }
01305 }
01306
01307
01308
01309 if( TRKplanes>0 ){
01310 var = fabs(TRKFITqpsigmaqp);
01311 bin = hTrackQPsigmaQP_1D_CC->FindBin(var);
01312 area = hTrackQPsigmaQP_1D_CC->GetBinWidth(bin);
01313
01314 if( area!=0.0 ){
01315 weight = osc/area;
01316 if( ccnc==1 ) hTrackQPsigmaQP_1D_CC->Fill(var,weight);
01317 if( ccnc==0 ) hTrackQPsigmaQP_1D_NC->Fill(var,weight);
01318 }
01319 }
01320
01321
01322
01323 if( TRKplanes>0 ){
01324 var = TRKtrackplanes;
01325 bin = hTrackLikePlanes_1D_CC->FindBin(var);
01326 area = hTrackLikePlanes_1D_CC->GetBinWidth(bin);
01327
01328 if( area!=0.0 ){
01329 weight = osc/area;
01330 if( ccnc==1 ) hTrackLikePlanes_1D_CC->Fill(var,weight);
01331 if( ccnc==0 ) hTrackLikePlanes_1D_NC->Fill(var,weight);
01332 }
01333 }
01334
01335
01336
01337 if( TRKplanes>0 ){
01338 trkplanes = TRKplanes;
01339 var = TRKtotph/RECOtotph;
01340
01341 mean = trkplanes;
01342 sigma = sqrt(mean);
01343 min = (Int_t)(mean-2.0*sigma); if(min<5) min=5;
01344 max = (Int_t)(mean+2.0*sigma); if(max>485) max=485;
01345
01346 ptot = 0.0;
01347 for( Int_t pln=min; pln<=max; pln++ ){
01348 p = exp(-0.5*pow((pln-mean)/(sigma),2.0));
01349 ptot += p;
01350 }
01351
01352 for( Int_t pln=min; pln<=max; pln++ ){
01353 binx = hTrackPHfrac_2D_CC->GetXaxis()->FindBin(pln);
01354 widthx = hTrackPHfrac_2D_CC->GetXaxis()->GetBinWidth(binx);
01355 biny = hTrackPHfrac_2D_CC->GetYaxis()->FindBin(var);
01356 widthy = hTrackPHfrac_2D_CC->GetYaxis()->GetBinWidth(biny);
01357 area = widthx*widthy;
01358 if( area!=0.0 ){
01359 p = exp(-0.5*pow((pln-mean)/(sigma),2.0));
01360 weight = (osc/area)*(p/ptot);
01361 if( ccnc==1 ) hTrackPHfrac_2D_CC->Fill(pln,var,weight);
01362 if( ccnc==0 ) hTrackPHfrac_2D_NC->Fill(pln,var,weight);
01363 }
01364 }
01365 }
01366
01367
01368
01369 if( TRKplanes>0 ){
01370 trkplanes = TRKplanes;
01371 if(TRKtrackplanes>2) var = TRKtrkph/TRKtrackplanes;
01372 else var = TRKtotph/TRKplanes;
01373
01374 mean = trkplanes;
01375 sigma = sqrt(mean);
01376 min = (Int_t)(mean-2.0*sigma); if(min<5) min=5;
01377 max = (Int_t)(mean+2.0*sigma); if(max>485) max=485;
01378
01379 ptot = 0.0;
01380 for( Int_t pln=min; pln<=max; pln++ ){
01381 p = exp(-0.5*pow((pln-mean)/(sigma),2.0));
01382 ptot += p;
01383 }
01384
01385 for( Int_t pln=min; pln<=max; pln++ ){
01386 binx = hTrackPHmean_2D_CC->GetXaxis()->FindBin(pln);
01387 widthx = hTrackPHmean_2D_CC->GetXaxis()->GetBinWidth(binx);
01388 biny = hTrackPHmean_2D_CC->GetYaxis()->FindBin(var);
01389 widthy = hTrackPHmean_2D_CC->GetYaxis()->GetBinWidth(biny);
01390 area = widthx*widthy;
01391
01392 if( area!=0.0 ){
01393 p = exp(-0.5*pow((pln-mean)/(sigma),2.0));
01394 weight = (osc/area)*(p/ptot);
01395 if( ccnc==1 ) hTrackPHmean_2D_CC->Fill(pln,var,weight);
01396 if( ccnc==0 ) hTrackPHmean_2D_NC->Fill(pln,var,weight);
01397 }
01398 }
01399 }
01400
01401
01402
01403 if( TRKplanes>0 ){
01404 trkplanes = TRKplanes;
01405 var = fabs(TRKFITqpsigmaqp);
01406
01407 mean = trkplanes;
01408 sigma = sqrt(mean);
01409 min = (Int_t)(mean-2.0*sigma); if(min<5) min=5;
01410 max = (Int_t)(mean+2.0*sigma); if(max>485) max=485;
01411
01412 ptot = 0.0;
01413 for( Int_t pln=min; pln<=max; pln++ ){
01414 p = exp(-0.5*pow((pln-mean)/(sigma),2.0));
01415 ptot += p;
01416 }
01417
01418 for( Int_t pln=min; pln<=max; pln++ ){
01419 binx = hTrackQPsigmaQP_2D_CC->GetXaxis()->FindBin(pln);
01420 widthx = hTrackQPsigmaQP_2D_CC->GetXaxis()->GetBinWidth(binx);
01421 biny = hTrackQPsigmaQP_2D_CC->GetYaxis()->FindBin(var);
01422 widthy = hTrackQPsigmaQP_2D_CC->GetYaxis()->GetBinWidth(biny);
01423 area = widthx*widthy;
01424
01425 if( area!=0.0 ){
01426 p = exp(-0.5*pow((pln-mean)/(sigma),2.0));
01427 weight = (osc/area)*(p/ptot);
01428 if( ccnc==1 ) hTrackQPsigmaQP_2D_CC->Fill(pln,var,weight);
01429 if( ccnc==0 ) hTrackQPsigmaQP_2D_NC->Fill(pln,var,weight);
01430 }
01431 }
01432 }
01433
01434
01435
01436 if( TRKplanes>0 ){
01437 var = TRKtrackplanes;
01438 binx = hTrackLikePlanes_2D_CC->GetXaxis()->FindBin(TRKplanes);
01439 widthx = hTrackLikePlanes_2D_CC->GetXaxis()->GetBinWidth(binx);
01440 biny = hTrackLikePlanes_2D_CC->GetYaxis()->FindBin(var);
01441 widthy = hTrackLikePlanes_2D_CC->GetYaxis()->GetBinWidth(biny);
01442 area = widthx*widthy;
01443
01444 if( area!=0.0 ){
01445 weight = osc/area;
01446 if( ccnc==1 ) hTrackLikePlanes_2D_CC->Fill(TRKplanes,var,weight);
01447 if( ccnc==0 ) hTrackLikePlanes_2D_NC->Fill(TRKplanes,var,weight);
01448 }
01449 }
01450
01451
01452
01453 if( TRKplanes>0 ){
01454 bin = hTrackPlanes_2D_CC_1->FindBin(TRKplanes);
01455 area = hTrackPlanes_2D_CC_1->GetBinWidth(bin);
01456
01457 if( area!=0.0 ){
01458 weight = osc/area;
01459 if( ccnc==1 ) hTrackPlanes_2D_CC_1->Fill(TRKplanes,weight);
01460 if( ccnc==0 ) hTrackPlanes_2D_NC_1->Fill(TRKplanes,weight);
01461 }
01462 }
01463
01464
01465
01466 if( TRKplanes>0 ){
01467 trkplanes = TRKplanes;
01468
01469 mean = trkplanes;
01470 sigma = sqrt(mean);
01471 min = (Int_t)(mean-2.0*sigma); if(min<5) min=5;
01472 max = (Int_t)(mean+2.0*sigma); if(max>485) max=485;
01473
01474 ptot = 0.0;
01475 for( Int_t pln=min; pln<=max; pln++ ){
01476 p = exp(-0.5*pow((pln-mean)/(sigma),2.0));
01477 ptot += p;
01478 }
01479
01480 for( Int_t pln=min; pln<=max; pln++ ){
01481 bin = hTrackPlanes_2D_CC_2->FindBin(pln);
01482 area = hTrackPlanes_2D_CC_2->GetBinWidth(bin);
01483
01484 if( area!=0.0 ){
01485 p = exp(-0.5*pow((pln-mean)/(sigma),2.0));
01486 weight = (osc/area)*(p/ptot);
01487 if( ccnc==1 ) hTrackPlanes_2D_CC_2->Fill(pln,weight);
01488 if( ccnc==0 ) hTrackPlanes_2D_NC_2->Fill(pln,weight);
01489 }
01490 }
01491 }
01492
01493
01494
01495 if( TRKplanes>0 ){
01496 var = TRKFITemcharge;
01497 bin = hTrackCharge_1D_CC->FindBin(var);
01498 area = hTrackCharge_1D_CC->GetBinWidth(bin);
01499
01500 if( area!=0 ){
01501 weight = osc/area;
01502 if( ccnc==1 ) hTrackCharge_1D_CC->Fill(var,weight);
01503 if( ccnc==0 ) hTrackCharge_1D_NC->Fill(var,weight);
01504 }
01505 }
01506
01507
01508
01509 if( TRKplanes>0 ){
01510 var = TRKspanplanes;
01511 bin = hTrackEnergy_1D_CC->FindBin(var);
01512 area = hTrackEnergy_1D_CC->GetBinWidth(bin);
01513
01514 if( area!=0.0 ){
01515 weight = osc/area;
01516 if( ccnc==1 ) hTrackEnergy_1D_CC->Fill(var,weight);
01517 if( ccnc==0 ) hTrackEnergy_1D_NC->Fill(var,weight);
01518 }
01519 }
01520
01521
01522
01523 if( TRKplanes>0 ){
01524 var = RECOenu;
01525 bin = hEventEnergy_1D_CC->FindBin(var);
01526 area = hEventEnergy_1D_CC->GetBinWidth(bin);
01527
01528 if( area!=0.0 ){
01529 weight = osc/area;
01530 if( ccnc==1 ) hEventEnergy_1D_CC->Fill(var,weight);
01531 if( ccnc==0 ) hEventEnergy_1D_NC->Fill(var,weight);
01532 }
01533 }
01534
01535
01536
01537 if( TRKplanes>0 ){
01538 weight = osc;
01539 hNormalization->Fill(ccnc,weight);
01540 }
01541 }
01542
01543 }
01544
01545 void MadAbID::MakePDFs()
01546 {
01547 MSG("MadAb",Msg::kDebug) << " *** MadAbID::MakePDFs() *** " << endl;
01548
01549 fMakePdfs = 1;
01550
01551 if( fMakePdfs ){
01552
01553
01554
01555 if(hCosTheta_1D_CC) delete hCosTheta_1D_CC;
01556 if(hX_1D_CC) delete hX_1D_CC;
01557 if(hY_1D_CC) delete hY_1D_CC;
01558 if(hQ2_1D_CC) delete hQ2_1D_CC;
01559 if(hW2_1D_CC) delete hW2_1D_CC;
01560 if(hCosTheta_2D_CC) delete hCosTheta_2D_CC;
01561 if(hX_2D_CC) delete hX_2D_CC;
01562 if(hY_2D_CC) delete hY_2D_CC;
01563 if(hQ2_2D_CC) delete hQ2_2D_CC;
01564 if(hW2_2D_CC) delete hW2_2D_CC;
01565 if(hE_2D_CC) delete hE_2D_CC;
01566 if(hTrackPHfrac_1D_CC) delete hTrackPHfrac_1D_CC;
01567 if(hTrackPHmean_1D_CC) delete hTrackPHmean_1D_CC;
01568 if(hTrackQPsigmaQP_1D_CC) delete hTrackQPsigmaQP_1D_CC;
01569 if(hTrackLikePlanes_1D_CC) delete hTrackLikePlanes_1D_CC;
01570 if(hTrackPHfrac_2D_CC) delete hTrackPHfrac_2D_CC;
01571 if(hTrackPHmean_2D_CC) delete hTrackPHmean_2D_CC;
01572 if(hTrackQPsigmaQP_2D_CC) delete hTrackQPsigmaQP_2D_CC;
01573 if(hTrackLikePlanes_2D_CC) delete hTrackLikePlanes_2D_CC;
01574 if(hTrackPlanes_2D_CC_1) delete hTrackPlanes_2D_CC_1;
01575 if(hTrackPlanes_2D_CC_2) delete hTrackPlanes_2D_CC_2;
01576 if(hTrackCharge_1D_CC) delete hTrackCharge_1D_CC;
01577 if(hTrackEnergy_1D_CC) delete hTrackEnergy_1D_CC;
01578 if(hEventEnergy_1D_CC) delete hEventEnergy_1D_CC;
01579
01580 if(hCosTheta_1D_NC) delete hCosTheta_1D_NC;
01581 if(hX_1D_NC) delete hX_1D_NC;
01582 if(hY_1D_NC) delete hY_1D_NC;
01583 if(hQ2_1D_NC) delete hQ2_1D_NC;
01584 if(hW2_1D_NC) delete hW2_1D_NC;
01585 if(hCosTheta_2D_NC) delete hCosTheta_2D_NC;
01586 if(hX_2D_NC) delete hX_2D_NC;
01587 if(hY_2D_NC) delete hY_2D_NC;
01588 if(hQ2_2D_NC) delete hQ2_2D_NC;
01589 if(hW2_2D_NC) delete hW2_2D_NC;
01590 if(hE_2D_NC) delete hE_2D_NC;
01591 if(hTrackPHfrac_1D_NC) delete hTrackPHfrac_1D_NC;
01592 if(hTrackPHmean_1D_NC) delete hTrackPHmean_1D_NC;
01593 if(hTrackQPsigmaQP_1D_NC) delete hTrackQPsigmaQP_1D_NC;
01594 if(hTrackLikePlanes_1D_NC) delete hTrackLikePlanes_1D_NC;
01595 if(hTrackPHfrac_2D_NC) delete hTrackPHfrac_2D_NC;
01596 if(hTrackPHmean_2D_NC) delete hTrackPHmean_2D_NC;
01597 if(hTrackQPsigmaQP_2D_NC) delete hTrackQPsigmaQP_2D_NC;
01598 if(hTrackLikePlanes_2D_NC) delete hTrackLikePlanes_2D_NC;
01599 if(hTrackPlanes_2D_NC_1) delete hTrackPlanes_2D_NC_1;
01600 if(hTrackPlanes_2D_NC_2) delete hTrackPlanes_2D_NC_2;
01601 if(hTrackCharge_1D_NC) delete hTrackCharge_1D_NC;
01602 if(hTrackEnergy_1D_NC) delete hTrackEnergy_1D_NC;
01603 if(hEventEnergy_1D_NC) delete hEventEnergy_1D_NC;
01604
01605 if(hNormalization) delete hNormalization;
01606
01607 fFoundPdfs = 0;
01608
01609
01610
01611
01612 Double_t* xcos = new Double_t[9];
01613 xcos[0]=-0.10;
01614 xcos[1]=+0.30;
01615 xcos[2]=+0.50;
01616 xcos[3]=+0.65;
01617 xcos[4]=+0.76;
01618 xcos[5]=+0.85;
01619 xcos[6]=+0.92;
01620 xcos[7]=+0.97;
01621 xcos[8]=+1.00;
01622
01623 hCosTheta_1D_CC = new TH1D("hCosTheta_1D_CC","hCosTheta_1D_CC",8,xcos);
01624 hCosTheta_1D_NC = new TH1D("hCosTheta_1D_NC","hCosTheta_1D_NC",8,xcos);
01625 hCosTheta_1D_CC->SetLineColor(2);
01626 hCosTheta_1D_CC->SetLineWidth(2);
01627 hCosTheta_1D_NC->SetLineColor(4);
01628 hCosTheta_1D_NC->SetLineWidth(2);
01629
01630 hX_1D_CC = new TH1D("hX_1D_CC","hX_1D_CC",1,0.0,1.0);
01631 hX_1D_NC = new TH1D("hX_1D_NC","hX_1D_NC",1,0.0,1.0);
01632 hX_1D_CC->SetLineColor(2);
01633 hX_1D_CC->SetLineWidth(2);
01634 hX_1D_NC->SetLineColor(4);
01635 hX_1D_NC->SetLineWidth(2);
01636
01637 hY_1D_CC = new TH1D("hY_1D_CC","hY_1D_CC",50,0.0,1.0);
01638 hY_1D_NC = new TH1D("hY_1D_NC","hY_1D_NC",50,0.0,1.0);
01639 hY_1D_CC->SetLineColor(2);
01640 hY_1D_CC->SetLineWidth(2);
01641 hY_1D_NC->SetLineColor(4);
01642 hY_1D_NC->SetLineWidth(2);
01643
01644 hQ2_1D_CC = new TH1D("hQ2_1D_CC","hQ2_1D_CC",1,0.0,1000.0);
01645 hQ2_1D_NC = new TH1D("hQ2_1D_NC","hQ2_1D_NC",1,0.0,1000.0);
01646 hQ2_1D_CC->SetLineColor(2);
01647 hQ2_1D_CC->SetLineWidth(2);
01648 hQ2_1D_NC->SetLineColor(4);
01649 hQ2_1D_NC->SetLineWidth(2);
01650
01651 hW2_1D_CC = new TH1D("hW2_1D_CC","hW2_1D_CC",1,0.0,1000.0);
01652 hW2_1D_NC = new TH1D("hW2_1D_NC","hW2_1D_NC",1,0.0,1000.0);
01653 hW2_1D_CC->SetLineColor(2);
01654 hW2_1D_CC->SetLineWidth(2);
01655 hW2_1D_NC->SetLineColor(4);
01656 hW2_1D_NC->SetLineWidth(2);
01657
01658 hTrackPHfrac_1D_CC = new TH1D("hTrackPHfrac_1D_CC","hTrackPHfrac_1D_CC",50,0.001,1.001);
01659 hTrackPHfrac_1D_NC = new TH1D("hTrackPHfrac_1D_NC","hTrackPHfrac_1D_NC",50,0.001,1.001);
01660 hTrackPHfrac_1D_CC->SetLineColor(2);
01661 hTrackPHfrac_1D_CC->SetLineWidth(2);
01662 hTrackPHfrac_1D_NC->SetLineColor(4);
01663 hTrackPHfrac_1D_NC->SetLineWidth(2);
01664
01665 Double_t* xph = new Double_t[37];
01666 xph[0] = 0.0;
01667 xph[1] = 100.0;
01668 xph[2] = 150.0;
01669 for(Int_t i=3; i<31; i++){
01670 xph[i] = 180.0 + 20.0*(i-3.0);
01671 }
01672 xph[31] = 750.0;
01673 xph[32] = 800.0;
01674 xph[33] = 1000.0;
01675 xph[34] = 2000.0;
01676 xph[35] = 5000.0;
01677 xph[36] = 10000.0;
01678
01679 hTrackPHmean_1D_CC = new TH1D("hTrackPHmean_1D_CC","hTrackPHmean_1D_CC",36,xph);
01680 hTrackPHmean_1D_NC = new TH1D("hTrackPHmean_1D_NC","hTrackPHmean_1D_NC",36,xph);
01681 hTrackPHmean_1D_CC->SetLineColor(2);
01682 hTrackPHmean_1D_CC->SetLineWidth(2);
01683 hTrackPHmean_1D_NC->SetLineColor(4);
01684 hTrackPHmean_1D_NC->SetLineWidth(2);
01685
01686 Double_t* xqp = new Double_t[46];
01687
01688 for(Int_t i=0; i<20; i++){
01689 xqp[i] = 0.0 + 1.0*(i-0.0);
01690 }
01691
01692 for(Int_t i=20; i<35; i++){
01693 xqp[i] = +20.0 + 2.0*(i-20.0);
01694 }
01695
01696 for(Int_t i=35; i<42; i++){
01697 xqp[i] = +50.0 + 5.0*(i-35.0);
01698 }
01699
01700 xqp[42] = +90.0;
01701 xqp[43] = +100.0;
01702 xqp[44] = +250.0;
01703 xqp[45] = +1000.0;
01704
01705 hTrackQPsigmaQP_1D_CC = new TH1D("hTrackQPsigmaQP_1D_CC","hTrackQPsigmaQP_1D_CC",45,xqp);
01706 hTrackQPsigmaQP_1D_NC = new TH1D("hTrackQPsigmaQP_1D_NC","hTrackQPsigmaQP_1D_NC",45,xqp);
01707 hTrackQPsigmaQP_1D_CC->SetLineColor(2);
01708 hTrackQPsigmaQP_1D_CC->SetLineWidth(2);
01709 hTrackQPsigmaQP_1D_NC->SetLineColor(4);
01710 hTrackQPsigmaQP_1D_NC->SetLineWidth(2);
01711
01712 Double_t* xtrk = new Double_t[151];
01713
01714 for(Int_t i=0; i<100; i++){
01715 xtrk[i] = -0.5+0.0+1.0*(i-0.0);
01716 }
01717
01718 for(Int_t i=100; i<120; i++){
01719 xtrk[i] = -0.5+100.0+5.0*(i-100.0);
01720 }
01721
01722 for(Int_t i=120; i<150; i++){
01723 xtrk[i] = -0.5+200.0+10.0*(i-120.0);
01724 }
01725
01726 xtrk[150] = -0.5+500.0;
01727
01728 hTrackLikePlanes_1D_CC = new TH1D("hTrackLikePlanes_1D_CC","hTrackLikePlanes_1D_CC",150,xtrk);
01729 hTrackLikePlanes_1D_NC = new TH1D("hTrackLikePlanes_1D_NC","hTrackLikePlanes_1D_NC",150,xtrk);
01730
01731 hTrackLikePlanes_1D_CC->SetLineColor(2);
01732 hTrackLikePlanes_1D_CC->SetLineWidth(2);
01733 hTrackLikePlanes_1D_NC->SetLineColor(4);
01734 hTrackLikePlanes_1D_NC->SetLineWidth(2);
01735
01736 hTrackCharge_1D_CC = new TH1D("hTrackCharge_1D_CC","hTrackCharge_1D_CC",3,-1.5,+1.5);
01737 hTrackCharge_1D_NC = new TH1D("hTrackCharge_1D_NC","hTrackCharge_1D_NC",3,-1.5,+1.5);
01738 hTrackCharge_1D_CC->SetLineColor(2);
01739 hTrackCharge_1D_CC->SetLineWidth(2);
01740 hTrackCharge_1D_NC->SetLineColor(4);
01741 hTrackCharge_1D_NC->SetLineWidth(2);
01742
01743 hTrackEnergy_1D_CC = new TH1D("hTrackEnergy_1D_CC","hTrackEnergy_1D_CC",150,xtrk);
01744 hTrackEnergy_1D_NC = new TH1D("hTrackEnergy_1D_NC","hTrackEnergy_1D_NC",150,xtrk);
01745 hTrackEnergy_1D_CC->SetLineColor(2);
01746 hTrackEnergy_1D_CC->SetLineWidth(2);
01747 hTrackEnergy_1D_NC->SetLineColor(4);
01748 hTrackEnergy_1D_NC->SetLineWidth(2);
01749
01750 Double_t* xevt = new Double_t[53];
01751
01752 for(Int_t i=0; i<20; i++){
01753 xevt[i] = 0.0 + 0.5*(i-0.0);
01754 }
01755
01756 for(Int_t i=20; i<30; i++){
01757 xevt[i] = 10.0 + 1.0*(i-20.0);
01758 }
01759
01760 for(Int_t i=30; i<45; i++){
01761 xevt[i] = 20.0 + 2.0*(i-30.0);
01762 }
01763
01764 for(Int_t i=45; i<50; i++){
01765 xevt[i] = 50.0 + 5.0*(i-45.0);
01766 }
01767
01768 xevt[50] = 80.0;
01769 xevt[51] = 100.0;
01770 xevt[52] = 1000.0;
01771
01772 hEventEnergy_1D_CC = new TH1D("hEventEnergy_1D_CC","hEventEnergy_1D_CC",52,xevt);
01773 hEventEnergy_1D_NC = new TH1D("hEventEnergy_1D_NC","hEventEnergy_1D_NC",52,xevt);
01774 hEventEnergy_1D_CC->SetLineColor(2);
01775 hEventEnergy_1D_CC->SetLineWidth(2);
01776 hEventEnergy_1D_NC->SetLineColor(4);
01777 hEventEnergy_1D_NC->SetLineWidth(2);
01778
01779 hNormalization = new TH1D("hNormalization","hNormalization",2,-0.5,1.5);
01780
01781
01782
01783
01784 Double_t* xenu = new Double_t[27];
01785
01786 for(Int_t i=0; i<10; i++){
01787 xenu[i] = 0.0 + 1.0*(i-0.0);
01788 }
01789
01790 for(Int_t i=10; i<15; i++){
01791 xenu[i] = 10.0 + 2.0*(i-10.0);
01792 }
01793
01794 for(Int_t i=15; i<21; i++){
01795 xenu[i] = 20.0 + 5.0*(i-15.0);
01796 }
01797
01798 for(Int_t i=21; i<24; i++){
01799 xenu[i] = 50.0 + 10.0*(i-21.0);
01800 }
01801
01802 xenu[24] = 80.0;
01803 xenu[25] = 100.0;
01804 xenu[26] = 1000.0;
01805
01806 hCosTheta_2D_CC = new TH2D("hCosTheta_2D_CC","hCosTheta_2D_CC",26,xenu,8,xcos);
01807 hCosTheta_2D_NC = new TH2D("hCosTheta_2D_NC","hCosTheta_2D_NC",26,xenu,8,xcos);
01808 hCosTheta_2D_CC->SetMarkerColor(2);
01809 hCosTheta_2D_NC->SetMarkerColor(4);
01810
01811 hX_2D_CC = new TH2D("hX_2D_CC","hX_2D_CC",26,xenu,1,0.0,1.0);
01812 hX_2D_NC = new TH2D("hX_2D_NC","hX_2D_NC",26,xenu,1,0.0,1.0);
01813 hX_2D_CC->SetMarkerColor(2);
01814 hX_2D_NC->SetMarkerColor(4);
01815
01816 hY_2D_CC = new TH2D("hY_2D_CC","hY_2D_CC",26,xenu,25,0.0,1.0);
01817 hY_2D_NC = new TH2D("hY_2D_NC","hY_2D_NC",26,xenu,25,0.0,1.0);
01818 hY_2D_CC->SetMarkerColor(2);
01819 hY_2D_NC->SetMarkerColor(4);
01820
01821 hQ2_2D_CC = new TH2D("hQ2_2D_CC","hQ2_2D_CC",26,xenu,1,0.0,1000.0);
01822 hQ2_2D_NC = new TH2D("hQ2_2D_NC","hQ2_2D_NC",26,xenu,1,0.0,1000.0);
01823 hQ2_2D_CC->SetMarkerColor(2);
01824 hQ2_2D_NC->SetMarkerColor(4);
01825
01826 hW2_2D_CC = new TH2D("hW2_2D_CC","hW2_2D_CC",26,xenu,1,0.0,1000.0);
01827 hW2_2D_NC = new TH2D("hW2_2D_NC","hW2_2D_NC",26,xenu,1,0.0,1000.0);
01828 hW2_2D_CC->SetMarkerColor(2);
01829 hW2_2D_NC->SetMarkerColor(4);
01830
01831 hE_2D_CC = new TH1D("hE_2D_CC","hE_2D_CC",26,xenu);
01832 hE_2D_NC = new TH1D("hE_2D_NC","hE_2D_NC",26,xenu);
01833 hE_2D_CC->SetLineColor(2);
01834 hE_2D_CC->SetLineWidth(2);
01835 hE_2D_NC->SetLineColor(4);
01836 hE_2D_NC->SetLineWidth(2);
01837
01838 Double_t* xemu = new Double_t[61];
01839
01840 for(Int_t i=0; i<25; i++){
01841 xemu[i] = -0.5+0.0+2.0*(i-0.0);
01842 }
01843
01844 for(Int_t i=25; i<35; i++){
01845 xemu[i] = -0.5+50.0+5.0*(i-25.0);
01846 }
01847
01848 for(Int_t i=35; i<45; i++){
01849 xemu[i] = -0.5+100.0+10.0*(i-35.0);
01850 }
01851
01852 for(Int_t i=45; i<61; i++){
01853 xemu[i] = -0.5+200.0+20.0*(i-45.0);
01854 }
01855
01856 Double_t* xemu2 = new Double_t[5];
01857 xemu2[0] = -0.5;
01858 xemu2[1] = 19.5;
01859 xemu2[2] = 39.5;
01860 xemu2[3] = 59.5;
01861 xemu2[4] = 499.5;
01862
01863 hTrackPHfrac_2D_CC = new TH2D("hTrackPHfrac_2D_CC","hTrackPHfrac_2D_CC",60,xemu,25,0.001,1.001);
01864 hTrackPHfrac_2D_NC = new TH2D("hTrackPHfrac_2D_NC","hTrackPHfrac_2D_NC",60,xemu,25,0.001,1.001);
01865 hTrackPHfrac_2D_CC->SetMarkerColor(2);
01866 hTrackPHfrac_2D_NC->SetMarkerColor(4);
01867
01868 Double_t* xph2 = new Double_t[22];
01869
01870 xph2[0] = 0.0;
01871 xph2[1] = 100.0;
01872
01873 for(Int_t i=2; i<15; i++){
01874 xph2[i] = 150.0 + 50.0*(i-2.0);
01875 }
01876
01877 xph2[15] = 800.0;
01878 xph2[16] = 900.0;
01879 xph2[17] = 1000.0;
01880 xph2[18] = 1500.0;
01881 xph2[19] = 2500.0;
01882 xph2[20] = 5000.0;
01883 xph2[21] = 10000.0;
01884
01885 hTrackPHmean_2D_CC = new TH2D("hTrackPHmean_2D_CC","hTrackPHmean_2D_CC",60,xemu,21,xph2);
01886 hTrackPHmean_2D_NC = new TH2D("hTrackPHmean_2D_NC","hTrackPHmean_2D_NC",60,xemu,21,xph2);
01887 hTrackPHmean_2D_CC->SetMarkerColor(2);
01888 hTrackPHmean_2D_NC->SetMarkerColor(4);
01889
01890 Double_t* xqp2 = new Double_t[35];
01891
01892 for(Int_t i=0; i<25; i++){
01893 xqp2[i] = 0.0 + 2.0*(i-0.0);
01894 }
01895
01896 for(Int_t i=25; i<30; i++){
01897 xqp2[i] = +50.0 + 5.0*(i-25.0);
01898 }
01899
01900 xqp2[30] = +80.0;
01901 xqp2[31] = +90.0;
01902 xqp2[32] = +100.0;
01903 xqp2[33] = +250.0;
01904 xqp2[34] = +1000.0;
01905
01906 hTrackQPsigmaQP_2D_CC = new TH2D("hTrackQPsigmaQP_2D_CC","hTrackQPsigmaQP_2D_CC",60,xemu,34,xqp2);
01907 hTrackQPsigmaQP_2D_NC = new TH2D("hTrackQPsigmaQP_2D_NC","hTrackQPsigmaQP_2D_NC",60,xemu,34,xqp2);
01908 hTrackQPsigmaQP_2D_CC->SetMarkerColor(2);
01909 hTrackQPsigmaQP_2D_NC->SetMarkerColor(4);
01910
01911 hTrackLikePlanes_2D_CC = new TH2D("hTrackLikePlanes_2D_CC","hTrackLikePlanes_2D_CC",60,xemu,60,xemu);
01912 hTrackLikePlanes_2D_NC = new TH2D("hTrackLikePlanes_2D_NC","hTrackLikePlanes_2D_NC",60,xemu,60,xemu);
01913 hTrackLikePlanes_2D_CC->SetMarkerColor(2);
01914 hTrackLikePlanes_2D_NC->SetMarkerColor(4);
01915
01916 hTrackPlanes_2D_CC_1 = new TH1D("hTrackPlanes_2D_CC_1","hTrackPlanes_2D_CC_1",60,xemu);
01917 hTrackPlanes_2D_NC_1 = new TH1D("hTrackPlanes_2D_NC_1","hTrackPlanes_2D_NC_1",60,xemu);
01918 hTrackPlanes_2D_CC_1->SetLineColor(2);
01919 hTrackPlanes_2D_CC_1->SetLineWidth(2);
01920 hTrackPlanes_2D_NC_1->SetLineColor(4);
01921 hTrackPlanes_2D_NC_1->SetLineWidth(2);
01922
01923 hTrackPlanes_2D_CC_2 = new TH1D("hTrackPlanes_2D_CC_2","hTrackPlanes_2D_CC_2",60,xemu);
01924 hTrackPlanes_2D_NC_2 = new TH1D("hTrackPlanes_2D_NC_2","hTrackPlanes_2D_NC_2",60,xemu);
01925 hTrackPlanes_2D_CC_2->SetLineColor(2);
01926 hTrackPlanes_2D_CC_2->SetLineWidth(2);
01927 hTrackPlanes_2D_NC_2->SetLineColor(4);
01928 hTrackPlanes_2D_NC_2->SetLineWidth(2);
01929
01930 delete [] xcos;
01931 delete [] xph;
01932 delete [] xqp;
01933 delete [] xtrk;
01934 delete [] xevt;
01935 delete [] xenu;
01936 delete [] xemu;
01937 delete [] xemu2;
01938 delete [] xph2;
01939 delete [] xqp2;
01940
01941 fFoundPdfs = 1;
01942 }
01943 }
01944
01945 void MadAbID::MakePidVariables(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
01946 {
01947
01948
01949
01950 const RecCandHeader* ntpHeader = &(ntpRecord->GetHeader());
01951 VldContext vldc = ntpHeader->GetVldContext();
01952
01953 Int_t newRun = ntpHeader->GetRun();
01954 Int_t newSnarl = ntpHeader->GetSnarl();
01955 Int_t newEvent = ntpEvent->index;
01956
01957 if( newRun==PidRun
01958 && newSnarl==PidSnarl
01959 && newEvent==PidEvent ){
01960 return;
01961 }
01962
01963
01964
01965
01966 MSG("MadAb",Msg::kDebug) << " *** MadAbPid::MakePidVariables(...) *** " << endl;
01967 MSG("MadAb",Msg::kDebug) << " Run=" << newRun << " Snarl=" << newSnarl << " Event=" << newEvent << endl;
01968
01969 PidRun = newRun;
01970 PidSnarl = newSnarl;
01971 PidEvent = newEvent;
01972
01973 MSG("MadAb",Msg::kVerbose) << endl
01974 << " PID ANALYSIS MODE: " << endl
01975 << " CosTheta: [" << fCosTheta_1D << "," << fCosTheta_2D << "]" << endl
01976 << " X: [" << fX_1D << "," << fX_2D << "]" << endl
01977 << " Y: [" << fY_1D << "," << fY_2D << "]" << endl
01978 << " Q2: [" << fQ2_1D << "," << fQ2_2D << "]" << endl
01979 << " W2: [" << fW2_1D << "," << fW2_2D << "]" << endl
01980 << " TrackPHfrac: [" << fTrackPHfrac_1D << "," << fTrackPHfrac_2D << "]" << endl
01981 << " TrackPHmean: [" << fTrackPHmean_1D << "," << fTrackPHmean_2D << "]" << endl
01982 << " TrackQPsigmaQP: [" << fTrackQPsigmaQP_1D << "," << fTrackQPsigmaQP_2D << "]" << endl
01983 << " TrackLikePlanes: [" << fTrackLikePlanes_1D << "," << fTrackLikePlanes_2D << "]" << endl
01984 << " TrackCharge: [" << fTrackCharge_1D << "]" << endl
01985 << " TrackEnergy: [" << fTrackEnergy_1D << "]" << endl
01986 << " EventEnergy: [" << fEventEnergy_1D << "]" << endl
01987 << " Normalization: [" << fNormalization << "]" << endl;
01988
01989
01990
01991 if( !fFoundPdfs ) this->MakePDFs();
01992
01993
01994
01995 PIDPASSFAIL = 0;
01996 PROBCC = 0.0;
01997 PROBNC = 0.0;
01998 PID = -10.0;
01999
02000 PID_COSTHETA_PROBCC = 0.0;
02001 PID_X_PROBCC = 0.0;
02002 PID_Y_PROBCC = 0.0;
02003 PID_Q2_PROBCC = 0.0;
02004 PID_W2_PROBCC = 0.0;
02005 PID_TRACKPHFRAC_PROBCC = 0.0;
02006 PID_TRACKPHMEAN_PROBCC = 0.0;
02007 PID_TRACKQPSIGMAQP_PROBCC = 0.0;
02008 PID_TRACKLIKEPLANES_PROBCC = 0.0;
02009 PID_TRACKCHARGE_PROBCC = 0.0;
02010 PID_TRACKENERGY_PROBCC = 0.0;
02011 PID_EVENTENERGY_PROBCC = 0.0;
02012 PID_NORMALIZATION_PROBCC = 0.0;
02013
02014 PID_COSTHETA_PROBNC = 0.0;
02015 PID_X_PROBNC = 0.0;
02016 PID_Y_PROBNC = 0.0;
02017 PID_Q2_PROBNC = 0.0;
02018 PID_W2_PROBNC = 0.0;
02019 PID_TRACKPHFRAC_PROBNC = 0.0;
02020 PID_TRACKPHMEAN_PROBNC = 0.0;
02021 PID_TRACKQPSIGMAQP_PROBNC = 0.0;
02022 PID_TRACKLIKEPLANES_PROBNC = 0.0;
02023 PID_TRACKCHARGE_PROBNC = 0.0;
02024 PID_TRACKENERGY_PROBNC = 0.0;
02025 PID_EVENTENERGY_PROBNC = 0.0;
02026 PID_NORMALIZATION_PROBNC = 0.0;
02027
02028 if( fFoundPdfs ){
02029
02030
02031
02032 this->MakeRecoVariables(ntpEvent,ntpRecord);
02033
02034 Double_t var = 0.0;
02035 Double_t probmin = 0.0;
02036 Double_t probcc = 0.0;
02037 Double_t probnc = 0.0;
02038 Double_t probCC = -999.9;
02039 Double_t probNC = -999.9;
02040 Double_t tempprobCC = 0.0;
02041 Double_t tempprobNC = 0.0;
02042
02043
02044
02045 if( fCosTheta_1D ){
02046 probcc = 0.0; probnc = 0.0;
02047 if( TRKplanes>0 ){
02048 var = RECOdircosneu;
02049 probcc = hCosTheta_1D_CC->GetBinContent(hCosTheta_1D_CC->FindBin(var));
02050 probnc = hCosTheta_1D_NC->GetBinContent(hCosTheta_1D_NC->FindBin(var));
02051 }
02052 if( probcc>probmin || probnc>probmin ){
02053 if( probCC<0.0 ) probCC = 1.0;
02054 if( probNC<0.0 ) probNC = 1.0;
02055 probCC *= probcc; probNC *= probnc;
02056 PID_COSTHETA_PROBCC = probcc; PID_COSTHETA_PROBNC = probnc;
02057 }
02058 MSG("MadAb",Msg::kVerbose) << " CosTheta_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02059 }
02060
02061
02062
02063 if( fX_1D ){
02064 probcc = 0.0; probnc = 0.0;
02065 if( TRKplanes>0 ){
02066 var = RECOx;
02067 probcc = hX_1D_CC->GetBinContent(hX_1D_CC->FindBin(var));
02068 probnc = hX_1D_NC->GetBinContent(hX_1D_NC->FindBin(var));
02069 }
02070 if( probcc>probmin || probnc>probmin ){
02071 if( probCC<0.0 ) probCC = 1.0;
02072 if( probNC<0.0 ) probNC = 1.0;
02073 probCC *= probcc; probNC *= probnc;
02074 PID_X_PROBCC = probcc; PID_X_PROBNC = probnc;
02075 }
02076 MSG("MadAb",Msg::kVerbose) << " X_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02077 }
02078
02079
02080
02081 if( fY_1D ){
02082 probcc = 0.0; probnc = 0.0;
02083 if( TRKplanes>0 ){
02084 var = RECOy;
02085 probcc = hY_1D_CC->GetBinContent(hY_1D_CC->FindBin(var));
02086 probnc = hY_1D_NC->GetBinContent(hY_1D_NC->FindBin(var));
02087 }
02088 if( probcc>probmin || probnc>probmin ){
02089 if( probCC<0.0 ) probCC = 1.0;
02090 if( probNC<0.0 ) probNC = 1.0;
02091 probCC *= probcc; probNC *= probnc;
02092 PID_Y_PROBCC = probcc; PID_Y_PROBNC = probnc;
02093 }
02094 MSG("MadAb",Msg::kVerbose) << " Y_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02095 }
02096
02097
02098
02099 if( fQ2_1D ){
02100 probcc = 0.0; probnc = 0.0;
02101 if( TRKplanes>0 ){
02102 var = RECOq2;
02103 probcc = hQ2_1D_CC->GetBinContent(hQ2_1D_CC->FindBin(var));
02104 probnc = hQ2_1D_NC->GetBinContent(hQ2_1D_NC->FindBin(var));
02105 }
02106 if( probcc>probmin || probnc>probmin ){
02107 if( probCC<0.0 ) probCC = 1.0;
02108 if( probNC<0.0 ) probNC = 1.0;
02109 probCC *= probcc; probNC *= probnc;
02110 PID_Q2_PROBCC = probcc; PID_Q2_PROBNC = probnc;
02111 }
02112 MSG("MadAb",Msg::kVerbose) << " Q2_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02113 }
02114
02115
02116
02117 if( fW2_1D ){
02118 probcc = 0.0; probnc = 0.0;
02119 if( TRKplanes>0 ){
02120 var = RECOw2;
02121 probcc = hW2_1D_CC->GetBinContent(hW2_1D_CC->FindBin(var));
02122 probnc = hW2_1D_NC->GetBinContent(hW2_1D_NC->FindBin(var));
02123 }
02124 if( probcc>probmin || probnc>probmin ){
02125 if( probCC<0.0 ) probCC = 1.0;
02126 if( probNC<0.0 ) probNC = 1.0;
02127 probCC *= probcc; probNC *= probnc;
02128 PID_W2_PROBCC = probcc; PID_W2_PROBNC = probnc;
02129 }
02130 MSG("MadAb",Msg::kVerbose) << " W2_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02131 }
02132
02133
02134
02135 if( fCosTheta_2D ){
02136 probcc = 0.0; probnc = 0.0;
02137 if( TRKplanes>0 ){
02138 var = RECOdircosneu;
02139 if( hE_2D_CC->GetBinContent(hE_2D_CC->FindBin(RECOenu))>0.0 )
02140 probcc = hCosTheta_2D_CC->GetBinContent(hCosTheta_2D_CC->FindBin(RECOenu,var))/
02141 hE_2D_CC->GetBinContent(hE_2D_CC->FindBin(RECOenu));
02142 if( hE_2D_NC->GetBinContent(hE_2D_NC->FindBin(RECOenu))>0.0 )
02143 probnc = hCosTheta_2D_NC->GetBinContent(hCosTheta_2D_NC->FindBin(RECOenu,var))/
02144 hE_2D_NC->GetBinContent(hE_2D_NC->FindBin(RECOenu));
02145 }
02146 if( probcc>probmin || probnc>probmin ){
02147 if( probCC<0.0 ) probCC = 1.0;
02148 if( probNC<0.0 ) probNC = 1.0;
02149 probCC *= probcc; probNC *= probnc;
02150 PID_COSTHETA_PROBCC = probcc; PID_COSTHETA_PROBNC = probnc;
02151 }
02152 MSG("MadAb",Msg::kVerbose) << " CosTheta_2D: pCC=" << probcc << " pNC=" << probnc << endl;
02153 }
02154
02155
02156
02157 if( fX_2D ){
02158 probcc = 0.0; probnc = 0.0;
02159 if( TRKplanes>0 ){
02160 var = RECOx;
02161 if( hE_2D_CC->GetBinContent(hE_2D_CC->FindBin(RECOenu))>0.0 )
02162 probcc = hX_2D_CC->GetBinContent(hX_2D_CC->FindBin(RECOenu,var))/
02163 hE_2D_CC->GetBinContent(hE_2D_CC->FindBin(RECOenu));
02164 if( hE_2D_NC->GetBinContent(hE_2D_NC->FindBin(RECOenu))>0.0 )
02165 probnc = hX_2D_NC->GetBinContent(hX_2D_NC->FindBin(RECOenu,var))/
02166 hE_2D_NC->GetBinContent(hE_2D_NC->FindBin(RECOenu));
02167 }
02168 if( probcc>probmin || probnc>probmin ){
02169 if( probCC<0.0 ) probCC = 1.0;
02170 if( probNC<0.0 ) probNC = 1.0;
02171 probCC *= probcc; probNC *= probnc;
02172 PID_X_PROBCC = probcc; PID_X_PROBNC = probnc;
02173 }
02174 MSG("MadAb",Msg::kVerbose) << " X_2D: pCC=" << probcc << " pNC=" << probnc << endl;
02175 }
02176
02177
02178
02179 if( fY_2D ){
02180 probcc = 0.0; probnc = 0.0;
02181 if( TRKplanes>0 ){
02182 var = RECOy;
02183 if( hE_2D_CC->GetBinContent(hE_2D_CC->FindBin(RECOenu))>0.0 )
02184 probcc = hY_2D_CC->GetBinContent(hY_2D_CC->FindBin(RECOenu,var))/
02185 hE_2D_CC->GetBinContent(hE_2D_CC->FindBin(RECOenu));
02186 if( hE_2D_NC->GetBinContent(hE_2D_NC->FindBin(RECOenu))>0.0 )
02187 probnc = hY_2D_NC->GetBinContent(hY_2D_NC->FindBin(RECOenu,var))/
02188 hE_2D_NC->GetBinContent(hE_2D_NC->FindBin(RECOenu));
02189 }
02190 if( probcc>probmin || probnc>probmin ){
02191 if( probCC<0.0 ) probCC = 1.0;
02192 if( probNC<0.0 ) probNC = 1.0;
02193 probCC *= probcc; probNC *= probnc;
02194 PID_Y_PROBCC = probcc; PID_Y_PROBNC = probnc;
02195 }
02196 MSG("MadAb",Msg::kVerbose) << " Y_2D: pCC=" << probcc << " pNC=" << probnc << endl;
02197 }
02198
02199
02200
02201 if( fQ2_2D ){
02202 probcc = 0.0; probnc = 0.0;
02203 if( TRKplanes>0 ){
02204 var = RECOq2;
02205 if( hE_2D_CC->GetBinContent(hE_2D_CC->FindBin(RECOenu))>0.0 )
02206 probcc = hQ2_2D_CC->GetBinContent(hQ2_2D_CC->FindBin(RECOenu,var))/
02207 hE_2D_CC->GetBinContent(hE_2D_CC->FindBin(RECOenu));
02208 if( hE_2D_NC->GetBinContent(hE_2D_NC->FindBin(RECOenu))>0.0 )
02209 probnc = hQ2_2D_NC->GetBinContent(hQ2_2D_NC->FindBin(RECOenu,var))/
02210 hE_2D_NC->GetBinContent(hE_2D_NC->FindBin(RECOenu));
02211 }
02212 if( probcc>probmin || probnc>probmin ){
02213 if( probCC<0.0 ) probCC = 1.0;
02214 if( probNC<0.0 ) probNC = 1.0;
02215 probCC *= probcc; probNC *= probnc;
02216 PID_Q2_PROBCC = probcc; PID_Q2_PROBNC = probnc;
02217 }
02218 MSG("MadAb",Msg::kVerbose) << " Q2_2D: pCC=" << probcc << " pNC=" << probnc << endl;
02219 }
02220
02221
02222
02223 if( fW2_2D ){
02224 probcc = 0.0; probnc = 0.0;
02225 if( TRKplanes>0 ){
02226 var = RECOw2;
02227 if( hE_2D_CC->GetBinContent(hE_2D_CC->FindBin(RECOenu))>0.0 )
02228 probcc = hW2_2D_CC->GetBinContent(hQ2_2D_CC->FindBin(RECOenu,var))/
02229 hE_2D_CC->GetBinContent(hE_2D_CC->FindBin(RECOenu));
02230 if( hE_2D_NC->GetBinContent(hE_2D_NC->FindBin(RECOenu))>0.0 )
02231 probnc = hW2_2D_NC->GetBinContent(hQ2_2D_NC->FindBin(RECOenu,var))/
02232 hE_2D_NC->GetBinContent(hE_2D_NC->FindBin(RECOenu));
02233 }
02234 if( probcc>probmin || probnc>probmin ){
02235 if( probCC<0.0 ) probCC = 1.0;
02236 if( probNC<0.0 ) probNC = 1.0;
02237 probCC *= probcc; probNC *= probnc;
02238 PID_W2_PROBCC = probcc; PID_W2_PROBNC = probnc;
02239 }
02240 MSG("MadAb",Msg::kVerbose) << " W2_2D: pCC=" << probcc << " pNC=" << probnc << endl;
02241 }
02242
02243
02244
02245 if( fTrackPHfrac_1D ){
02246 probcc = 0.0; probnc = 0.0;
02247 if( TRKplanes>0 ){
02248 var = TRKtotph/RECOtotph;
02249 probcc = hTrackPHfrac_1D_CC->GetBinContent(hTrackPHfrac_1D_CC->FindBin(var));
02250 probnc = hTrackPHfrac_1D_NC->GetBinContent(hTrackPHfrac_1D_NC->FindBin(var));
02251 }
02252 if( probcc>probmin || probnc>probmin ){
02253 if( probCC<0.0 ) probCC = 1.0;
02254 if( probNC<0.0 ) probNC = 1.0;
02255 probCC *= probcc; probNC *= probnc;
02256 PID_TRACKPHFRAC_PROBCC = probcc; PID_TRACKPHFRAC_PROBNC = probnc;
02257 }
02258 MSG("MadAb",Msg::kVerbose) << " TrackPHfrac_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02259 }
02260
02261
02262
02263 if( fTrackPHmean_1D ){
02264 probcc = 0.0; probnc = 0.0;
02265 if( TRKplanes>0 ){
02266 if( TRKtrackplanes>2 ) var = TRKtrkph/TRKtrackplanes;
02267 else var = TRKtotph/TRKplanes;
02268 probcc = hTrackPHmean_1D_CC->GetBinContent(hTrackPHmean_1D_CC->FindBin(var));
02269 probnc = hTrackPHmean_1D_NC->GetBinContent(hTrackPHmean_1D_NC->FindBin(var));
02270 }
02271 if( probcc>probmin || probnc>probmin ){
02272 if( probCC<0.0 ) probCC = 1.0;
02273 if( probNC<0.0 ) probNC = 1.0;
02274 probCC *= probcc; probNC *= probnc;
02275 PID_TRACKPHMEAN_PROBCC = probcc; PID_TRACKPHMEAN_PROBNC = probnc;
02276 }
02277 MSG("MadAb",Msg::kVerbose) << " TrackPHmean_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02278 }
02279
02280
02281
02282 if( fTrackQPsigmaQP_1D ){
02283 probcc = 0.0; probnc = 0.0;
02284 if( TRKplanes>0 ){
02285 var = fabs(TRKFITqpsigmaqp);
02286 probcc = hTrackQPsigmaQP_1D_CC->GetBinContent(hTrackQPsigmaQP_1D_CC->FindBin(var));
02287 probnc = hTrackQPsigmaQP_1D_NC->GetBinContent(hTrackQPsigmaQP_1D_NC->FindBin(var));
02288 }
02289 if( probcc>probmin || probnc>probmin ){
02290 if( probCC<0.0 ) probCC = 1.0;
02291 if( probNC<0.0 ) probNC = 1.0;
02292 probCC *= probcc; probNC *= probnc;
02293 PID_TRACKQPSIGMAQP_PROBCC = probcc; PID_TRACKQPSIGMAQP_PROBNC = probnc;
02294 }
02295 MSG("MadAb",Msg::kVerbose) << " QPsigmaQP_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02296 }
02297
02298
02299
02300 if( fTrackLikePlanes_1D ){
02301 probcc = 0.0; probnc = 0.0;
02302 if( TRKplanes>0 ){
02303 var = TRKtrackplanes;
02304 probcc = hTrackLikePlanes_1D_CC->GetBinContent(hTrackLikePlanes_1D_CC->FindBin(var));
02305 probnc = hTrackLikePlanes_1D_NC->GetBinContent(hTrackLikePlanes_1D_NC->FindBin(var));
02306 }
02307 if( probcc>probmin || probnc>probmin ){
02308 if( probCC<0.0 ) probCC = 1.0;
02309 if( probNC<0.0 ) probNC = 1.0;
02310 probCC *= probcc; probNC *= probnc;
02311 PID_TRACKLIKEPLANES_PROBCC = probcc; PID_TRACKLIKEPLANES_PROBNC = probnc;
02312 }
02313 MSG("MadAb",Msg::kVerbose) << " TrackLikePlanes_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02314 }
02315
02316
02317
02318 if( fTrackPHfrac_2D ){
02319 probcc = 0.0; probnc = 0.0;
02320 if( TRKplanes>0.0 ){
02321 var = TRKtotph/RECOtotph;
02322 if( hTrackPlanes_2D_CC_2->GetBinContent(hTrackPlanes_2D_CC_2->FindBin(TRKplanes))>0.0 )
02323 probcc = hTrackPHfrac_2D_CC->GetBinContent(hTrackPHfrac_2D_CC->FindBin(TRKplanes,var))/
02324 hTrackPlanes_2D_CC_2->GetBinContent(hTrackPlanes_2D_CC_2->FindBin(TRKplanes));
02325 if( hTrackPlanes_2D_NC_2->GetBinContent(hTrackPlanes_2D_NC_2->FindBin(TRKplanes))>0.0 )
02326 probnc = hTrackPHfrac_2D_NC->GetBinContent(hTrackPHfrac_2D_NC->FindBin(TRKplanes,var))/
02327 hTrackPlanes_2D_NC_2->GetBinContent(hTrackPlanes_2D_NC_2->FindBin(TRKplanes));
02328 }
02329 if( probcc>probmin || probnc>probmin ){
02330 if( probCC<0.0 ) probCC = 1.0;
02331 if( probNC<0.0 ) probNC = 1.0;
02332 probCC *= probcc; probNC *= probnc;
02333 PID_TRACKPHFRAC_PROBCC = probcc; PID_TRACKPHFRAC_PROBNC = probnc;
02334 }
02335 MSG("MadAb",Msg::kVerbose) << " TrackPHfrac_2D: pCC=" << probcc << " pNC=" << probnc << endl;
02336 }
02337
02338
02339
02340 if( fTrackPHmean_2D ){
02341 probcc = 0.0; probnc = 0.0;
02342 if( TRKplanes>0 ){
02343 if( TRKtrackplanes>2 ) var = TRKtrkph/TRKtrackplanes;
02344 else var = TRKtotph/TRKplanes;
02345 if( hTrackPlanes_2D_CC_2->GetBinContent(hTrackPlanes_2D_CC_2->FindBin(TRKplanes))>0.0 )
02346 probcc = hTrackPHmean_2D_CC->GetBinContent(hTrackPHmean_2D_CC->FindBin(TRKplanes,var))/
02347 hTrackPlanes_2D_CC_2->GetBinContent(hTrackPlanes_2D_CC_2->FindBin(TRKplanes));
02348 if( hTrackPlanes_2D_NC_2->GetBinContent(hTrackPlanes_2D_NC_2->FindBin(TRKplanes))>0.0 )
02349 probnc = hTrackPHmean_2D_NC->GetBinContent(hTrackPHmean_2D_NC->FindBin(TRKplanes,var))/
02350 hTrackPlanes_2D_NC_2->GetBinContent(hTrackPlanes_2D_NC_2->FindBin(TRKplanes));
02351 }
02352 if( probcc>probmin || probnc>probmin ){
02353 if( probCC<0.0 ) probCC = 1.0;
02354 if( probNC<0.0 ) probNC = 1.0;
02355 probCC *= probcc; probNC *= probnc;
02356 PID_TRACKPHMEAN_PROBCC = probcc; PID_TRACKPHMEAN_PROBNC = probnc;
02357 }
02358 MSG("MadAb",Msg::kVerbose) << " TrackPHmean_2D: pCC=" << probcc << " pNC=" << probnc << endl;
02359 }
02360
02361
02362
02363 if( fTrackQPsigmaQP_2D ){
02364 probcc = 0.0; probnc = 0.0;
02365 if( TRKplanes>0 ){
02366 var = fabs(TRKFITqpsigmaqp);
02367 if( hTrackPlanes_2D_CC_2->GetBinContent(hTrackPlanes_2D_CC_2->FindBin(TRKplanes))>0.0 )
02368 probcc = hTrackQPsigmaQP_2D_CC->GetBinContent(hTrackQPsigmaQP_2D_CC->FindBin(TRKplanes,var))/
02369 hTrackPlanes_2D_CC_2->GetBinContent(hTrackPlanes_2D_CC_2->FindBin(TRKplanes));
02370 if( hTrackPlanes_2D_NC_2->GetBinContent(hTrackPlanes_2D_NC_2->FindBin(TRKplanes))>0.0 )
02371 probnc = hTrackQPsigmaQP_2D_NC->GetBinContent(hTrackQPsigmaQP_2D_NC->FindBin(TRKplanes,var))/
02372 hTrackPlanes_2D_NC_2->GetBinContent(hTrackPlanes_2D_NC_2->FindBin(TRKplanes));
02373 }
02374 if( probcc>probmin || probnc>probmin ){
02375 if( probCC<0.0 ) probCC = 1.0;
02376 if( probNC<0.0 ) probNC = 1.0;
02377 probCC *= probcc; probNC *= probnc;
02378 PID_TRACKQPSIGMAQP_PROBCC = probcc; PID_TRACKQPSIGMAQP_PROBNC = probnc;
02379 }
02380 MSG("MadAb",Msg::kVerbose) << " QPsigmaQP_2D: pCC=" << probcc << " pNC=" << probnc << endl;
02381 }
02382
02383
02384
02385 if( fTrackLikePlanes_2D ){
02386 probcc = 0.0; probnc = 0.0;
02387 if( TRKplanes>0 ){
02388 var = TRKtrackplanes;
02389 if( hTrackPlanes_2D_CC_1->GetBinContent(hTrackPlanes_2D_CC_1->FindBin(TRKplanes))>0.0 )
02390 probcc = hTrackLikePlanes_2D_CC->GetBinContent(hTrackLikePlanes_2D_CC->FindBin(TRKplanes,var))/
02391 hTrackPlanes_2D_CC_1->GetBinContent(hTrackPlanes_2D_CC_1->FindBin(TRKplanes));
02392 if( hTrackPlanes_2D_NC_1->GetBinContent(hTrackPlanes_2D_NC_1->FindBin(TRKplanes))>0.0 )
02393 probnc = hTrackLikePlanes_2D_NC->GetBinContent(hTrackLikePlanes_2D_NC->FindBin(TRKplanes,var))/
02394 hTrackPlanes_2D_NC_1->GetBinContent(hTrackPlanes_2D_NC_1->FindBin(TRKplanes));
02395 }
02396 if( probcc>probmin || probnc>probmin ){
02397 if( probCC<0.0 ) probCC = 1.0;
02398 if( probNC<0.0 ) probNC = 1.0;
02399 probCC *= probcc; probNC *= probnc;
02400 PID_TRACKLIKEPLANES_PROBCC = probcc; PID_TRACKLIKEPLANES_PROBNC = probnc;
02401 }
02402 MSG("MadAb",Msg::kVerbose) << " TrackLikePlanes_2D: pCC=" << probcc << " pNC=" << probnc << endl;
02403 }
02404
02405
02406
02407 if( fTrackCharge_1D ){
02408 probcc = 0.0; probnc = 0.0;
02409 if( TRKplanes>0 ){
02410 var = TRKFITemcharge;
02411 probcc = hTrackCharge_1D_CC->GetBinContent(hTrackCharge_1D_CC->FindBin(var));
02412 probnc = hTrackCharge_1D_NC->GetBinContent(hTrackCharge_1D_NC->FindBin(var));
02413 }
02414 if( probcc>probmin || probnc>probmin ){
02415 if( probCC<0.0 ) probCC = 1.0;
02416 if( probNC<0.0 ) probNC = 1.0;
02417 probCC *= probcc; probNC *= probnc;
02418 PID_TRACKCHARGE_PROBCC = probcc; PID_TRACKCHARGE_PROBNC = probnc;
02419 }
02420 MSG("MadAb",Msg::kVerbose) << " TrackCharge_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02421 }
02422
02423
02424
02425 if( fTrackEnergy_1D ){
02426 probcc = 0.0; probnc = 0.0;
02427 if( TRKplanes>0 ){
02428 probcc = hTrackEnergy_1D_CC->GetBinContent(hTrackEnergy_1D_CC->FindBin(TRKspanplanes));
02429 probnc = hTrackEnergy_1D_NC->GetBinContent(hTrackEnergy_1D_NC->FindBin(TRKspanplanes));
02430 }
02431 if( probcc>probmin || probnc>probmin ){
02432 if( probCC<0.0 ) probCC = 1.0;
02433 if( probNC<0.0 ) probNC = 1.0;
02434 probCC *= probcc; probNC *= probnc;
02435 PID_TRACKENERGY_PROBCC = probcc; PID_TRACKENERGY_PROBNC = probnc;
02436 }
02437 MSG("MadAb",Msg::kVerbose) << " TrackEnergy_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02438 }
02439
02440
02441
02442 if( fEventEnergy_1D ){
02443 probcc = 0.0; probnc = 0.0;
02444 if( TRKplanes>0 ){
02445 probcc = hEventEnergy_1D_CC->GetBinContent(hEventEnergy_1D_CC->FindBin(RECOenu));
02446 probnc = hEventEnergy_1D_NC->GetBinContent(hEventEnergy_1D_NC->FindBin(RECOenu));
02447 }
02448 if( probcc>probmin || probnc>probmin ){
02449 if( probCC<0.0 ) probCC = 1.0;
02450 if( probNC<0.0 ) probNC = 1.0;
02451 probCC *= probcc; probNC *= probnc;
02452 PID_EVENTENERGY_PROBCC = probcc; PID_EVENTENERGY_PROBNC = probnc;
02453 }
02454 MSG("MadAb",Msg::kVerbose) << " EventEnergy_1D: pCC=" << probcc << " pNC=" << probnc << endl;
02455 }
02456
02457
02458
02459 if( fNormalization ){
02460 probcc = 0.0; probnc = 0.0;
02461 if( TRKplanes>0 ){
02462 probcc = hNormalization->GetBinContent(hNormalization->FindBin(1));
02463 probnc = hNormalization->GetBinContent(hNormalization->FindBin(0));
02464 }
02465 if( probcc>probmin || probnc>probmin ){
02466 if( probCC<0.0 ) probCC = 1.0;
02467 if( probNC<0.0 ) probNC = 1.0;
02468 probCC *= probcc; probNC *= probnc;
02469 PID_NORMALIZATION_PROBCC = probcc; PID_NORMALIZATION_PROBNC = probnc;
02470 }
02471 MSG("MadAb",Msg::kVerbose) << " Normalization: pCC=" << probcc << " pNC=" << probnc << endl;
02472 }
02473
02474
02475
02476 if( probCC>=0 && probNC>=0 ){
02477 if( probCC+probNC>1.0 ){
02478 tempprobCC = (probCC)/(probCC+probNC);
02479 tempprobNC = (probNC)/(probCC+probNC);
02480 probCC = tempprobCC; probNC = tempprobNC;
02481 }
02482 PROBCC = probCC; PROBNC = probNC;
02483 }
02484
02485 MSG("MadAb",Msg::kVerbose) << " Overall: pCC=" << PROBCC << " pNC=" << PROBNC << endl;
02486
02487
02488
02489 if( PROBCC+PROBNC>0.0 ){
02490 PID = PROBCC/(PROBCC+PROBNC);
02491 PIDPASSFAIL = 1;
02492 }
02493 }
02494
02495 MSG("MadAb",Msg::kDebug) << " PASSFAIL = " << PIDPASSFAIL << endl;
02496 MSG("MadAb",Msg::kDebug) << " PID = " << PID << endl;
02497
02498 MSG("MadAb",Msg::kVerbose) << " *** MadAbID::MakePidVariables(...) FINISHED *** " << endl;
02499
02500 }
02501
02502 void MadAbID::MakeRecoVariables(const NtpSREvent* ntpEvent, const NtpStRecord* ntpRecord)
02503 {
02504
02505
02506
02507 const RecCandHeader* ntpHeader = &(ntpRecord->GetHeader());
02508 VldContext vldc = ntpHeader->GetVldContext();
02509
02510 Int_t newRun = ntpHeader->GetRun();
02511 Int_t newSnarl = ntpHeader->GetSnarl();
02512 Int_t newEvent = ntpEvent->index;
02513
02514 if( newRun==RecoRun
02515 && newSnarl==RecoSnarl
02516 && newEvent==RecoEvent ){
02517 return;
02518 }
02519
02520
02521
02522
02523 MSG("MadAb",Msg::kDebug) << " *** MadAbPid::MakeRecoVariables(...) *** " << endl;
02524 MSG("MadAb",Msg::kDebug) << " Run=" << newRun << " Snarl=" << newSnarl << " Event=" << newEvent << endl;
02525
02526 RecoRun = newRun;
02527 RecoSnarl = newSnarl;
02528 RecoEvent = newEvent;
02529
02530 Ndigits = 0;
02531 Nstrips = 0;
02532 Ntracks = 0;
02533 Nshowers = 0;
02534 RECOminplane = -999;
02535 RECOmaxplane = -999;
02536 RECOtotph = 0.0;
02537
02538 RECOenu = 0.0;
02539 RECOemu = 0.0;
02540 RECOehad = 0.0;
02541 RECOx = 0.0;
02542 RECOy = 0.0;
02543 RECOq2 = 0.0;
02544 RECOw2 = 0.0;
02545 RECOdircosU = 0.0;
02546 RECOdircosV = 0.0;
02547 RECOdircosX = 0.0;
02548 RECOdircosY = 0.0;
02549 RECOdircosZ = 0.0;
02550 RECOdircosneu = 0.0;
02551 RECOQEenu = 0.0;
02552 RECOQEemu = 0.0;
02553 RECOQEehad = 0.0;
02554 RECOQEx = 0.0;
02555 RECOQEy = 0.0;
02556 RECOQEq2 = 0.0;
02557 RECOQEw2 = 0.0;
02558 RECOvtxplane = -999;
02559 RECOvtxU = -999.9;
02560 RECOvtxV = -999.9;
02561 RECOvtxX = -999.9;
02562 RECOvtxY = -999.9;
02563 RECOvtxR = -999.9;
02564 RECOvtxZ = -999.9;
02565 RECOvtxDistToEdge = -999.9;
02566 RECOvtxDistToEndBack = -999.9;
02567 RECOvtxDistToEndForward = -999.9;
02568 RECOcontained = 0;
02569 RECOvtxcontained = 0;
02570 RECOendcontained = 0;
02571
02572 TRKvtxplane = -999;
02573 TRKvtxU = -999.9;
02574 TRKvtxV = -999.9;
02575 TRKvtxX = -999.9;
02576 TRKvtxY = -999.9;
02577 TRKvtxR = -999.9;
02578 TRKvtxZ = -999.9;
02579 TRKvtxdircosU = 0.0;
02580 TRKvtxdircosV = 0.0;
02581 TRKvtxdircosX = 0.0;
02582 TRKvtxdircosY = 0.0;
02583 TRKvtxdircosZ = 0.0;
02584 TRKvtxDistToEdge = -999.9;
02585 TRKvtxDistToEndBack = -999.9;
02586 TRKvtxDistToEndForward = -999.9;
02587 TRKvtxtrace = -999.9;
02588 TRKvtxtraceZ = -999.9;
02589 TRKendplane = -999;
02590 TRKendU = -999.9;
02591 TRKendV = -999.9;
02592 TRKendX = -999.9;
02593 TRKendY = -999.9;
02594 TRKendR = -999.9;
02595 TRKendZ = -999.9;
02596 TRKenddircosU = 0.0;
02597 TRKenddircosV = 0.0;
02598 TRKenddircosX = 0.0;
02599 TRKenddircosY = 0.0;
02600 TRKenddircosZ = 0.0;
02601 TRKendDistToEdge = -999.9;
02602 TRKendDistToEndBack = -999.9;
02603 TRKendDistToEndForward = -999.9;
02604 TRKendtrace = -999.9;
02605 TRKendtraceZ = -999.9;
02606 TRKstrips = 0;
02607 TRKplanes = 0;
02608 TRKtotph = 0.0;
02609 TRKtrkph = 0.0;
02610 TRKshwph = 0.0;
02611 TRKspanplanes = 0;
02612 TRKtrackplanes = 0;
02613 TRKshowerplanes = 0;
02614 TRKtracklikeplanes = 0;
02615 TRKshowerlikeplanes = 0;
02616 TRKcontained = 0;
02617 TRKmomentumRange = 0.0;
02618 TRKmomentumCurve = 0.0;
02619 TRKforwardRMS = 0.0;
02620 TRKforwardNDOF = 0;
02621 TRKbackwardRMS = 0.0;
02622 TRKbackwardNDOF = 0;
02623 TRKFITchi2 = -999.9;
02624 TRKFITndof = 0;
02625 TRKFITpass = 0;
02626 TRKFITqpsigmaqp = 0.0;
02627 TRKFITemcharge = 0;
02628
02629 SHWenergy = 0.0;
02630 SHWcontained = 0;
02631
02632 if( ntpRecord ){
02633
02634
02635
02636 const RecCandHeader* ntpHeader = &(ntpRecord->GetHeader());
02637 VldContext vldc = ntpHeader->GetVldContext();
02638
02639 TClonesArray* stripArray = (TClonesArray*)(ntpRecord->stp);
02640 TClonesArray* showerArray = (TClonesArray*)(ntpRecord->shw);
02641 TClonesArray* trackArray = (TClonesArray*)(ntpRecord->trk);
02642
02643 Ndigits = 0;
02644 Nstrips = 0;
02645 Ntracks = 0;
02646 Nshowers = 0;
02647
02648 Int_t ctr = 0;
02649 Int_t last = 1+stripArray->GetLast();
02650
02651 Int_t* flag = new Int_t[500];
02652 Int_t* shwflag = new Int_t[last];
02653 Int_t* trkflag = new Int_t[last];
02654
02655 for( Int_t n=0; n<last; n++ ){
02656 shwflag[n]=0;
02657 trkflag[n]=0;
02658 }
02659
02660 for( Int_t n=0; n<500; n++ ){
02661 flag[n]=0;
02662 fStripList[n].Clear();
02663 fTrackStripList[n].Clear();
02664 fTrkShowerStripList[n].Clear();
02665 fShwShowerStripList[n].Clear();
02666 }
02667
02668
02669
02670 Double_t temp,dr,dzm,dzp;
02671 RECOminplane = -999;
02672 RECOmaxplane = -999;
02673 RECOtotph = 0.0;
02674
02675
02676 ctr = 0;
02677 Int_t* evtstp = ntpEvent->stp;
02678 for( Int_t k=0; k<ntpEvent->nstrip; k++ ){
02679 Int_t sptr = evtstp[k];
02680 if( sptr>=0 ){
02681 NtpSRStrip* strip = (NtpSRStrip*)(stripArray->At(sptr));
02682 Int_t plane = strip->plane;
02683 if( plane>=0 && plane<500 ){
02684 fStripList[plane].Add(strip);
02685 if( RECOminplane<0 || plane<RECOminplane ) RECOminplane=plane;
02686 if( RECOmaxplane<0 || plane>RECOmaxplane ) RECOmaxplane=plane;
02687 RECOtotph += strip->ph0.sigcor+strip->ph1.sigcor;
02688 Ndigits += strip->ndigit;
02689 Nstrips++;
02690 ctr++;
02691 }
02692 }
02693 }
02694 MSG("MadAb",Msg::kDebug) << " found " << ctr << " strips *** " << endl;
02695
02696
02697
02698 NtpSRVertex vertex = ntpEvent->vtx;
02699 RECOvtxplane = vertex.plane;
02700 RECOvtxU = vertex.u;
02701 RECOvtxV = vertex.v;
02702 RECOvtxX = vertex.x;
02703 RECOvtxY = vertex.y;
02704 RECOvtxZ = vertex.z;
02705
02706
02707
02708 Ntracks = ntpEvent->ntrack;
02709
02710 if( ntpEvent->ntrack>0 ){
02711 MSG("MadAb",Msg::kDebug) << " found reconstructed track(s) *** " << endl;
02712
02713
02714
02715
02716
02717
02718 Int_t trkptr = -1;
02719 Int_t trackplanes = 0;
02720 Int_t* trk = ntpEvent->trk;
02721 for( Int_t n=0; n<ntpEvent->ntrack; n++ ){
02722 Int_t ptr=trk[n];
02723 if( ptr>=0 ){
02724 NtpSRTrack* track = (NtpSRTrack*)(trackArray->At(ptr));
02725 if( abs(track->end.plane-track->vtx.plane)>trackplanes ){
02726 trackplanes=abs(track->end.plane-track->vtx.plane);
02727 trkptr=ptr;
02728 }
02729 }
02730 }
02731
02732
02733 if( trkptr>=0 ){
02734 NtpSRTrack* track = (NtpSRTrack*)(trackArray->At(trkptr));
02735
02736 Double_t p = track->momentum.range;
02737 Double_t qp = track->momentum.qp;
02738 Double_t eqp = track->momentum.eqp;
02739 Double_t chi2 = track->fit.chi2;
02740 Int_t ndof = track->fit.ndof;
02741 Int_t pass = track->fit.pass;
02742 if( qp==0.0 || eqp==0.0 || ndof<=1 ) pass = 0;
02743
02744 TRKmomentumRange = p;
02745 if( qp!=0.0 ) TRKmomentumCurve = fabs(1/qp);
02746 TRKFITchi2 = chi2;
02747 TRKFITndof = ndof;
02748 TRKFITpass = pass;
02749 if( pass ) TRKFITqpsigmaqp = qp/eqp;
02750 if( pass ) TRKFITemcharge = (Int_t)(qp/fabs(qp));
02751
02752 TRKvtxplane = track->vtx.plane;
02753 TRKvtxU = track->vtx.u;
02754 TRKvtxV = track->vtx.v;
02755 TRKvtxX = track->vtx.x;
02756 TRKvtxY = track->vtx.y;
02757 TRKvtxZ = track->vtx.z;
02758
02759 if( vldc.GetDetector()==Detector::kFar ){
02760 TRKvtxR = sqrt(TRKvtxX*TRKvtxX+TRKvtxY*TRKvtxY);
02761 }
02762 if( vldc.GetDetector()==Detector::kNear ){
02763 TRKvtxR = sqrt( (TRKvtxX-1.4828)*(TRKvtxX-1.4828)
02764 + (TRKvtxY-0.2384)*(TRKvtxY-0.2384) );
02765 }
02766
02767 TRKvtxdircosU = track->vtx.dcosu;
02768 TRKvtxdircosV = track->vtx.dcosv;
02769 TRKvtxdircosX = track->vtx.dcosx;
02770 TRKvtxdircosY = track->vtx.dcosy;
02771 TRKvtxdircosZ = track->vtx.dcosz;
02772
02773 TRKendplane = track->end.plane;
02774 TRKendU = track->end.u;
02775 TRKendV = track->end.v;
02776 TRKendX = track->end.x;
02777 TRKendY = track->end.y;
02778 TRKendR = sqrt(TRKendX*TRKendX+TRKendY*TRKendY);
02779 TRKendZ = track->end.z;
02780
02781 TRKenddircosU = track->end.dcosu;
02782 TRKenddircosV = track->end.dcosv;
02783 TRKenddircosX = track->end.dcosx;
02784 TRKenddircosY = track->end.dcosy;
02785 TRKenddircosZ = track->end.dcosz;
02786
02787 TRKcontained = track->contained;
02788
02789 TRKforwardRMS = track->time.forwardRMS;
02790 TRKforwardNDOF = track->time.forwardNDOF;
02791 TRKbackwardRMS = track->time.backwardRMS;
02792 TRKbackwardNDOF = track->time.backwardNDOF;
02793
02794 TRKstrips = 0;
02795
02796
02797 ctr = 0;
02798 Int_t* trkstp = track->stp;
02799 for( Int_t k=0; k<track->nstrip; k++ ){
02800 Int_t sptr = trkstp[k];
02801 if( sptr>=0 ){
02802 NtpSRStrip* strip = (NtpSRStrip*)(stripArray->At(sptr));
02803 Int_t plane = strip->plane;
02804 if( plane>=0 && plane<500 ){
02805 fTrackStripList[plane].Add(strip);
02806 TRKstrips++;
02807 ctr++;
02808 }
02809 if( sptr<last ){
02810 trkflag[sptr] = 1;
02811 }
02812 }
02813 }
02814 MSG("MadAb",Msg::kDebug) << " found " << ctr << " track strips *** " << endl;
02815
02816
02817 if( vldc.GetDetector()==Detector::kFar ){
02818
02819
02820 Double_t pm;
02821 Double_t gradu,gradv,grads;
02822 Double_t dU,dV,dX,dY,dR;
02823 Double_t posU,posV,posX,posY,posZ;
02824 Double_t dirU,dirV,dirX,dirY,dirZ;
02825 Double_t trace,traceZ;
02826 Int_t trkflag;
02827
02828
02829 Double_t vtxU = TRKvtxU;
02830 Double_t vtxV = TRKvtxV;
02831 Double_t vtxX = TRKvtxX;
02832 Double_t vtxY = TRKvtxY;
02833 Double_t vtxZ = TRKvtxZ;
02834
02835 dr = 999.9;
02836 temp=4.0-vtxU; if(temp<dr) dr=temp;
02837 temp=4.0+vtxU; if(temp<dr) dr=temp;
02838 temp=4.0-vtxV; if(temp<dr) dr=temp;
02839 temp=4.0+vtxV; if(temp<dr) dr=temp;
02840 temp=4.0-vtxX; if(temp<dr) dr=temp;
02841 temp=4.0+vtxX; if(temp<dr) dr=temp;
02842 temp=4.0-vtxY; if(temp<dr) dr=temp;
02843 temp=4.0+vtxY; if(temp<dr) dr=temp;
02844 if( dr<999.9 ) TRKvtxDistToEdge = dr;
02845
02846 if( TRKvtxdircosZ>=0.0 ){
02847 if( TRKvtxplane<249 ){ dzm = vtxZ-0.110; dzp = 14.793-vtxZ; }
02848 if( TRKvtxplane>249 ){ dzm = vtxZ-15.986; dzp = 29.955-vtxZ; }
02849 }
02850 if( TRKvtxdircosZ<0.0 ){
02851 if( TRKvtxplane<249 ){ dzm = 14.793-vtxZ; dzp = vtxZ-0.110; }
02852 if( TRKvtxplane>249 ){ dzm = 29.955-vtxZ; dzp = vtxZ-15.986; }
02853 }
02854
02855 TRKvtxDistToEndBack = dzm;
02856 TRKvtxDistToEndForward = dzp;
02857
02858
02859 posU=vtxU; posV=vtxV; posZ=vtxZ;
02860 posX=0.7071*(posU-posV); posY=0.7071*(posU+posV);
02861
02862 gradu=(TRKvtxdircosU/TRKvtxdircosZ);
02863 gradv=(TRKvtxdircosV/TRKvtxdircosZ);
02864 grads=sqrt(gradu*gradu+gradv*gradv+1.0);
02865 if(TRKvtxdircosZ>0.0) pm=-1.0; else pm=+1.0;
02866 dirU=pm*gradu/grads; dirV=pm*gradv/grads; dirZ=pm*1.0/grads;
02867 dirX=0.7071*(dirU-dirV); dirY=0.7071*(dirU+dirV);
02868 trkflag=0; dR=0.0;
02869
02870 if(trkflag==0){
02871 if(dirX>0.0){
02872 dX = 4.0-posX; dY = dX*(dirY/dirX);
02873 if( posY+dY<4.0*(sqrt(2.0)-1.0) && posY+dY>-4.0*(sqrt(2.0)-1.0) ){
02874 trkflag=1; if(dX>0.0) dR=sqrt(dX*dX+dY*dY); if(dX<0.0) dR=-sqrt(dX*dX+dY*dY);
02875 }
02876 }
02877 if(dirX<0.0){
02878 dX = 4.0+posX; dY = -dX*(dirY/dirX);
02879 if( posY+dY<4.0*(sqrt(2.0)-1.0) && posY+dY>-4.0*(sqrt(2.0)-1.0) ){
02880 trkflag=1; if(dX>0.0) dR=sqrt(dX*dX+dY*dY); if(dX<0.0) dR=-sqrt(dX*dX+dY*dY);
02881 }
02882 }
02883 }
02884
02885 if(trkflag==0){
02886 if(dirY>0.0){
02887 dY = 4.0-posY; dX = dY*(dirX/dirY);
02888 if( posX+dX<4.0*(sqrt(2.0)-1.0) && posX+dX>-4.0*(sqrt(2.0)-1.0) ){
02889 trkflag=1; if(dY>0.0) dR=sqrt(dX*dX+dY*dY); if(dY<0.0) dR=-sqrt(dX*dX+dY*dY);
02890 }
02891 }
02892 if(dirY<0.0){
02893 dY = 4.0+posY; dX = -dY*(dirX/dirY);
02894 if( posX+dX<4.0*(sqrt(2.0)-1.0) && posX+dX>-4.0*(sqrt(2.0)-1.0) ){
02895 trkflag=1; if(dY>0.0) dR=sqrt(dX*dX+dY*dY); if(dY<0.0) dR=-sqrt(dX*dX+dY*dY);
02896 }
02897 }
02898 }
02899
02900 if(trkflag==0){
02901 if(dirU>0.0){
02902 dU = 4.0-posU; dV = dU*(dirV/dirU);
02903 if( posV+dV<4.0*(sqrt(2.0)-1.0) && posV+dV>-4.0*(sqrt(2.0)-1.0) ){
02904 trkflag=1; if(dU>0.0) dR=sqrt(dV*dV+dU*dU); if(dU<0.0) dR=-sqrt(dV*dV+dU*dU);
02905 }
02906 }
02907 if(dirU<0.0){
02908 dU = 4.0+posU; dV = -dU*(dirV/dirU);
02909 if( posV+dV<4.0*(sqrt(2.0)-1.0) && posV+dV>-4.0*(sqrt(2.0)-1.0) ){
02910 trkflag=1; if(dU>0.0) dR=sqrt(dV*dV+dU*dU); if(dU<0.0) dR=-sqrt(dV*dV+dU*dU);
02911 }
02912 }
02913 }
02914
02915 if(trkflag==0){
02916 if(dirV>0.0){
02917 dV = 4.0-posV; dU = dV*(dirU/dirV);
02918 if( posU+dU<4.0*(sqrt(2.0)-1.0) && posU+dU>-4.0*(sqrt(2.0)-1.0) ){
02919 trkflag=1; if(dV>0.0) dR=sqrt(dV*dV+dU*dU); if(dV<0.0) dR=-sqrt(dV*dV+dU*dU);
02920 }
02921 }
02922 if(dirV<0.0){
02923 dV = 4.0+posV; dU = -dV*(dirU/dirV);
02924 if( posU+dU<4.0*(sqrt(2.0)-1.0) && posU+dU>-4.0*(sqrt(2.0)-1.0) ){
02925 trkflag=1; if(dV>0.0) dR=sqrt(dV*dV+dU*dU); if(dV<0.0) dR=-sqrt(dV*dV+dU*dU);
02926 }
02927 }
02928 }
02929
02930 traceZ=0.0; trace=0.0;
02931 if(sqrt(dirX*dirX+dirY*dirY)>0.0){
02932 traceZ=(pm*dR*dirZ)/sqrt(dirX*dirX+dirY*dirY);
02933 }
02934 trace=(traceZ)/(pm*dirZ);
02935
02936 TRKvtxtrace = trace;
02937 TRKvtxtraceZ = traceZ;
02938
02939
02940 Double_t endU = TRKendU;
02941 Double_t endV = TRKendV;
02942 Double_t endX = TRKendX;
02943 Double_t endY = TRKendY;
02944 Double_t endZ = TRKendZ;
02945
02946 dr = 999.9;
02947 temp=4.0-endU; if(temp<dr) dr=temp;
02948 temp=4.0+endU; if(temp<dr) dr=temp;
02949 temp=4.0-endV; if(temp<dr) dr=temp;
02950 temp=4.0+endV; if(temp<dr) dr=temp;
02951 temp=4.0-endX; if(temp<dr) dr=temp;
02952 temp=4.0+endX; if(temp<dr) dr=temp;
02953 temp=4.0-endY; if(temp<dr) dr=temp;
02954 temp=4.0+endY; if(temp<dr) dr=temp;
02955 if( dr<999.9 ) TRKendDistToEdge = dr;
02956
02957 if( TRKenddircosZ>0.0 ){
02958 if( TRKendplane<249 ){ dzm = endZ-0.110; dzp = 14.793-endZ; }
02959 if( TRKendplane>249 ){ dzm = endZ-15.986; dzp = 29.955-endZ; }
02960 }
02961 if( TRKenddircosZ<0.0 ){
02962 if( TRKendplane<249 ){ dzm = 14.793-endZ; dzp = endZ-0.110; }
02963 if( TRKendplane>249 ){ dzm = 29.955-endZ; dzp = endZ-15.986; }
02964 }
02965
02966 TRKendDistToEndBack = dzm;
02967 TRKendDistToEndForward = dzp;
02968
02969
02970 posU=endU; posV=endV; posZ=endZ;
02971 posX=0.7071*(posU-posV); posY=0.7071*(posU+posV);
02972
02973 gradu = (TRKenddircosU/TRKenddircosZ);
02974 gradv = (TRKenddircosV/TRKenddircosZ);
02975 grads = sqrt(gradu*gradu+gradv*gradv+1.0);
02976 if(TRKenddircosZ>0.0) pm=+1.0; else pm=-1.0;
02977 dirU=pm*gradu/grads; dirV=pm*gradv/grads; dirZ=pm*1.0/grads;
02978 dirX=0.7071*(dirU-dirV); dirY=0.7071*(dirU+dirV);
02979 trkflag=0; dR=0.0;
02980
02981 if(trkflag==0){
02982 if(dirX>0.0){
02983 dX = 4.0-posX; dY = dX*(dirY/dirX);
02984 if( posY+dY<4.0*(sqrt(2.0)-1.0) && posY+dY>-4.0*(sqrt(2.0)-1.0) ){
02985 trkflag=1; if(dX>0.0) dR=sqrt(dX*dX+dY*dY); if(dX<0.0) dR=-sqrt(dX*dX+dY*dY);
02986 }
02987 }
02988 if(dirX<0.0){
02989 dX = 4.0+posX; dY = -dX*(dirY/dirX);
02990 if( posY+dY<4.0*(sqrt(2.0)-1.0) && posY+dY>-4.0*(sqrt(2.0)-1.0) ){
02991 trkflag=1; if(dX>0.0) dR=sqrt(dX*dX+dY*dY); if(dX<0.0) dR=-sqrt(dX*dX+dY*dY);
02992 }
02993 }
02994 }
02995
02996 if(trkflag==0){
02997 if(dirY>0.0){
02998 dY = 4.0-posY; dX = dY*(dirX/dirY);
02999 if( posX+dX<4.0*(sqrt(2.0)-1.0) && posX+dX>-4.0*(sqrt(2.0)-1.0) ){
03000 trkflag=1; if(dY>0.0) dR=sqrt(dX*dX+dY*dY); if(dY<0.0) dR=-sqrt(dX*dX+dY*dY);
03001 }
03002 }
03003 if(dirY<0.0){
03004 dY = 4.0+posY; dX = -dY*(dirX/dirY);
03005 if( posX+dX<4.0*(sqrt(2.0)-1.0) && posX+dX>-4.0*(sqrt(2.0)-1.0) ){
03006 trkflag=1; if(dY>0.0) dR=sqrt(dX*dX+dY*dY); if(dY<0.0) dR=-sqrt(dX*dX+dY*dY);
03007 }
03008 }
03009 }
03010
03011 if(trkflag==0){
03012 if(dirU>0.0){
03013 dU = 4.0-posU; dV = dU*(dirV/dirU);
03014 if( posV+dV<4.0*(sqrt(2.0)-1.0) && posV+dV>-4.0*(sqrt(2.0)-1.0) ){
03015 trkflag=1; if(dU>0.0) dR=sqrt(dV*dV+dU*dU); if(dU<0.0) dR=-sqrt(dV*dV+dU*dU);
03016 }
03017 }
03018 if(dirU<0.0){
03019 dU = 4.0+posU; dV = -dU*(dirV/dirU);
03020 if( posV+dV<4.0*(sqrt(2.0)-1.0) && posV+dV>-4.0*(sqrt(2.0)-1.0) ){
03021 trkflag=1; if(dU>0.0) dR=sqrt(dV*dV+dU*dU); if(dU<0.0) dR=-sqrt(dV*dV+dU*dU);
03022 }
03023 }
03024 }
03025
03026 if(trkflag==0){
03027 if(dirV>0.0){
03028 dV = 4.0-posV; dU = dV*(dirU/dirV);
03029 if( posU+dU<4.0*(sqrt(2.0)-1.0) && posU+dU>-4.0*(sqrt(2.0)-1.0) ){
03030 trkflag=1; if(dV>0.0) dR=sqrt(dV*dV+dU*dU); if(dV<0.0) dR=-sqrt(dV*dV+dU*dU);
03031 }
03032 }
03033 if(dirV<0.0){
03034 dV = 4.0+posV; dU = -dV*(dirU/dirV);
03035 if( posU+dU<4.0*(sqrt(2.0)-1.0) && posU+dU>-4.0*(sqrt(2.0)-1.0) ){
03036 trkflag=1; if(dV>0.0) dR=sqrt(dV*dV+dU*dU); if(dV<0.0) dR=-sqrt(dV*dV+dU*dU);
03037 }
03038 }
03039 }
03040
03041 traceZ=0.0; trace=0.0;
03042 if(sqrt(dirX*dirX+dirY*dirY)>0.0){
03043 traceZ=(pm*dR*dirZ)/sqrt(dirX*dirX+dirY*dirY);
03044 }
03045 trace=(traceZ)/(pm*dirZ);
03046
03047 TRKendtrace = trace;
03048 TRKendtraceZ = traceZ;
03049 }
03050
03051
03052 if( vldc.GetDetector()==Detector::kNear ){
03053 TRKvtxDistToEdge = track->fidvtx.dr;
03054 TRKvtxDistToEndBack = track->fidvtx.dz;
03055 TRKvtxDistToEndForward = track->fidvtx.dz;
03056 TRKvtxtrace = track->fidvtx.trace;
03057 TRKvtxtraceZ = track->fidvtx.tracez;
03058 TRKendDistToEdge = track->fidend.dr;
03059 TRKendDistToEndBack = track->fidend.dz;
03060 TRKendDistToEndForward = track->fidend.dz;
03061 TRKendtrace = track->fidend.trace;
03062 TRKendtraceZ = track->fidend.tracez;
03063 }
03064
03065
03066 Int_t bstrp,estrp;
03067 Double_t trkph,evtph;
03068
03069 Double_t sumtrkph = 0.0;
03070 Double_t sumevtph = 0.0;
03071 Int_t ntrackplanes = 0;
03072 Int_t ntracklikeplanes = 0;
03073 Int_t nshowerlikeplanes = 0;
03074
03075 for( Int_t n=0; n<500; n++ ){
03076
03077 Int_t dpln = 1;
03078 if( vldc.GetDetector()==Detector::kNear ){
03079 if( n>120 ) dpln = 5;
03080 }
03081
03082 if( 1+fStripList[n].GetLast()>0 ){
03083 for( Int_t ns=0; ns<1+fStripList[n].GetLast(); ns++ ){
03084 NtpSRStrip* strip = (NtpSRStrip*)(fStripList[n].At(ns));
03085 sumevtph += strip->ph0.sigcor+strip->ph1.sigcor;
03086 }
03087 }
03088
03089 if( 1+fTrackStripList[n].GetLast()>0 ){
03090 for( Int_t nt=0; nt<1+fTrackStripList[n].GetLast(); nt++ ){
03091 NtpSRStrip* strip = (NtpSRStrip*)(fTrackStripList[n].At(nt));
03092 sumtrkph += strip->ph0.sigcor+strip->ph1.sigcor;
03093 }
03094 }
03095
03096 if( 1+fStripList[n].GetLast()>0
03097 && 1+fTrackStripList[n].GetLast()>0 ){
03098
03099 bstrp=-1; estrp=-1;
03100 for( Int_t nt=0; nt<1+fTrackStripList[n].GetLast(); nt++ ){
03101 NtpSRStrip* strip = (NtpSRStrip*)(fTrackStripList[n].At(nt));
03102 if( bstrp<0 || strip->strip<bstrp ) bstrp=strip->strip;
03103 if( estrp<0 || strip->strip>estrp ) estrp=strip->strip;
03104 }
03105
03106 if(bstrp>=0 && estrp>=0 && estrp-bstrp>=0){
03107 trkph=0.0; evtph=0.0;
03108 for( Int_t ns=0; ns<1+fStripList[n].GetLast(); ns++ ){
03109 NtpSRStrip* strip = (NtpSRStrip*)(fStripList[n].At(ns));
03110 evtph += strip->ph0.sigcor+strip->ph1.sigcor;
03111 if( strip->strip>=bstrp-1 && strip->strip<=estrp+1 ){
03112 trkph += strip->ph0.sigcor+strip->ph1.sigcor;
03113 }
03114 }
03115 if( evtph>0.0 ){
03116 if( trkph>0.0 ) ntrackplanes += 1;
03117 if( trkph/evtph>0.8 ) ntracklikeplanes += 1;
03118 else nshowerlikeplanes += 1;
03119 }
03120 }
03121 }
03122 }
03123
03124 TRKtotph = sumtrkph;
03125 TRKplanes = ntrackplanes;
03126 TRKtracklikeplanes = ntracklikeplanes;
03127 TRKshowerlikeplanes = nshowerlikeplanes;
03128
03129 }
03130 }
03131
03132
03133
03134 Nshowers = ntpEvent->nshower;
03135
03136 if( ntpEvent->nshower>0 ){
03137 MSG("MadAb",Msg::kDebug) << " found reconstructed shower(s) *** " << endl;
03138
03139 Double_t eshw = 0.0;
03140 Int_t shwcontained = 0;
03141
03142 Int_t* shw = ntpEvent->shw;
03143 for( Int_t n=0; n<ntpEvent->nshower; n++ ){
03144 Int_t ptr=shw[n];
03145
03146 if( ptr>=0 ){
03147 NtpSRShower* shower = (NtpSRShower*)(showerArray->At(ptr));
03148
03149
03150 Bool_t vtxshw = 0;
03151 if( ( shower->vtx.z-vertex.z<0.5 )
03152 || ( shower->vtx.z-vertex.z>=0.5
03153 && shower->shwph.wtCCgev>1.0 ) ){
03154 vtxshw = 1;
03155 }
03156
03157
03158 if( vtxshw ){
03159 eshw += shower->shwph.linCCgev;
03160 }
03161
03162
03163 if( shwcontained==0 && shower->contained ){
03164 shwcontained = shower->contained;
03165 }
03166
03167
03168 Int_t* shwstp = shower->stp;
03169 for( Int_t k=0; k<shower->nstrip; k++ ){
03170 Int_t sptr = shwstp[k];
03171 if( sptr>=0 ){
03172 NtpSRStrip* strip = (NtpSRStrip*)(stripArray->At(sptr));
03173 Int_t plane = strip->plane;
03174 if( plane>=0 && plane<500 ){
03175 if( vtxshw ) fShwShowerStripList[plane].Add(strip);
03176 else fTrkShowerStripList[plane].Add(strip);
03177 }
03178 if( sptr<last ){
03179 shwflag[sptr] = 1;
03180 }
03181 }
03182 }
03183
03184 }
03185 }
03186
03187 SHWenergy = eshw;
03188 SHWcontained = shwcontained;
03189 }
03190
03191 Int_t sumplanes = 0;
03192 Int_t sumtrkplanes = 0;
03193 Int_t sumshwplanes = 0;
03194 Double_t sumtrktrkph = 0.0;
03195 Double_t sumtrkshwph = 0.0;
03196
03197 for( Int_t n=0; n<last; n++ ){
03198 Int_t ptr=n;
03199 NtpSRStrip* strip = (NtpSRStrip*)(stripArray->At(ptr));
03200 Int_t plane = strip->plane;
03201 if( plane>=0 && plane<500 ){
03202
03203 Int_t dpln = 1;
03204 if( vldc.GetDetector()==Detector::kNear ){
03205 if( plane>120 ) dpln = 5;
03206 }
03207
03208 if( trkflag[n]==1 && shwflag[n]==0 ){
03209 sumtrktrkph += (strip->ph0.sigcor+strip->ph1.sigcor);
03210 if( flag[plane]==0 ){
03211 flag[plane]=1; sumtrkplanes += 1; sumplanes += dpln;
03212 }
03213 }
03214
03215 if( trkflag[n]==1 && shwflag[n]==1 ){
03216 sumtrkshwph += (strip->ph0.sigcor+strip->ph1.sigcor);
03217 if( flag[plane]==0 ){
03218 flag[plane]=1; sumshwplanes += 1; sumplanes += dpln;
03219 }
03220 }
03221 }
03222 }
03223
03224
03225
03226
03227
03228 if( vldc.GetDetector()==Detector::kFar ){
03229 for( Int_t n=0; n<last; n++ ){
03230 Int_t ptr=n;
03231 NtpSRStrip* strip = (NtpSRStrip*)(stripArray->At(ptr));
03232
03233 if( strip ){
03234 Int_t plane = strip->plane;
03235 if( plane>=0 && plane<500 ){
03236
03237 if( trkflag[n]==0 && shwflag[n]==0 ){
03238 Bool_t ontrack = 0;
03239 for( Int_t s=0; s<1+fTrackStripList[plane].GetLast(); s++ ){
03240 NtpSRStrip* trackstrip = (NtpSRStrip*)(fTrackStripList[plane].At(s));
03241 if( abs(trackstrip->strip-strip->strip)<=1 ) ontrack=1;
03242 }
03243 if( ontrack==0 ){
03244 if( strip->ph0.raw+strip->ph1.raw>200.0 ){
03245
03246 }
03247 }
03248 }
03249 }
03250 }
03251
03252 }
03253 }
03254
03255
03256 if( vldc.GetDetector()==Detector::kNear ){
03257
03258 }
03259
03260 TRKtrkph = sumtrktrkph;
03261 TRKshwph = sumtrkshwph;
03262 TRKtrackplanes = sumtrkplanes;
03263 TRKshowerplanes = sumshwplanes;
03264 TRKspanplanes = sumplanes;
03265
03266
03267
03268
03269
03270 if( vldc.GetDetector()==Detector::kFar ){
03271 RECOvtxR = sqrt(RECOvtxX*RECOvtxX+RECOvtxY*RECOvtxY);
03272 }
03273 if( vldc.GetDetector()==Detector::kNear ){
03274 RECOvtxR = sqrt( (RECOvtxX-1.4828)*(RECOvtxX-1.4828)
03275 + (RECOvtxY-0.2384)*(RECOvtxY-0.2384) );
03276 }
03277
03278 RECOcontained = ntpEvent->contained;
03279
03280
03281 if( vldc.GetDetector()==Detector::kFar ){
03282
03283 Bool_t vtxcontained = 0;
03284 Bool_t endcontained = 0;
03285
03286 Double_t evtxU = RECOvtxU;
03287 Double_t evtxV = RECOvtxV;
03288 Double_t evtxX = RECOvtxX;
03289 Double_t evtxY = RECOvtxY;
03290 Double_t evtxZ = RECOvtxZ;
03291
03292 dr = 999.9;
03293 temp=4.0-evtxU; if(temp<dr) dr=temp;
03294 temp=4.0+evtxU; if(temp<dr) dr=temp;
03295 temp=4.0-evtxV; if(temp<dr) dr=temp;
03296 temp=4.0+evtxV; if(temp<dr) dr=temp;
03297 temp=4.0-evtxX; if(temp<dr) dr=temp;
03298 temp=4.0+evtxX; if(temp<dr) dr=temp;
03299 temp=4.0-evtxY; if(temp<dr) dr=temp;
03300 temp=4.0+evtxY; if(temp<dr) dr=temp;
03301 if(dr<999.9) RECOvtxDistToEdge = dr;
03302
03303
03304 Double_t evtxdircosZ = 1.0;
03305
03306 if( evtxdircosZ>=0.0 ){
03307 if( RECOvtxplane<249 ){ dzm = evtxZ-0.110; dzp = 14.793-evtxZ; }
03308 if( RECOvtxplane>249 ){ dzm = evtxZ-15.986; dzp = 29.955-evtxZ; }
03309 }
03310 if( evtxdircosZ<0.0 ){
03311 if( RECOvtxplane<249 ){ dzm = 14.793-evtxZ; dzp = evtxZ-0.110; }
03312 if( RECOvtxplane>249 ){ dzm = 29.955-evtxZ; dzp = evtxZ-15.986; }
03313 }
03314
03315 RECOvtxDistToEndBack = dzm;
03316 RECOvtxDistToEndForward = dzp;
03317
03318
03319 vtxcontained = 0;
03320
03321
03322
03323
03324
03325
03326
03327
03328
03329
03330
03331 if (Ntracks>0 && TRKvtxR<sqrt(14.)
03332 && ((TRKvtxplane>4 && TRKvtxplane<241) ||
03333 (TRKvtxplane>253 && TRKvtxplane<466))) {
03334 vtxcontained = 1;
03335 }
03336
03337
03338
03339
03340 endcontained = 1;
03341 if( TRKvtxplane>=0 ){
03342 endcontained = 0;
03343 if( TRKendDistToEdge>0.2 && TRKendR>0.4
03344 && TRKendDistToEndBack>0.25 && TRKendDistToEndForward>0.25 ){
03345 endcontained = 1;
03346 }
03347 }
03348
03349 RECOvtxcontained = vtxcontained;
03350 RECOendcontained = endcontained;
03351 }
03352
03353
03354 if( vldc.GetDetector()==Detector::kNear ){
03355 RECOvtxcontained = 0;
03356
03357
03358
03359
03360
03361
03362
03363
03364
03365 if(Ntracks>0 && TRKvtxR<1.0 && TRKvtxZ>1.0 && TRKvtxZ<5.0) {
03366 RECOvtxcontained = 1;
03367 }
03368
03369
03370 RECOendcontained = RECOcontained;
03371 RECOvtxDistToEdge = TRKvtxDistToEdge;
03372 RECOvtxDistToEndBack = TRKvtxDistToEndBack;
03373 RECOvtxDistToEndForward = TRKvtxDistToEndForward;
03374 }
03375
03376
03377
03378 RECOdircosU = TRKvtxdircosU;
03379 RECOdircosV = TRKvtxdircosV;
03380 RECOdircosX = TRKvtxdircosX;
03381 RECOdircosY = TRKvtxdircosY;
03382 RECOdircosZ = TRKvtxdircosZ;
03383 RECOdircosneu = 0.99834*RECOdircosZ+0.05759*RECOdircosY;
03384
03385
03386
03387
03388 if( RECOendcontained || TRKmomentumCurve==0){
03389 RECOemu = TRKmomentumRange;
03390 }
03391 else{
03392 RECOemu = TRKmomentumCurve;
03393 }
03394 RECOehad = SHWenergy;
03395
03396 RECOenu = RECOemu+RECOehad;
03397
03398 if( RECOenu>0.0 && RECOemu>0.0 ){
03399 Double_t Mp = 0.5*(0.93827+0.93957);
03400 RECOy = RECOehad/RECOenu;
03401 RECOq2 = 2.0*RECOenu*RECOemu*(1.0-RECOdircosneu);
03402 RECOw2 = Mp*Mp-RECOq2+2*Mp*RECOehad;
03403 if( RECOehad>0.0 ) RECOx = RECOq2/(2.0*Mp*RECOehad);
03404 else RECOx = 1.0;
03405
03406 RECOQEenu = (Mp*RECOemu)/(Mp-RECOemu*(1.0-RECOdircosneu));
03407 RECOQEemu = RECOemu;
03408 if( RECOQEenu>0.0 ){
03409 RECOQEw2 = Mp*Mp;
03410 RECOQEq2 = 2.0*RECOQEenu*RECOQEemu*(1.0-RECOdircosneu);
03411 RECOQEehad = RECOQEq2/2*Mp;
03412 RECOQEy = RECOQEehad/RECOQEenu;
03413 RECOQEx = 1.0;
03414 }
03415 }
03416
03417 delete [] shwflag;
03418 delete [] trkflag;
03419 delete [] flag;
03420 }
03421
03422 MSG("MadAb",Msg::kVerbose) << " *** MadAbPid::MakeRecoVariables(...) FINISHED *** " << endl;
03423
03424 return;
03425 }
