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