loon: NueData.cxx Source File

00001 
00002 //$Id: NueData.cxx,v 1.12 2009/08/27 15:20:01 scavan Exp $
00003 //
00004 //NueData
00005 //
00006 //J Boehm 12/2007
00008                                                                                                                                                       
00009 #include "NueAna/Extrapolation/NueData.h"
00010 #include "MessageService/MsgService.h"
00011 #include "TMath.h"
00012 #include "NueAna/NueRecord.h"
00013 #include "NueAna/NueMini.h"
00014 #include "NueAna/NueMiniAna.h"
00015 
00016 CVSID("$Id: NueData.cxx,v 1.12 2009/08/27 15:20:01 scavan Exp $");                                                                                                                                                      
00017 //-------------------------------------------------------------------
00018 NueData::NueData() :
00019   fBeam(BeamType::kL010z185i),
00020   fDet(Detector::kFar),
00021   //  fRelease(ReleaseType::kCedarPhyDaikon)
00022   fRelease(ReleaseType::kDogwood1Daikon)
00023 {
00024  fPOT = 0.0;
00025  kIsMR = false;
00026  fNM = new NueMini(fBeam, fDet, fRelease);
00027 }
00028 
00029 NueData::NueData(BeamType::BeamType_t beam,
00030              Detector::Detector_t det,
00031              ReleaseType::Release_t rel, bool isNueData) :
00032   fBeam(beam),
00033   fDet(det),
00034   fRelease(rel),
00035   kIsNueData(isNueData)
00036 {
00037   fPOT = 0.0;
00038   kIsMR = false;
00039   fNM = new NueMini(fBeam, fDet, fRelease);
00040 }
00041 
00042 void NueData::SetupNueHeader(NueHeader &nh)
00043 {
00044    const VldTimeStamp vts;
00045    const VldContext Vld(fDet, SimFlag::kMC, vts);
00046    nh.SetVldContext(Vld);
00047    nh.SetRelease(fRelease);
00048    nh.SetBeamType(fBeam);
00049    nh.SetFoundBits(0,0,0,kIsMR);
00050 }
00051 
00052 void NueData::AddEvent(NueRecord *nr)
00053 {
00054    if(!nr){
00055      MSG("NueData", Msg::kError) << "No NueRecord Passed"<<endl;
00056      return;
00057    }
00058 
00059    if(nr->GetHeader().GetRelease() != fRelease){
00060      MSG("NueData", Msg::kError) << "ReleaseType MisMatch" 
00061             << nr->GetHeader().GetRelease()<<" != "<<fRelease<<endl;
00062    }
00063    if(nr->GetHeader().GetBeamType() != fBeam){
00064      MSG("NueData", Msg::kError) << "BeamType MisMatch"
00065             << nr->GetHeader().GetBeamType()<<" != "<<fBeam<<endl;
00066    }
00067    if(nr->GetHeader().GetVldContext().GetDetector() != fDet){
00068      MSG("NueData", Msg::kError) << "Detector MisMatch"
00069             << nr->GetHeader().GetVldContext().GetDetector() 
00070             <<" != "<<fDet<<endl;
00071    }
00072 
00073    fNMA.FillMini(nr, fNM);
00074    AddEvent(fNM);
00075 
00076 /* 
00077    evtRecoNueEnergy.push_back(nr->srevent.phNueGeV);
00078    evtRecoMEUEnergy.push_back(nr->srevent.phMip);
00079 
00080    int tP = 0;    int tLP = 0;     float trkECC = 0;
00081    if(nr->srtrack.phCCGeV > -10){
00082      trkECC = nr->srtrack.phCCGeV;
00083      tP = TMath::Abs(nr->srtrack.endPlane - nr->srtrack.begPlane);
00084      tLP = nr->srtrack.trklikePlanes;
00085    }
00086    trkRecoCCEnergy.push_back(trkECC);
00087    trkPlanes.push_back(tP);
00088    trkEndPlane.push_back(nr->srtrack.endPlane);
00089    trkBegPlane.push_back(nr->srtrack.begPlane);
00090    trkLikePlanes.push_back(tLP);
00091 
00092    nshower.push_back(nr->srevent.showers);
00093    contPlanes.push_back(nr->shwfit.contPlaneCount050);
00094 
00095    if(nr->srshower.phCCGeV > -10)
00096      shwRecoCCEnergy.push_back(nr->srshower.phCCGeV);
00097    else shwRecoCCEnergy.push_back(0);
00098 
00099    ann2pe.push_back(nr->ann.pid_11inp);
00100    ann30.push_back(nr->ann.pid_30inp);
00101    ann6.push_back(nr->ann.pid_6inp);
00102 
00103    cutPID.push_back(nr->treepid.fCutPID);
00104    ssPID.push_back(nr->subshowervars.pid);
00105    mcnnPID.push_back(nr->mcnnv.mcnn_var1);
00106    mcnnMatch.push_back(nr->mcnnv.bestmatches);
00107 
00108    //truth values
00109    shiEpi0.push_back(nr->shi.epi0);
00110    shiEmEnergy.push_back(nr->shi.emenergy);
00111    nuEnergy.push_back(nr->mctrue.nuEnergy);
00112    nuFlavor.push_back(nr->mctrue.nuFlavor);
00113    nonOscNuFlavor.push_back(nr->mctrue.nonOscNuFlavor);
00114    nueClass.push_back(nr->mctrue.fNueClass);
00115 
00116    nueOscProb.push_back(nr->mctrue.fOscProb);
00117    interactionType.push_back(nr->mctrue.interactionType);
00118    nuDCosX.push_back(nr->mctrue.nuDCosX);
00119    nuDCosY.push_back(nr->mctrue.nuDCosY);
00120    nuDCosZ.push_back(nr->mctrue.nuDCosZ);
00121    hadronicY.push_back(nr->mctrue.hadronicY);
00122    hadronicFinalState.push_back(nr->mctrue.hadronicFinalState);
00123    w2.push_back(nr->mctrue.w2);
00124    q2.push_back(nr->mctrue.q2);
00125    bjorkenX.push_back(nr->mctrue.bjorkenX);
00126    targetPX.push_back(nr->mctrue.targetPX);
00127    targetPY.push_back(nr->mctrue.targetPY);
00128    targetPZ.push_back(nr->mctrue.targetPZ);
00129    targetEnergy.push_back(nr->mctrue.targetEnergy);
00130    atomicNumber.push_back(nr->mctrue.atomicNumber);
00131    atomicWeight.push_back(nr->mctrue.atomicWeight);
00132    initialState.push_back(nr->mctrue.initialState);
00133    resonanceCode.push_back(nr->mctrue.resonanceCode);
00134 
00135    skzpWeight.push_back(nr->fluxweights.totbeamweight);
00136 
00137 
00138    ntrack.push_back(nr->srevent.tracks);
00139    trkPass.push_back(nr->srtrack.passedFit);
00140    endPlaneU.push_back(nr->srtrack.endPlaneU);
00141    endPlaneV.push_back(nr->srtrack.endPlaneV);
00142    deltaUVVtx.push_back(nr->srtrack.deltaUVVtx);
00143    abCCPID.push_back(nr->anainfo.abCCPID);
00144    roCCPID.push_back(nr->anainfo.roCCPID);
00145 
00146    mri_roCCPID.push_back(nr->mri.orig_roCCPID);
00147    mri_abCCPID.push_back(nr->mri.orig_abCCPID);
00148    mri_trkPass.push_back(nr->mri.fitp);
00149    gapPlanes.push_back(nr->srshower.gapPlanes);
00150    neugenStdXsec.push_back(1);
00151    cosmicCut.push_back((int) nr->dtree.bt_var1);
00152    largestEvent.push_back(nr->srevent.largestEvent);
00153 */
00154 }
00155 
00156 void NueData::AddEvent(NueMini *nm)
00157 {
00158    if(!nm){
00159      MSG("NueData", Msg::kError) << "No NueRecord Passed"<<endl;
00160      return;
00161    }
00162                                                                                                                             
00163    if(nm->fRelease != fRelease){
00164      MSG("NueData", Msg::kError) << "ReleaseType MisMatch"
00165             << nm->fRelease<<" != "<<fRelease<<endl;
00166    }
00167    if(nm->fBeam != fBeam){
00168      MSG("NueData", Msg::kError) << "BeamType MisMatch"
00169             << nm->fBeam<<" != "<<fBeam<<endl;
00170    }
00171    if(nm->fDet != fDet){
00172      MSG("NueData", Msg::kError) << "Detector MisMatch"
00173             << nm->fDet <<" != "<<fDet<<endl;
00174    }                                                                                                                          
00175    evtRecoNueEnergy.push_back(nm->evtRecoNueEnergy);
00176    evtRecoMEUEnergy.push_back(nm->evtRecoMEUEnergy);
00177    trkRecoCCEnergy.push_back(nm->trkRecoCCEnergy);
00178    trkPlanes.push_back(nm->trkPlanes);
00179    trkEndPlane.push_back(nm->trkEndPlane);
00180    trkBegPlane.push_back(nm->trkBegPlane);
00181    trkLikePlanes.push_back(nm->trkLikePlanes);
00182    shwRecoCCEnergy.push_back(nm->shwRecoCCEnergy);
00183 
00184    nshower.push_back(nm->nshower);
00185    contPlanes.push_back(nm->contPlanes);
00186 
00187    //ann2pe_daikon04.push_back(nm->ann2pe_daikon04);
00188    //ann2pe.push_back(nm->ann2pe);
00189    
00190    printf("we've adjusted the NueMini structure (changed ann var names)... but didn't fix NueData.cxx... do that now... (ask Steve)\n");
00191    exit(1);
00192    
00193    //ann30.push_back(nm->ann30);
00194    //ann6.push_back(nm->ann6);
00195    ssPID.push_back(nm->ssPID);
00196    //cutPID.push_back(nm->cutPID);                                                                                                                            
00197    mcnnPID.push_back(nm->mcnnPID);
00198    mcnnMatch.push_back(nm->mcnnMatch);
00199 
00200    //truth values
00201    shiEpi0.push_back(nm->shiEpi0);
00202    shiEmEnergy.push_back(nm->shiEmEnergy);
00203    nuEnergy.push_back(nm->nuEnergy);
00204    nuFlavor.push_back(nm->nuFlavor);
00205    nonOscNuFlavor.push_back(nm->nonOscNuFlavor);
00206    nueClass.push_back(nm->nueClass);
00207                                                                                                                             
00208    nueOscProb.push_back(nm->nueOscProb);
00209    interactionType.push_back(nm->interactionType);
00210    nuDCosX.push_back(nm->nuDCosX);
00211    nuDCosY.push_back(nm->nuDCosY);
00212    nuDCosZ.push_back(nm->nuDCosZ);
00213    hadronicY.push_back(nm->hadronicY);
00214    hadronicFinalState.push_back(nm->hadronicFinalState);
00215    w2.push_back(nm->w2);
00216    q2.push_back(nm->q2);
00217    bjorkenX.push_back(nm->bjorkenX);
00218    targetPX.push_back(nm->targetPX);
00219    targetPY.push_back(nm->targetPY);
00220    targetPZ.push_back(nm->targetPZ);
00221    targetEnergy.push_back(nm->targetEnergy);
00222    atomicNumber.push_back(nm->atomicNumber);
00223    atomicWeight.push_back(nm->atomicWeight);
00224    initialState.push_back(nm->initialState);
00225    resonanceCode.push_back(nm->resonanceCode);
00226                 
00227 //the meaning of skzpWeight might have changed in nuemini!
00228 
00229 //it seems to have been RPWeight....
00230 // now it is actually skzpWeight .....
00231 // RPWeight is now nm->RPWeight                                                                                                            
00232    //skzpWeight.push_back(nm->skzpWeight);
00233         exit(1);
00234    ntrack.push_back(nm->ntrack);
00235    trkPass.push_back(nm->trkPass);
00236    endPlaneU.push_back(nm->endPlaneU);
00237    endPlaneV.push_back(nm->endPlaneV);
00238    deltaUVVtx.push_back(nm->deltaUVVtx);
00239    abCCPID.push_back(nm->abCCPID);
00240    roCCPID.push_back(nm->roCCPID);
00241    mri_roCCPID.push_back(nm->mri_orig_roCCPID);
00242    mri_abCCPID.push_back(nm->mri_orig_abCCPID);
00243    mri_trkPass.push_back(nm->mri_trkPass);   
00244    gapPlanes.push_back(nm->gapPlanes);
00245    neugenStdXsec.push_back(1);
00246    cosmicCut.push_back((int) nm->cosmicCut);
00247    largestEvent.push_back(nm->largestEvent);
00248 
00249 }
00250 
00251 
00252 
00253 void NueData::FillRecord(NueRecord *nr, int i)
00254 {  
00255    if(!nr){
00256      MSG("NueData", Msg::kError) << "No NueRecord Passed"<<endl;
00257      return;
00258    }
00259 
00260    if(nr->GetHeader().GetVldContext().GetDetector() != fDet){
00261      MSG("NueData", Msg::kError) << "Detector MisMatch"
00262             << nr->GetHeader().GetVldContext().GetDetector()
00263             <<" != "<<fDet<<endl;
00264    }
00265    if(nr->GetHeader().GetRelease() != fRelease){
00266      MSG("NueData", Msg::kError) << "ReleaseType MisMatch"
00267             << nr->GetHeader().GetRelease()<<" != "<<fRelease<<endl;
00268    }
00269    if(nr->GetHeader().GetBeamType() != fBeam){
00270      MSG("NueData", Msg::kError) << "BeamType MisMatch"
00271             << nr->GetHeader().GetBeamType()<<" != "<<fBeam<<endl;
00272    }
00273 
00274    nr->srevent.phNueGeV = evtRecoNueEnergy[i];
00275    nr->srevent.phMip = evtRecoMEUEnergy[i];
00276    nr->srtrack.phCCGeV = trkRecoCCEnergy[i];
00277    nr->srtrack.trklikePlanes = trkLikePlanes[i];
00278    nr->srtrack.begPlane = trkBegPlane[i];
00279    nr->srtrack.endPlane = trkEndPlane[i];
00280    nr->srshower.phCCGeV = shwRecoCCEnergy[i];
00281    nr->srevent.showers = nshower[i];
00282    nr->shwfit.contPlaneCount050 = contPlanes[i];
00283    nr->ann.pid_11inp_daikon04 = ann2pe_daikon04[i];
00284    nr->ann.pid_11inp = ann2pe[i];
00285                                                                                 
00286    nr->ann.pid_30inp = ann30[i];
00287    nr->ann.pid_6inp = ann6[i];
00288 
00289    nr->subshowervars.pid = ssPID[i];
00290    nr->treepid.fCutPID = cutPID[i];   
00291    nr->mcnnv.mcnn_var1 = mcnnPID[i];
00292    nr->mcnnv.bestmatches = mcnnMatch[i];
00293 
00294                                                                              
00295    //truth values
00296    nr->shi.epi0 = shiEpi0[i];
00297    nr->shi.emenergy = shiEmEnergy[i];
00298    nr->mctrue.nuEnergy = nuEnergy[i];
00299    nr->mctrue.nuFlavor = nuFlavor[i];
00300    nr->mctrue.nonOscNuFlavor = nonOscNuFlavor[i];
00301    nr->mctrue.fNueClass = nueClass[i];
00302                                                                                 
00303    nr->mctrue.fOscProb = nueOscProb[i];
00304    nr->mctrue.interactionType = interactionType[i];
00305    nr->mctrue.nuDCosX = nuDCosX[i];
00306    nr->mctrue.nuDCosY = nuDCosY[i];
00307    nr->mctrue.nuDCosZ = nuDCosZ[i];
00308    nr->mctrue.hadronicY = hadronicY[i];
00309    nr->mctrue.hadronicFinalState = hadronicFinalState[i];
00310    nr->mctrue.w2 = w2[i];
00311    nr->mctrue.q2 = q2[i];
00312    nr->mctrue.bjorkenX = bjorkenX[i];
00313    nr->mctrue.targetPX = targetPX[i];
00314    nr->mctrue.targetPY = targetPY[i];
00315    nr->mctrue.targetPZ = targetPZ[i];
00316    nr->mctrue.targetEnergy = targetEnergy[i];
00317    nr->mctrue.atomicNumber = atomicNumber[i];
00318    nr->mctrue.atomicWeight = atomicWeight[i];
00319    nr->mctrue.initialState = initialState[i];
00320    nr->mctrue.resonanceCode = resonanceCode[i];
00321                                                                                 
00322    nr->fluxweights.totbeamweight = skzpWeight[i];
00323 
00324    nr->srevent.tracks = ntrack[i];
00325    nr->srtrack.passedFit = trkPass[i];
00326    nr->srtrack.endPlaneU = endPlaneU[i];
00327    nr->srtrack.endPlaneV = endPlaneV[i];
00328    nr->srtrack.deltaUVVtx = deltaUVVtx[i];
00329    nr->anainfo.abCCPID = abCCPID[i];
00330    nr->xsecweights.xsecweight = neugenStdXsec[i];
00331    nr->mri.orig_abCCPID = mri_abCCPID[i];
00332    nr->mri.fitp = mri_trkPass[i];
00333    nr->srshower.gapPlanes = gapPlanes[i];
00334    nr->eventq.passCosmicCut = 1; // cosmicCut[i]; this is a hack for now due to cheer problems
00335    nr->srevent.largestEvent = largestEvent[i];
00336 
00337    nr->anainfo.roCCPID = roCCPID[i];
00338    nr->mri.orig_roCCPID = mri_roCCPID[i];
00339 
00340 }
00341 
00342 double NueData::GetNeugenStdXsec(int i)
00343 {
00344    return neugenStdXsec[i];
00345 }
00346 void   NueData::SetNeugenStdXsec(double xsec, int i){
00347   neugenStdXsec[i] = xsec;
00348 }
00349 
00350 void NueData::Clear()
00351 {
00352   trkRecoCCEnergy.clear();
00353   shwRecoCCEnergy.clear();
00354   evtRecoNueEnergy.clear();
00355   evtRecoMEUEnergy.clear();
00356   skzpWeight.clear();
00357   weight.clear();
00358   trkPlanes.clear();
00359   trkEndPlane.clear();
00360   trkBegPlane.clear();
00361   trkLikePlanes.clear();
00362   shiEpi0.clear();
00363   shiEmEnergy.clear();
00364   nuEnergy.clear();
00365   nuFlavor.clear();
00366   nonOscNuFlavor.clear();
00367   nueOscProb.clear();
00368   nueClass.clear();
00369   ParentType.clear();
00370   interactionType.clear();
00371   nuDCosX.clear();
00372   nuDCosY.clear();
00373   nuDCosZ.clear();
00374   hadronicY.clear();
00375   hadronicFinalState.clear();
00376   w2.clear();
00377   q2.clear();
00378   bjorkenX.clear();
00379   targetPX.clear();
00380   targetPY.clear();
00381   targetPZ.clear();
00382   targetEnergy.clear();
00383   atomicNumber.clear();
00384   atomicWeight.clear();
00385   initialState.clear();
00386   resonanceCode.clear();
00387   ann30.clear();
00388   ann6.clear();
00389   ssPID.clear();
00390   mcnnPID.clear();
00391   mcnnMatch.clear();
00392   cutPID.clear();
00393   abCCPID.clear();
00394   ntrack.clear();
00395   trkPass.clear();
00396   endPlaneU.clear();
00397   endPlaneV.clear();
00398   deltaUVVtx.clear();
00399   mri_abCCPID.clear();
00400   mri_trkPass.clear();
00401   gapPlanes.clear();
00402   neugenStdXsec.clear();
00403   cosmicCut.clear();
00404   largestEvent.clear();
00405 
00406   contPlanes.clear();
00407   ann2pe_daikon04.clear();
00408   ann2pe.clear();
00409   nshower.clear();
00410   roCCPID.clear();
00411   mri_roCCPID.clear();
00412 
00413 
00414   vector<double>(trkRecoCCEnergy).swap(trkRecoCCEnergy);
00415   vector<double>(shwRecoCCEnergy).swap(shwRecoCCEnergy);
00416   vector<double>(evtRecoNueEnergy).swap(evtRecoNueEnergy);
00417   vector<double>(evtRecoMEUEnergy).swap(evtRecoMEUEnergy);
00418   vector<double>(skzpWeight).swap(skzpWeight);
00419   vector<double>(weight).swap(weight);
00420   vector<int>(trkPlanes).swap(trkPlanes);
00421   vector<int>(trkEndPlane).swap(trkEndPlane);
00422   vector<int>(trkBegPlane).swap(trkBegPlane);
00423   vector<int>(trkLikePlanes).swap(trkLikePlanes);      
00424   vector<double>(shiEpi0).swap(shiEpi0);
00425   vector<double>(shiEmEnergy).swap(shiEmEnergy);
00426   vector<double>(nuEnergy).swap(nuEnergy);
00427   vector<int>(nuFlavor).swap(nuFlavor);
00428   vector<int>(nonOscNuFlavor).swap(nonOscNuFlavor);
00429   vector<double>(nueOscProb).swap(nueOscProb);
00430   vector<int>(nueClass).swap(nueClass);
00431   vector<double>(ParentType).swap(ParentType);
00432   vector<int>(interactionType).swap(interactionType);
00433   vector<double>(nuDCosX).swap(nuDCosX);
00434   vector<double>(nuDCosY).swap(nuDCosY);
00435   vector<double>(nuDCosZ).swap(nuDCosZ);
00436   vector<double>(hadronicY).swap(hadronicY);
00437   vector<int>(hadronicFinalState).swap(hadronicFinalState);
00438 
00439   vector<double>(w2).swap(w2);
00440   vector<double>(q2).swap(q2);
00441   vector<double>(bjorkenX).swap(bjorkenX);
00442   vector<double>(targetPX).swap(targetPX);
00443   vector<double>(targetPY).swap(targetPY);
00444   vector<double>(targetPZ).swap(targetPZ);
00445   vector<double>(targetEnergy).swap(targetEnergy);
00446   vector<double>(atomicNumber).swap(atomicNumber);
00447   vector<double>(atomicWeight).swap(atomicWeight);
00448   vector<int>(initialState).swap(initialState);
00449   vector<int>(resonanceCode).swap(resonanceCode);
00450   vector<double>(ann30).swap(ann30);
00451   vector<double>(ann6).swap(ann6);
00452   vector<double>(ssPID).swap(ssPID);
00453   vector<double>(mcnnPID).swap(mcnnPID);
00454   vector<int>(mcnnMatch).swap(mcnnMatch);
00455   vector<int>(cutPID).swap(cutPID);
00456   vector<double>(abCCPID).swap(abCCPID);
00457   vector<int>(ntrack).swap(ntrack);
00458   vector<int>(trkPass).swap(trkPass);
00459   vector<int>(endPlaneU).swap(endPlaneU);
00460   vector<int>(endPlaneV).swap(endPlaneV);
00461   vector<int>(deltaUVVtx).swap(deltaUVVtx);
00462   vector<double>(mri_abCCPID).swap(mri_abCCPID);
00463   vector<int>(mri_trkPass).swap(mri_trkPass);
00464   vector<int>(gapPlanes).swap(gapPlanes);
00465   vector<double>(neugenStdXsec).swap(neugenStdXsec);
00466   vector<int>(cosmicCut).swap(cosmicCut);
00467   vector<int>(largestEvent).swap(largestEvent);
00468 
00469   vector<int>(contPlanes).swap(contPlanes); 
00470   vector<double>(ann2pe_daikon04).swap(ann2pe_daikon04); 
00471   vector<double>(ann2pe).swap(ann2pe); 
00472   vector<int>(nshower).swap(nshower); 
00473   vector<double>(roCCPID).swap(roCCPID); 
00474   vector<double>(mri_roCCPID).swap(mri_roCCPID); 
00475  
00476 
00477 
00478 
00479 }
00480