#include <MdaDiscrimAna.h>
Public Types | |
| typedef std::deque< std::deque< Int_t > > | DeqDeqInt_t |
| typedef DeqDeqInt_t::iterator | IterDeqDeqInt_t |
| typedef DeqDeqInt_t::size_type | SizeDeqDeqInt_t |
| typedef std::deque< Double_t > | DeqDouble_t |
| typedef DeqDouble_t::iterator | IterDeqDouble_t |
| typedef std::deque< TMatrixD > | DeqTMatrixD |
| typedef DeqTMatrixD::iterator | IterDeqTMatrixD |
| typedef std::deque< TVectorD > | DeqTVectorD |
| typedef DeqTVectorD::iterator | IterDeqTVectorD |
Public Member Functions | |
| MdaDiscrimAna (NueRecord &nr, MdaDiscrim &md) | |
| ~MdaDiscrimAna () | |
| void | Analyze () |
| void | Reset () |
| Bool_t | NeedsSpecialAttention (TString name, NueRecord *nr, Double_t &value) |
| void | SetMdaParams (Double_t threshC, std::string sasFN) |
| bool | isFillDone () const |
Public Attributes | |
| DeqTVectorD | varPVector |
| DeqTVectorD | discrimDiff |
| DeqDouble_t | probClass |
Static Public Attributes | |
| bool | isFilled |
Private Member Functions | |
| MdaDiscrimAna (const MdaDiscrimAna &rhs) | |
| MdaDiscrimAna & | operator= (const MdaDiscrimAna &rhs) |
| Bool_t | FillCalibArrays () |
| void | FillPVector () |
| void | PIDCalc () |
| void | MdaClassify () |
Private Attributes | |
| NueRecord & | nueRec |
| MdaDiscrim & | fMdaDiscrim |
| Double_t | threshCut |
| std::string | sasFileName |
Static Private Attributes | |
| Int_t | nClass |
| std::deque< TMatrixD > | quadCoeff |
| std::deque< TVectorD > | linCoeff |
| std::deque< TVectorD > | constCoeff |
| std::deque< TVectorD > | meanVec |
| std::deque< std::deque< std::string > > | typeClassVec |
| std::deque< std::string > | typeCoeffVec |
| std::deque< std::string > | varNameVec |
|
|
Definition at line 34 of file MdaDiscrimAna.h. |
|
|
Definition at line 38 of file MdaDiscrimAna.h. Referenced by PIDCalc(). |
|
|
Definition at line 41 of file MdaDiscrimAna.h. |
|
|
Definition at line 44 of file MdaDiscrimAna.h. |
|
|
Definition at line 35 of file MdaDiscrimAna.h. |
|
|
Definition at line 39 of file MdaDiscrimAna.h. Referenced by MdaClassify(). |
|
|
Definition at line 42 of file MdaDiscrimAna.h. |
|
|
Definition at line 45 of file MdaDiscrimAna.h. |
|
|
Definition at line 36 of file MdaDiscrimAna.h. |
|
||||||||||||
|
Definition at line 51 of file MdaDiscrimAna.cxx. 00051 : 00052 nueRec(nr), 00053 fMdaDiscrim(md) 00054 { 00055 00056 }
|
|
|
Definition at line 58 of file MdaDiscrimAna.cxx. 00059 {
00060
00061 }
|
|
|
|
|
|
Definition at line 64 of file MdaDiscrimAna.cxx. References FillCalibArrays(), FillPVector(), isFillDone(), isFilled, MdaClassify(), MSG, and PIDCalc(). Referenced by NueRecordAna::Analyze(). 00065 {
00066
00067 MSG("MdaDiscrimAna",Msg::kDebug)<<"In MdaDiscrimAna::Analyze"<<endl;
00068
00069 //Fill SAS Calibration Arrays
00070 if (!isFillDone()) isFilled = FillCalibArrays();
00071
00072 //Fill P-Vector values from NueRecord
00073 FillPVector();
00074
00075 //Calculate nue PID from discriminator function
00076 PIDCalc();
00077
00078 //Classify events according to a probability threshold cut
00079 MdaClassify();
00080
00081 }
|
|
|
Definition at line 83 of file MdaDiscrimAna.cxx. References constCoeff, linCoeff, meanVec, MSG, nClass, quadCoeff, sasFileName, StringUtil::SplitString(), typeClassVec, typeCoeffVec, and varNameVec. Referenced by Analyze(). 00084 {
00085
00086 //Reads the SAS coefficient input file
00087 //(containing quadratic, linear and constant coeffs)
00088 //and fills the necessary matrices and vectors
00089 //Also fills a vector
00090
00091 static const Int_t kLineSize=2000;
00092
00093 ifstream sasInput;
00094 sasInput.open(sasFileName.c_str());
00095 if (sasInput.fail()) {
00096
00097 MSG("MdaDiscrimAna",Msg::kInfo)<< "SAS coefficient file not found! "
00098 << "Not calculating MDA PID..."
00099 << endl;
00100 return false;
00101 }
00102
00103 Char_t oneLine[kLineSize];
00104
00105 string sepaCol=" ";
00106 vector<string> oneVec;
00107
00108 //Read the first line (mean values) and get the number
00109 //of columns in the file
00110 sasInput.getline(oneLine,kLineSize);
00111 string topLine=oneLine;
00112 StringUtil::SplitString(topLine,sepaCol,oneVec);
00113
00114 //Read the mean values into a vector
00115 TVectorD meanTemp(oneVec.size()-2);
00116 for(UInt_t col=2; col < oneVec.size(); col++){
00117 meanTemp(col-2)=(atof((oneVec.at(col)).c_str()));
00118 }
00119
00120 for(Int_t m=0; m < meanTemp.GetNoElements(); m++){
00121 MSG("MdaDiscrimAna",Msg::kDebug)<< "Mean:"
00122 << m
00123 << "="
00124 << meanTemp(m)
00125 << endl;
00126 }
00127
00128 meanVec.push_back(meanTemp);
00129
00130 //Declare matrices and vectors for coefficient reading
00131 TMatrixD quadTemp(oneVec.size()-2, oneVec.size()-2);
00132
00133 TVectorD linTemp(oneVec.size()-2);
00134
00135 TVectorD constTemp(oneVec.size()-2);
00136
00137 deque <string> tempClassVec;
00138
00139 //Count the rows for each matrix
00140 Int_t rowCounter=0;
00141
00142 //Reset the number of classes counter
00143 nClass=0;
00144
00145 //Reset the ifstream pointer to the beginning of the file
00146 // sasInput.seekg(0,ios::beg);
00147
00148 //Reset the line vector
00149 oneVec.clear();
00150
00151 while (! sasInput.eof() ){
00152 sasInput.getline(oneLine,kLineSize);
00153 topLine=oneLine;
00154 StringUtil::SplitString(topLine,sepaCol,oneVec);
00155
00156 string coeffType;
00157
00158 for(UInt_t col=0; col < oneVec.size(); col++){
00159
00160 if (col==0) {
00161
00162 tempClassVec.push_back(oneVec.at(col));
00163
00164 }else if (col==1){
00165
00166 coeffType=oneVec.at(col);
00167
00168 if(coeffType!="_LINEAR_" && coeffType !="_CONST_"){
00169 if (varNameVec.size()<(oneVec.size()-2))
00170 varNameVec.push_back(oneVec.at(col));
00171 typeCoeffVec.push_back("_QUAD_");
00172 }else typeCoeffVec.push_back(oneVec.at(col));
00173
00174 }else{
00175
00176 if (coeffType=="_LINEAR_"){
00177
00178 linTemp(col-2)=(atof((oneVec.at(col)).c_str()));
00179
00180 }else if (coeffType=="_CONST_"){
00181
00182 constTemp(col-2)=(atof((oneVec.at(col)).c_str()));
00183
00184 }else{
00185 quadTemp(rowCounter,(col-2))
00186 =(atof((oneVec.at(col)).c_str()));
00187 }
00188 }
00189 }
00190 oneVec.clear();
00191 rowCounter++;
00192
00193 if (coeffType=="_CONST_"){
00194 nClass++;
00195 rowCounter=0;
00196 typeClassVec.push_back(tempClassVec);
00197 quadCoeff.push_back(quadTemp);
00198 linCoeff.push_back(linTemp);
00199 constCoeff.push_back(constTemp);
00200
00201 tempClassVec.clear();
00202
00203 }else continue;
00204 }
00205
00206 if(sasInput) sasInput.close();
00207
00208 return true;
00209 }
|
|
|
Definition at line 211 of file MdaDiscrimAna.cxx. References count, ANtpDefaultValue::IsDefault(), meanVec, MSG, NeedsSpecialAttention(), varNameVec, and varPVector. Referenced by Analyze(). 00212 {
00213 //Now use varNameVec to fill the variable values p-vector
00214 //Shameless rip-off of CompareAll::FillFromList
00215
00216 NueRecord *nueR=&nueRec;
00217
00218 if(varNameVec.size() == 0) return;
00219 TString hname;
00220 UInt_t count = 0;
00221
00222 TClass *cl;
00223 TRealData *rd;
00224 string vName;
00225 TDataMember *member;
00226 TDataType *membertype;
00227 Double_t value = 0.0;
00228 Int_t valueI = 0;
00229
00230 cl=nueR->IsA();
00231 TIter next(cl->GetListOfRealData());
00232 TVectorD valueTemp(varNameVec.size());
00233
00234 while ((rd =dynamic_cast<TRealData*>(next()))) {
00235 member = rd->GetDataMember();
00236 membertype = member->GetDataType();
00237 vName=rd->GetName();
00238
00239 TString vNameOrig=vName;
00240
00241 //Replace "." with "_" in the ntuple variable names
00242 if (vName.length() > vName.find_first_of("."))
00243 vName.replace(vName.find_first_of("."),1,"_");
00244
00245 Int_t offset = rd->GetThisOffset();
00246 Char_t *pointer = reinterpret_cast<Char_t*>(nueR) + offset;
00247
00248 for(UInt_t i = 0; i < varNameVec.size();i++){
00249 MSG("MdaDiscrimAna",Msg::kDebug)<<"Found variable "
00250 << "NueRec: " << vName
00251 << " SAS: " << varNameVec.at(i)
00252 << endl;
00253 if(vName == varNameVec.at(i)){
00254 value = ANtpDefVal::kDouble;
00255 valueI = ANtpDefVal::kInt;
00256 if(!NeedsSpecialAttention(vNameOrig, nueR, value))
00257
00258 if (!strcmp(membertype->GetTypeName(),"float") ||
00259 !strcmp(membertype->GetTypeName(),"Float_t") ||
00260 !strcmp(membertype->GetTypeName(),"double") ||
00261 !strcmp(membertype->GetTypeName(),"Double_t")){
00262 value=atof(membertype->AsString(pointer));
00263 valueI=1;
00264 }
00265 else if(!strcmp(membertype->GetTypeName(),"int") ||
00266 !strcmp(membertype->GetTypeName(),"Int_t")){
00267 value=atoi(membertype->AsString(pointer));
00268 valueI=atoi(membertype->AsString(pointer));
00269 }
00270
00271 else MSG("MdaDiscrimAna",Msg::kWarning)<<"Found variable "
00272 << "NueRec: " << vName
00273 << " of unknown type "
00274 << membertype->GetTypeName()
00275 << endl ;
00276
00277 MSG("MdaDiscrimAna",Msg::kDebug)<<"Found variable "
00278 <<vName
00279 <<" with value "
00280 <<value
00281 <<endl;
00282
00283 if(!ANtpDefVal::IsDefault(value) &&
00284 !ANtpDefVal::IsDefault(valueI)){
00285
00286 valueTemp(i)=value;
00287
00288 }else {
00289 valueTemp(i)=(meanVec.at(0))(i);
00290 }
00291 MSG("MdaDiscrimAna",Msg::kDebug)<<"Found variable "
00292 <<vName<<" with value "
00293 <<value
00294 <<endl;
00295 count++;
00296 i = varNameVec.size();
00297 }
00298 }
00299 if(count == varNameVec.size()) break;
00300 }
00301
00302 varPVector.push_back(valueTemp);
00303
00304 return;
00305
00306 }
|
|
|
Definition at line 74 of file MdaDiscrimAna.h. Referenced by Analyze(). 00074 {return isFilled; }
|
|
|
Definition at line 310 of file MdaDiscrimAna.cxx. References MdaDiscrim::fMdaClass, fMdaDiscrim, IterDeqDouble_t, MSG, probClass, and typeClassVec. Referenced by Analyze(). 00311 {
00312
00313 //Use the calculated class probabilities and user supplied threshold cut
00314 //to perform posterior event classification
00315
00316 //Determine the maximum probability and its index
00317
00318 if(probClass.size() < static_cast<UInt_t>(nClass)) return;
00319
00320 IterDeqDouble_t maxIter =
00321 max_element(probClass.begin(),probClass.end());
00322
00323 Int_t maxPos= distance(probClass.begin(),maxIter);
00324 MSG("MdaDiscrimAna",Msg::kDebug)<< "Prob Max found for entry "
00325 << maxPos << ": "
00326 << probClass.at(maxPos)
00327 << " with class "
00328 << (typeClassVec.at(maxPos)).at(0)
00329 << endl;
00330
00331 //Classify into class ()
00332 Double_t probMax=probClass.at(maxPos);
00333 string classMax=(typeClassVec.at(maxPos)).at(0);
00334
00335 if (probMax < threshCut) fMdaDiscrim.fMdaClass=-1;
00336 else if (classMax=="nue") fMdaDiscrim.fMdaClass=ClassType::nue;
00337 else if (classMax=="ncu") fMdaDiscrim.fMdaClass=ClassType::NC;
00338 else if (classMax=="num") fMdaDiscrim.fMdaClass=ClassType::numu;
00339 else if (classMax=="nut") fMdaDiscrim.fMdaClass=ClassType::nutau;
00340 else {
00341 MSG("MdaDiscrimAna",Msg::kError)
00342 <<"Unknown class in SAS calibration set" << endl;
00343 fMdaDiscrim.fMdaClass=ANtpDefVal::kInt;
00344 }
00345 return;
00346
00347 }
|
|
||||||||||||||||
|
Definition at line 350 of file MdaDiscrimAna.cxx. References MSTCalc::eallm1, MSTCalc::eallw1, MSTCalc::enn1, NueHeader::GetEventNo(), NueHeader::GetEvents(), RecRecordImp< T >::GetHeader(), NueHeader::GetRun(), NueHeader::GetSnarl(), NueHeader::GetSubRun(), NueHeader::GetTrackLength(), NueRecord::mstvars, MSTCalc::oallm1, MSTCalc::oallw1, and MSTCalc::onn1. Referenced by FillPVector(). 00353 {
00354
00355 //All the fHeaders and four of the MST vars require special effort
00356 if(name == "fHeader.fSnarl") {
00357 value = nr->GetHeader().GetSnarl();
00358 }if(name == "fHeader.fRun") {
00359 value = nr->GetHeader().GetRun();
00360 }if(name == "fHeader.fSubRun") {
00361 value = nr->GetHeader().GetSubRun();
00362 }if(name == "fHeader.fEvtNo") {
00363 value = nr->GetHeader().GetEventNo();
00364 }if(name == "fHeader.fEvents") {
00365 value = nr->GetHeader().GetEvents();
00366 }if(name == "fHeader.fTrackLength") {
00367 value = nr->GetHeader().GetTrackLength();
00368 }
00369
00370 if(name == "mstvars.eallw1") {
00371 if(nr->mstvars.enn1 > 0) value = 0.0;
00372 for(int i=0;i<nr->mstvars.enn1;i++){
00373 value += nr->mstvars.eallw1[i];
00374 }
00375 }
00376 if(name == "mstvars.oallw1") {
00377 if(nr->mstvars.onn1 > 0) value = 0.0;
00378 for(int i=0;i<nr->mstvars.onn1;i++){
00379 value += nr->mstvars.oallw1[i];
00380 }
00381 }
00382 if(name == "mstvars.eallm1") {
00383 if(nr->mstvars.enn1 > 0) value = 0.0;
00384 for(int i=0;i<nr->mstvars.enn1;i++){
00385 value += nr->mstvars.eallm1[i];
00386 }
00387 }
00388 if(name == "mstvars.oallm1") {
00389 if(nr->mstvars.onn1 > 0) value = 0.0;
00390 for(int i=0;i<nr->mstvars.onn1;i++){
00391 value += nr->mstvars.oallm1[i];
00392 }
00393 }
00394
00395 if(value > -9999) return true;
00396 return false;
00397 }
|
|
|
|
|
|
Definition at line 400 of file MdaDiscrimAna.cxx. References constCoeff, DeqDouble_t, discrimDiff, fMdaDiscrim, MdaDiscrim::fMdaPIDnc, MdaDiscrim::fMdaPIDnue, MdaDiscrim::fMdaPIDnumu, MdaDiscrim::fMdaPIDnutau, linCoeff, MSG, probClass, quadCoeff, typeClassVec, and varPVector. Referenced by Analyze(). 00401 {
00402
00403 DeqDouble_t discrimValue;
00404
00405 //Calculate discriminant function Q for each class
00406 for (Int_t n=0; n < nClass; n++){
00407
00408 //Create an auxiliary copy of the Pvector
00409 TVectorD varPClone=(varPVector.at(0));
00410
00411 //Note that varPClone is modified by the "*=" operator
00412 Double_t quadTerm=((varPClone)*=(quadCoeff.at(n)))*(varPVector.at(0));
00413
00414 Double_t linTerm=(linCoeff.at(n))*(varPVector.at(0));
00415
00416 Double_t constTerm=(constCoeff.at(n))(0);
00417
00418 discrimValue.push_back(quadTerm+linTerm+constTerm);
00419
00420 MSG("MdaDiscrimAna",Msg::kDebug)<< "Q("
00421 << (typeClassVec.at(n)).at(0)
00422 <<") = "
00423 << discrimValue.at(n) <<endl;
00424 }
00425
00426 TVectorD tempDiff(nClass-1);
00427
00428 //Calculate discriminator difference vectors
00429 //P(a) = e^a/(e^(a)+e^(b)+e^(c)+...) = 1/(1+e^(b-a)+e^(c-a)+...)
00430 //Avoids e^X where X is big => machine-dep floating exception.
00431
00432 for (Int_t n=0; n < nClass; n++){
00433
00434 Int_t diffCount=0;
00435
00436 for (Int_t m=0; m < nClass; m++){
00437 if(m!=n){
00438 tempDiff(diffCount)=(discrimValue.at(m)-discrimValue.at(n));
00439 diffCount++;
00440 } else continue;
00441 }
00442
00443 if (TMath::Abs(tempDiff.Min())>600. || TMath::Abs(tempDiff.Max())>600){
00444 MSG("MdaDiscrimAna",Msg::kDebug)<<"Bypass floating exception "
00445 << endl;
00446 continue;
00447 }else discrimDiff.push_back(tempDiff);
00448
00449 tempDiff.Zero();
00450 }
00451
00452 if(discrimDiff.size()!=static_cast<UInt_t>(nClass)) return;
00453
00454 //Calculate probability value for each class
00455 for (Int_t n=0; n < nClass; n++){
00456
00457 Double_t denomSum=0.0;
00458
00459 for (Int_t m=0; m < (nClass-1); m++){
00460 denomSum+=TMath::Exp((discrimDiff.at(n))(m));
00461 }
00462
00463 probClass.push_back(1.0/(1.0+denomSum));
00464
00465 MSG("MdaDiscrimAna",Msg::kDebug) << "P("
00466 << (typeClassVec.at(n)).at(0)
00467 <<") = "
00468 << probClass.at(n) <<endl;
00469
00470 if((typeClassVec.at(n)).at(0)=="nue")
00471 fMdaDiscrim.fMdaPIDnue=probClass.at(n);
00472 if((typeClassVec.at(n)).at(0)=="ncu")
00473 fMdaDiscrim.fMdaPIDnc=probClass.at(n);
00474 if((typeClassVec.at(n)).at(0)=="num")
00475 fMdaDiscrim.fMdaPIDnumu=probClass.at(n);
00476 if((typeClassVec.at(n)).at(0)=="nut")
00477 fMdaDiscrim.fMdaPIDnutau=probClass.at(n);
00478 }
00479
00480 return;
00481 }
|
|
|
Definition at line 484 of file MdaDiscrimAna.cxx. 00485 {
00486
00487
00488 }
|
|
||||||||||||
|
Definition at line 59 of file MdaDiscrimAna.h. References sasFileName, and threshCut. Referenced by NueModule::Analyze(). 00059 {
00060 threshCut=threshC;
00061 sasFileName=sasFN;
00062 }
|
|
|
Definition at line 44 of file MdaDiscrimAna.cxx. Referenced by FillCalibArrays(), and PIDCalc(). |
|
|
Definition at line 67 of file MdaDiscrimAna.h. Referenced by PIDCalc(). |
|
|
Definition at line 79 of file MdaDiscrimAna.h. Referenced by MdaClassify(), and PIDCalc(). |
|
|
Definition at line 308 of file MdaDiscrimAna.cxx. Referenced by Analyze(). |
|
|
Definition at line 43 of file MdaDiscrimAna.cxx. Referenced by FillCalibArrays(), and PIDCalc(). |
|
|
Definition at line 45 of file MdaDiscrimAna.cxx. Referenced by FillCalibArrays(), and FillPVector(). |
|
|
Definition at line 41 of file MdaDiscrimAna.cxx. Referenced by FillCalibArrays(). |
|
|
Definition at line 78 of file MdaDiscrimAna.h. |
|
|
Definition at line 69 of file MdaDiscrimAna.h. Referenced by MdaClassify(), and PIDCalc(). |
|
|
Definition at line 42 of file MdaDiscrimAna.cxx. Referenced by FillCalibArrays(), and PIDCalc(). |
|
|
Definition at line 98 of file MdaDiscrimAna.h. Referenced by FillCalibArrays(), and SetMdaParams(). |
|
|
Definition at line 97 of file MdaDiscrimAna.h. Referenced by SetMdaParams(). |
|
|
Definition at line 47 of file MdaDiscrimAna.cxx. Referenced by FillCalibArrays(), MdaClassify(), and PIDCalc(). |
|
|
Definition at line 48 of file MdaDiscrimAna.cxx. Referenced by FillCalibArrays(). |
|
|
Definition at line 49 of file MdaDiscrimAna.cxx. Referenced by FillCalibArrays(), and FillPVector(). |
|
|
Definition at line 65 of file MdaDiscrimAna.h. Referenced by FillPVector(), and PIDCalc(). |
1.3.9.1