00001 #include <iostream>
00002
00003 #include <cstdlib>
00004
00005 #include "ParticleID/ParticleID.h"
00006
00007 #ifndef BEAN
00008 #include "EvtRecEvent/EvtRecTrack.h"
00009 #include "GaudiKernel/Bootstrap.h"
00010 #include "GaudiKernel/ISvcLocator.h"
00011 #include "GaudiKernel/ISvcLocator.h"
00012 #include "GaudiKernel/IDataProviderSvc.h"
00013 #include "GaudiKernel/SmartDataPtr.h"
00014 #include "EventModel/EventHeader.h"
00015 #endif
00016
00017
00018
00019
00020 ParticleID * ParticleID::m_pointer = 0;
00021
00022 ParticleID * ParticleID::instance() {
00023 if(!m_pointer) m_pointer = new ParticleID();
00024 return m_pointer;
00025 }
00026
00027 void ParticleID::init() {
00028
00029 if(IsDedxInfoUsed()) {
00030 if(!m_dedxpid) m_dedxpid = DedxPID::instance();
00031 m_dedxpid->init();
00032 }
00033
00034 if(IsTofInfoUsed()|IsTof1InfoUsed()|IsTof2InfoUsed()) {
00035 if(!m_tofpid) m_tofpid = TofPID::instance();
00036 m_tofpid->init();
00037 }
00038
00039 if(IsTofCorrInfoUsed()) {
00040 if(!m_tofcorrpid) m_tofcorrpid = TofCorrPID::instance();
00041
00042 }
00043
00044 if(IsEmcInfoUsed()) {
00045 if(!m_emcpid) m_emcpid = EmcPID::instance();
00046 m_emcpid->init();
00047 }
00048
00049 if(IsMucInfoUsed()) {
00050 if(!m_mucpid) m_mucpid = MucPID::instance();
00051 m_mucpid->init();
00052 }
00053
00054
00055 m_pidsys = 0;
00056 m_pidcase = 0;
00057 m_method = 0;
00058 m_TotalLikelihood =0;
00059 m_discard =1;
00060
00061 m_ndof = 0;
00062 m_nhitcut=5;
00063
00064 setChiMinCut( 4.0 );
00065 setChiMaxCut( 6.0 );
00066 for(int i = 0; i < 4; i++) {
00067 m_chisq[i] = 9999.;
00068 m_prob[i] = -1.0;
00069 }
00070 }
00071
00072
00073 ParticleID::ParticleID() : ParticleIDBase() {
00074 m_dedxpid = 0;
00075 m_tofpid = 0;
00076 m_tofepid = 0;
00077 m_tofqpid = 0;
00078 m_tofcpid = 0;
00079 m_tofcorrpid = 0;
00080 m_emcpid = 0;
00081 m_mucpid = 0;
00082
00083
00084
00085 }
00086
00087
00088 ParticleID::~ParticleID() {
00089
00090
00091
00092
00093
00094
00095 }
00096
00097 void ParticleID::calculate() {
00098 #ifdef BEAN
00099 cout << " please use ParticleID::calculate(run) ! " << endl;
00100 exit(1);
00101 }
00102
00103
00104 void ParticleID::calculate(int run)
00105 {
00106 #endif
00107 int nhitcutpid=getNhitCut();
00108
00109 #ifndef BEAN
00110 IDataProviderSvc* m_eventSvc;
00111 Gaudi::svcLocator()->service("EventDataSvc", m_eventSvc, true);
00112
00113 SmartDataPtr<Event::EventHeader> eventHeaderpid(m_eventSvc,"/Event/EventHeader");
00114 int runpid=eventHeaderpid->runNumber();
00115 int eventpid=eventHeaderpid->eventNumber();
00116
00117
00118 #else
00119 int runpid=run;
00120 #endif
00121
00122 EvtRecTrack* recTrk = PidTrk();
00123
00124
00125
00126 if(m_pidsys == 0) {
00127 m_pidsys = useDedx() | useTof() | useTofE() | useEmc() | useMuc() | useTofQ() | useTofC() | useTofCorr();
00128 }
00129
00130
00131 if(m_pidcase == 0 ) {
00132 m_pidcase = all();
00133 }
00134
00135 if(IsDedxInfoUsed()) {
00136 if(!m_dedxpid) m_dedxpid = DedxPID::instance();
00137 m_dedxpid->init();
00138 m_dedxpid->setRunNo(runpid);
00139 m_dedxpid->setNhitCutDx(nhitcutpid);
00140 m_dedxpid->setRecTrack(recTrk);
00141 m_dedxpid->setChiMinCut(chiMinCut());
00142 m_dedxpid->setPdfMinSigmaCut(pdfMinSigmaCut());
00143 m_dedxpid->calculate();
00144 }
00145
00146
00147 if(IsTofInfoUsed()|IsTof1InfoUsed()|IsTof2InfoUsed()|IsTofCInfoUsed())
00148 {
00149 if(IsTofCInfoUsed())
00150 {
00151 if(!m_tofcpid) m_tofcpid = TofCPID::instance();
00152 m_tofcpid->init();
00153 m_tofcpid->setRunNo(runpid);
00154 m_tofcpid->setRecTrack(recTrk);
00155 m_tofcpid->setChiMinCut(chiMinCut());
00156 m_tofcpid->setPdfMinSigmaCut(pdfMinSigmaCut());
00157 m_tofcpid->calculate();
00158 }
00159 else
00160 {
00161 if(!m_tofpid) m_tofpid = TofPID::instance();
00162 m_tofpid->init();
00163 m_tofpid->setRecTrack(recTrk);
00164 m_tofpid->setChiMinCut(chiMinCut());
00165 m_tofpid->setPdfMinSigmaCut(pdfMinSigmaCut());
00166 m_tofpid->calculate();
00167 }
00168
00169 }
00170
00171
00172
00173 if(IsTofCorrInfoUsed()) {
00174 if(!m_tofcorrpid) m_tofcorrpid = TofCorrPID::instance();
00175 m_tofcorrpid->setRunNo(runpid);
00176 m_tofcorrpid->init();
00177 m_tofcorrpid->setRecTrack(recTrk);
00178 m_tofcorrpid->setChiMinCut(chiMinCut());
00179 m_tofcorrpid->setChiMaxCut(chiMaxCut());
00180 m_tofcorrpid->setPdfMinSigmaCut(pdfMinSigmaCut());
00181 m_tofcorrpid->calculate();
00182 }
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207 if(IsTofQInfoUsed()) {
00208 if(!m_tofqpid) m_tofqpid = TofQPID::instance();
00209 m_tofqpid->init();
00210 m_tofqpid->setRecTrack(recTrk);
00211 m_tofqpid->setChiMinCut(chiMinCut());
00212 m_tofqpid->calculate();
00213 }
00214
00215
00216 if(IsTofEInfoUsed()&&(!IsTofCorrInfoUsed())) {
00217 if(!m_tofepid) m_tofepid = TofEPID::instance();
00218 m_tofepid->init();
00219 m_tofepid->setRecTrack(recTrk);
00220 m_tofepid->setChiMinCut(chiMinCut());
00221 m_tofepid->setPdfMinSigmaCut(pdfMinSigmaCut());
00222 m_tofepid->calculate();
00223 }
00224
00225 if(IsEmcInfoUsed()) {
00226 if(!m_emcpid) m_emcpid = EmcPID::instance();
00227 m_emcpid->init();
00228 m_emcpid->setRecTrack(recTrk);
00229 m_emcpid->setChiMinCut(chiMinCut());
00230 m_emcpid->calculate();
00231 }
00232
00233
00234 if(IsMucInfoUsed()) {
00235 if(!m_mucpid) m_mucpid = MucPID::instance();
00236 m_mucpid->init();
00237 m_mucpid->setRecTrack(recTrk);
00238 m_mucpid->setChiMinCut(chiMinCut());
00239 m_mucpid->calculate();
00240 }
00241
00242 if((m_method & methodProbability()) == methodProbability())
00243 if(particleIDCalculation() < 0) m_ndof = 0;
00244
00245
00246
00247 if((m_method & methodLikelihood()) == methodLikelihood())
00248 if(LikelihoodCalculation() < 0) m_discard =0;
00249
00250 if((m_method & methodNeuronNetwork()) == methodNeuronNetwork())
00251
00252 if(LikelihoodCalculation() < 0) m_discard =0;
00253
00254 }
00255
00256 int ParticleID ::particleIDCalculation() {
00257 int irc = -1;
00258 bool valid = IsDedxInfoValid() || IsTofInfoValid()||IsTofEInfoValid()
00259 || IsTofQInfoValid() || IsEmcInfoValid() || IsMucInfoValid()
00260 || IsTofCInfoValid() || IsTofCorrInfoValid();
00261
00262 if(!valid) return irc;
00263
00264 double chisq[5];
00265 bool pidcase[5];
00266 for(int i = 0; i < 5; i++) {
00267 chisq[i] = 0;
00268 pidcase[i] = false;
00269 }
00270
00271 if((m_pidcase & onlyElectron()) == onlyElectron()) pidcase[0] = true;
00272 if((m_pidcase & onlyMuon()) == onlyMuon()) pidcase[1] = true;
00273 if((m_pidcase & onlyPion()) == onlyPion()) pidcase[2] = true;
00274 if((m_pidcase & onlyKaon()) == onlyKaon()) pidcase[3] = true;
00275 if((m_pidcase & onlyProton()) == onlyProton()) pidcase[4] = true;
00276
00277
00278
00279
00280 if(IsDedxInfoUsed()) {
00281 if(IsDedxInfoValid()) {
00282 bool okpid = false;
00283 for(int i = 0; i < 5; i++) {
00284 if(pidcase[i] && (fabs(chiDedx(i)) < m_dedxpid->chiMinCut()))
00285 if(!okpid) okpid = true;
00286 }
00287 if(okpid) {
00288 m_ndof++;
00289 for(int i = 0; i < 5; i++) chisq[i] += chiDedx(i)*chiDedx(i);
00290
00291
00292 }
00293 }
00294 }
00295
00296
00297
00298
00299 if(IsTofInfoUsed()|IsTof1InfoUsed()|IsTof2InfoUsed() | IsTofCInfoUsed())
00300 { if(IsTofCInfoUsed())
00301 {
00302 if(IsTofCInfoValid()) {
00303 bool okpid = false;
00304 for(int i = 0; i < 5; i++) {
00305 if(pidcase[i] && (fabs(chiTofC(i)) < m_tofcpid->chiMinCut()))
00306 if(!okpid) okpid = true;
00307 }
00308 if(okpid) {
00309 m_ndof++;
00310 for(int i = 0; i < 5; i++) chisq[i] += chiTofC(i)*chiTofC(i);
00311 }
00312 }
00313 }
00314 else {
00315 if(IsTofInfoValid()) {
00316 bool okpid = false;
00317 for(int i = 0; i < 5; i++) {
00318 if(pidcase[i] && (fabs(chiTof(i)) < m_tofpid->chiMinCut()))
00319 if(!okpid) okpid = true;
00320 }
00321 if(okpid) {
00322 m_ndof++;
00323 for(int i = 0; i < 5; i++) chisq[i] += chiTof(i)*chiTof(i);
00324 }
00325 }
00326
00327
00328
00329
00330
00331
00332 if(IsTofEInfoUsed()) {
00333 if(IsTofEInfoValid()) {
00334 bool okpid = false;
00335 for(int i = 0; i < 5; i++) {
00336 if(pidcase[i] && (fabs(chiTofE(i)) < m_tofepid->chiMinCut()))
00337 if(!okpid) okpid = true;
00338 }
00339 if(okpid) {
00340 m_ndof++;
00341 for(int i = 0; i < 5; i++) chisq[i] += chiTofE(i)*chiTofE(i);
00342 }
00343 }
00344 }
00345
00346 }
00347 }
00348
00349
00350 if(IsTofCorrInfoUsed()) {
00351 if(IsTofCorrInfoValid()) {
00352 bool okpid = false;
00353 for(int i = 0; i < 5; i++) {
00354 if(pidcase[i] && ( chiTofCorr(i) < m_tofcorrpid->chiMaxCut() ) && ( chiTofCorr(i) > ( 0.0 - m_tofcorrpid->chiMinCut() ) ) )
00355
00356 if(!okpid) okpid = true;
00357 }
00358 if(okpid) {
00359 m_ndof++;
00360 for(int i = 0; i < 5; i++) chisq[i] += chiTofCorr(i)*chiTofCorr(i);
00361 }
00362 }
00363 }
00364
00365
00366
00367
00368
00369 if(IsTofQInfoUsed()) {
00370 if(IsTofQInfoValid()) {
00371 bool okpid = false;
00372 for(int i = 0; i < 5; i++) {
00373 if(pidcase[i] && (fabs(chiTofQ(i)) < m_tofqpid->chiMinCut()))
00374 if(!okpid) okpid = true;
00375 }
00376 if(okpid) {
00377 m_ndof++;
00378 for(int i = 0; i < 5; i++) chisq[i] += chiTofQ(i)*chiTofQ(i);
00379 }
00380 }
00381 }
00382
00383
00384 if(IsMucInfoUsed()) {
00385 if(IsMucInfoValid()) {
00386 bool okpid = false;
00387 for(int i = 0; i < 5; i++) {
00388 if(pidcase[i] && (fabs(chiMuc(i)) < m_mucpid->chiMinCut()))
00389 if(!okpid) okpid = true;
00390 }
00391 if(okpid) {
00392 m_ndof++;
00393 for(int i = 0; i < 5; i++) chisq[i] += chiMuc(i)*chiMuc(i);
00394 }
00395 }
00396 }
00397
00398
00399
00400 if(IsEmcInfoUsed()) {
00401 if(IsEmcInfoValid()) {
00402 bool okpid = false;
00403 for(int i = 0; i < 5; i++) {
00404 if(pidcase[i] && (fabs(chiEmc(i)) < m_emcpid->chiMinCut()))
00405 if(!okpid) okpid = true;
00406 }
00407 if(okpid) {
00408 m_ndof++;
00409 for(int i = 0; i < 5; i++) chisq[i] += chiEmc(i)*chiEmc(i);
00410 }
00411 }
00412 }
00413
00414
00415 if(m_ndof <= 0) return irc;
00416
00417
00418 for(int i = 0; i < 5; i++) {
00419 m_chisq[i] = chisq[i];
00420 m_prob[i] = probCalculate(chisq[i], m_ndof);
00421 }
00422
00423
00424 irc = 0;
00425 return irc;
00426 }
00427
00428
00429
00430
00431
00432 int ParticleID::LikelihoodCalculation() {
00433 int irc = -1;
00434
00435 bool valid = IsDedxInfoValid() || IsTofInfoValid() || IsTofEInfoValid() || IsTofQInfoValid() || IsTofCInfoValid() || IsTofCorrInfoValid() || IsEmcInfoValid()||IsMucInfoValid();
00436 if(!valid) return irc;
00437 double pdf[5];
00438 bool pidcase[5];
00439 for(int i = 0; i < 5; i++) {
00440 pdf[i] = 1;
00441 pidcase[i] = false;
00442 }
00443
00444 if((m_pidcase & onlyElectron()) == onlyElectron()) pidcase[0] = true;
00445 if((m_pidcase & onlyMuon()) == onlyMuon()) pidcase[1] = true;
00446 if((m_pidcase & onlyPion()) == onlyPion()) pidcase[2] = true;
00447 if((m_pidcase & onlyKaon()) == onlyKaon()) pidcase[3] = true;
00448 if((m_pidcase & onlyProton()) == onlyProton()) pidcase[4] = true;
00449
00450 for(int i = 0; i < 5; i++) {
00451 if(pidcase[i]==0)
00452 pdf[i]=0;
00453 }
00454
00455
00456
00457
00458 if(IsDedxInfoUsed()) {
00459 if(IsDedxInfoValid()) {
00460 bool okpid = false;
00461 for(int i = 0; i < 5; i++) {
00462 if(pidcase[i] && pdfCalculate(chiDedx(i),1) > pdfCalculate(m_dedxpid->pdfMinSigmaCut(),1.5))
00463 if(!okpid) okpid = true;
00464 }
00465 if(okpid) {
00466 m_ndof++;
00467 for(int i = 0; i < 5; i++) {
00468 pdf[i] *= pdfDedx(i);
00469 }
00470 }
00471 }
00472 }
00473
00474
00475
00476
00477
00478 if(IsTofInfoUsed()|IsTof1InfoUsed()|IsTof2InfoUsed()|IsTofCInfoUsed())
00479 { if(IsTofCInfoUsed())
00480 {
00481
00482 if(IsTofCInfoValid()) {
00483 bool okpid = false;
00484 for(int i = 0; i < 5; i++) {
00485 if(pidcase[i] && pdfCalculate(chiTof(i),1) > pdfCalculate(m_tofcpid->pdfMinSigmaCut(),1.5))
00486 if(!okpid) okpid = true;
00487 }
00488 if(okpid) {
00489 m_ndof++;
00490 for(int i = 0; i < 5; i++) {
00491 pdf[i] *= pdfTofC(i);
00492 }
00493 }
00494 }
00495 }
00496
00497 else {
00498 if(IsTofInfoValid()) {
00499 bool okpid = false;
00500 for(int i = 0; i < 5; i++) {
00501 if(pidcase[i] && pdfCalculate(chiTof(i),1) > pdfCalculate(m_tofpid->pdfMinSigmaCut(),1.5))
00502 if(!okpid) okpid = true;
00503 }
00504 if(okpid) {
00505 m_ndof++;
00506 for(int i = 0; i < 5; i++) {
00507 pdf[i] *= pdfTof(i);
00508 }
00509 }
00510 }
00511
00512
00513
00514
00515
00516
00517
00518 if(IsTofEInfoUsed()) {
00519 if(IsTofEInfoValid()) {
00520 bool okpid = false;
00521 for(int i = 0; i < 5; i++) {
00522 if(pidcase[i]&& pdfCalculate(chiTofE(i),1) > pdfCalculate(m_tofepid->pdfMinSigmaCut(),1.5))
00523 if(!okpid) okpid = true;
00524 }
00525 if(okpid) {
00526
00527
00528 }
00529 }
00530 }
00531 }
00532
00533 }
00534
00535
00536 if(IsTofCorrInfoUsed()) {
00537 if(IsTofCorrInfoValid()) {
00538 bool okpid = false;
00539 for(int i = 0; i < 5; i++) {
00540 if(pidcase[i] && pdfCalculate(chiTofCorr(i),1) > pdfCalculate(m_tofcorrpid->pdfMinSigmaCut(),1.5))
00541 if(!okpid) okpid = true;
00542 }
00543 if(okpid) {
00544 m_ndof++;
00545 for(int i = 0; i < 5; i++) {
00546 pdf[i] *= pdfTofCorr(i);
00547 }
00548 }
00549 }
00550 }
00551
00552
00553
00554
00555
00556
00557 if(IsTofQInfoValid()) {
00558 bool okpid = false;
00559 for(int i = 0; i < 5; i++) {
00560 if(pidcase[i])
00561 if(!okpid) okpid = true;
00562 }
00563 if(okpid) {
00564
00565 for(int i = 0; i < 5; i++) pdf[i] *= pdfTofQ(i);
00566 }
00567 }
00568
00569
00570
00571
00572 if(IsEmcInfoUsed()) {
00573 if(IsEmcInfoValid()) {
00574 bool okpid = false;
00575 for(int i = 0; i < 5; i++) {
00576 if(pidcase[i]&&pdfEmc(i)>0)
00577 if(!okpid) okpid = true;
00578 }
00579 if(okpid) {
00580 m_ndof++;
00581 for(int i = 0; i < 5; i++) {
00582 pdf[i] *= pdfEmc(i);
00583 }
00584 }
00585 }
00586 }
00587 if(IsMucInfoUsed()) {
00588 if(IsMucInfoValid()) {
00589 bool okpid = false;
00590 for(int i = 0; i < 5; i++) {
00591 if(pidcase[i]&&pdfMuc(i)>0)
00592 if(!okpid) okpid = true;
00593 }
00594 if(okpid) {
00595 m_ndof++;
00596 for(int i = 0; i < 5; i++) {
00597 pdf[i] *= pdfMuc(i);
00598 }
00599 }
00600 }
00601 }
00602
00603
00604
00605 if(m_ndof <= 0) return irc;
00606 for(int i = 0; i < 5; i++) {
00607 m_pdf[i] = pdf[i];
00608 m_TotalLikelihood += pdf[i];
00609 }
00610 for(int i = 0; i < 5; i++) {
00611 m_likelihoodfraction[i] = pdf[i]/m_TotalLikelihood;
00612 }
00613
00614
00615 irc = 0;
00616 return irc;
00617 }
00618