loon: FCPCFilterModule.cxx Source File

00001 
00002 //
00003 //FCPCFilterModule.cxx
00004 // A first pass Partially and fully contained Filter
00005 // John Chapman chapman@hep.phy.cam.ac.uk
00006 // (Original Code by Caius Howcroft)
00008 
00009 #include "FCPCFilter/FCPCFilterModule.h"
00010 #include "MessageService/MsgService.h"
00011 #include "MinosObjectMap/MomNavigator.h"
00012 #include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
00013 #include "UgliGeometry/UgliGeomHandle.h"
00014 #include "UgliGeometry/UgliStripHandle.h"
00015 #include "Validity/VldContext.h"
00016 #include "CandData/CandHeader.h"
00017 #include "CandData/CandRecord.h"
00018 #include "CandDigit/CandDigitListHandle.h"
00019 #include "CandDigit/CandDigitHandle.h"
00020 #include "CandDigit/CandDeMuxDigitHandle.h"
00021 #include "CandDigit/CandDeMuxDigitListHandle.h"
00022 #include "TFile.h"
00023 #include "TTree.h"
00024 #include "TDirectory.h"
00025 #include <string>
00026 using std::cin;
00027 using std::string;
00028 #include <cassert>
00029 #include <cmath>
00030 
00031 JOBMODULE(FCPCFilterModule, "FCPCFilterModule",
00032           "Filters out PC or FC events");
00033 CVSID("$Id: FCPCFilterModule.cxx,v 1.8 2007/03/01 17:42:09 rhatcher Exp $");
00034 //......................................................................
00035 
00036 FCPCFilterModule::FCPCFilterModule():
00037   cDoPC(0),
00038   cDoFC(1),
00039   cDoUC(0),
00040   cDoUp(0),
00041   cListIn("demuxdigitlist"),
00042   cMinDmxWeight(0.25),
00043   cMinDigitCharge(0.0),
00044   cFidRadius(0.3),
00045   cMinFidCharge(6.0),
00046 
00047   cSM1StartPlane(1),
00048   cSM1StopPlane(248),
00049   cSM2StartPlane(250),
00050   cSM2StopPlane(485),
00051 
00052   cMinPlaneCut(7),
00053   cMinPlaneCharge(3.0),
00054   cEndPlaneVeto(4),
00055   cEndPlaneVetoCharge(6.0),
00056   cMaxFidPctChargeCut(-1.0),
00057   cMaxRegionSep(0.5),
00058 
00059   cMakeTree(0),
00060   fFileName("FCPCFilterOutput.root")
00061 
00062 
00063 {
00064   ResetParams();
00065 //======================================================================
00066 // 
00067 //======================================================================
00068 }
00069 
00070 //......................................................................
00071 
00072 FCPCFilterModule::~FCPCFilterModule()
00073 {
00074 //======================================================================
00075 // FILL_IN: [Document your code!!]
00076 //======================================================================
00077 }
00078 
00079 //......................................................................
00080 
00081 void FCPCFilterModule::BeginJob()
00082 {
00083 //======================================================================
00084 // FILL_IN: [Document your code!!]
00085 //======================================================================
00086   if(cMakeTree==1)
00087     {
00088       TDirectory* current_dir = gDirectory;
00089       fFile = new TFile(fFileName.Data(),"recreate");
00090       fTree = new TTree("mu", "Muon Properties");
00091       fTree->Branch("run", &fRun, "run/I");
00092       fTree->Branch("snarl", &fSnarl, "snarl/I");
00093       fTree->Branch("Con",&fEventClass, "Con/I");
00094       fTree->Branch("TotalCharge", &fCharge, "TotalCharge/D");//total charge used (see above cuts)
00095       fTree->Branch("NPlane",&fNPlane, "NPlane/I");//Number of planes above some charge treshold
00096       fTree->Branch("SideCharge",fFidCharge,"SideCharge[8]/D"); 
00097       fTree->Branch("TotSideCharge", &fSumFidCharge, "TotSideCharge/D");//the total charge on all sides
00098       fTree->Branch("bsm1q",&fSM1begPlaneQ,"bsm1q/D");
00099       fTree->Branch("esm1q",&fSM1endPlaneQ,"esm1q/D");
00100       fTree->Branch("bsm2q",&fSM2begPlaneQ,"bsm2q/D");
00101       fTree->Branch("esm2q",&fSM2endPlaneQ,"esm2q/D");
00102       fTree->Branch("fidsm1q",&fChargeInSM1,"fidsm1q/D");
00103       fTree->Branch("fidsm2q",&fChargeInSM2,"fidsm2q/D");
00104       fTree->Branch("regU",fFidU,"regU[12]/D");
00105       fTree->Branch("regV",fFidV,"regV[12]/D");
00106       fTree->Branch("regZ",fFidZ,"regZ[12]/D");
00107 
00108       fTree->SetDirectory(fFile);
00109       gDirectory = current_dir;
00110     }
00111 }
00112 
00113 //......................................................................
00114 
00115 void FCPCFilterModule::EndJob()
00116 {
00117 //======================================================================
00118 // FILL_IN: [Document your code!!]
00119 //======================================================================
00120   if(cMakeTree==1)
00121     {
00122       TDirectory* current_dir = gDirectory;
00123       fFile->cd();
00124       fTree->SetDirectory(fFile);
00125       fTree->Write();
00126       fFile->Write();
00127       delete fTree;
00128       fFile->Close();
00129       current_dir->cd();
00130     }
00131 
00132 }
00133 
00134 //......................................................................
00135 
00136 void FCPCFilterModule::ResetParams()
00137 {
00138   fRun = -1;
00139   fSnarl = -1;
00140   fMinPlane =  999;
00141   fMaxPlane = -999;
00142   fCharge = 0.0;
00143   fNPlane = 0;
00144   memset(fFidCharge, 0, 8*sizeof(double));
00145   memset(fFidU, 0, 12*sizeof(double));//new
00146   memset(fFidV, 0, 12*sizeof(double));//new
00147   memset(fFidZ, 0, 12*sizeof(double));//new
00148   fSumFidCharge=0;
00149   fSM1begPlaneQ = 0;
00150   fSM1endPlaneQ = 0;
00151   fSM2begPlaneQ = 0;
00152   fSM2endPlaneQ = 0;
00153   fChargeInSM1 = 0;
00154   fChargeInSM2 = 0;
00155   fEventClass = 0;
00156   
00157   for(int i=0; i<2; ++i) {
00158     for(int j=0; j<2; ++j) {
00159       for(int k=0; k<2; ++k) {
00160         for(int l=0; l<2; ++l) {
00161           wPlaneContainment[i][j][k][l]=0;
00162         }
00163       }
00164     }
00165   }
00166   //Cases without an UC event from Planes.
00167   wPlaneContainment[0][0][0][0]=2;//FC -  within a SM
00168   wPlaneContainment[0][0][0][1]=1;//PC - in/out through end of SM2
00169   wPlaneContainment[0][0][1][0]=1;//PC - in/out through beg of SM2
00170   wPlaneContainment[0][1][0][0]=1;//PC - in/out through end of SM1
00171   wPlaneContainment[0][1][1][0]=2;//FC - passing through SM gap
00172   wPlaneContainment[0][1][1][1]=-1;//PC - through SM2, starts/stops in SM1
00173   wPlaneContainment[1][0][0][0]=1;//PC - in/out through beg of SM1
00174   wPlaneContainment[1][1][1][0]=-1;//PC - through SM1, starts/stops in SM2
00175 
00176   memset(wPlaneChargePE, 0, 500*sizeof(double));
00177   memset(wPlaneMeanTPos, 0, 500*sizeof(double));
00178   memset(wPlaneZPos,     0, 500*sizeof(double));
00179   //  memset(wMeanTime, 0, 12*sizeof(double)); //upward going mus
00180 
00181 }
00182 
00183 
00184 JobCResult FCPCFilterModule::Ana(const MomNavigator* mom)
00185 {
00186   MSG("FCPCFilter", Msg::kDebug)<< "FCPCFilterModule::Ana"<<endl;
00187   assert(mom);
00188   JobCResult result(JobCResult::kPassed);
00189 
00190   ResetParams();
00191 
00192   //==============================================================
00193   //Find all the records we need
00194   //
00195 
00196 
00197   CandRecord *candrec = dynamic_cast<CandRecord *>
00198     (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
00199   if (!candrec) {
00200     result.SetFailed();
00201     return result;
00202   }
00203   const CandHeader* ch =   candrec->GetCandHeader();
00204   fRun =  ch->GetRun();
00205   fSnarl = ch->GetSnarl();
00206   CandDigitListHandle *cdlh = dynamic_cast<CandDigitListHandle *>
00207     (candrec->FindCandHandle("CandDigitListHandle", cListIn.c_str()));
00208   if (!cdlh || cdlh->GetNDaughters() < 1) {
00209     MSG("FCPCFilter", Msg::kError)<< "FCPCFilter::Ana Null CandDigit list.  Bail out of event." << endl;
00210     result.SetWarning().SetFailed();
00211     return result;
00212   }
00213 
00214   MSG("FCPCFilter", Msg::kDebug)<< "Name of digit list: "<< cdlh->ClassName()<<cdlh->GetTitle()<<endl;
00215 
00216   //============================================================
00217   //Set up handles 
00218   //
00219   UgliGeomHandle ugh(*cdlh->GetVldContext());
00220   CandDigitHandleItr cdhItr(cdlh->GetDaughterIterator());
00221   
00222 
00223   //============================================================
00224   //Build some fast access local arrays
00225   //
00226   while (CandDigitHandle *cdh = (cdhItr())){
00227     if(cdh->GetPlexSEIdAltL().size()==0)continue;
00228     const PlexSEIdAltLItem psealtli = cdh->GetPlexSEIdAltL().GetBestItem();
00229     const PlexStripEndId  pseid = cdh->GetPlexSEIdAltL().GetBestSEId();    
00230     const int plane = pseid.GetPlane();
00231     if(plane<0 || plane>=500 || !pseid.IsValid())continue; //skip dark channels and shield
00232     UgliStripHandle ush = ugh.GetStripHandle(pseid);
00233     MSG("FCPCFilter", Msg::kDebug)<< "Trying Q1: "<<plane<<endl;
00234     const double  pecharge =cdh->GetCharge(CalDigitType::kPE);
00235     MSG("FCPCFilter", Msg::kDebug)<< "Done Q1: "<<plane<<endl;
00236     const CandDeMuxDigitHandle* demux_digit = dynamic_cast<const CandDeMuxDigitHandle*>(cdh);
00237     if(demux_digit)
00238       {
00239         if((demux_digit->GetDeMuxDigitFlagWord()&CandDeMuxDigit::kXTalk)!=0)continue;
00240       }
00241     const double  tpos = ush.GetTPos();
00242     wPlaneZPos[plane]=ush.GlobalPos(0).Z();
00243     //Do we use this digit?
00244     //
00245     if(psealtli.GetWeight()>cMinDmxWeight &&
00246        pecharge            >cMinDigitCharge )
00247       {
00248         wPlaneChargePE[plane] += pecharge;
00249         wPlaneMeanTPos[plane] += (pecharge*tpos);
00250         if(plane>fMaxPlane) fMaxPlane = plane;
00251         if(plane<fMinPlane) fMinPlane = plane;   
00252         fCharge+= pecharge;
00253       }
00254   }
00255 
00256   //============================================================
00257   //Clean up the work arrays
00258   //
00259   for(int ipln=fMinPlane; ipln<=fMaxPlane; ++ipln)
00260     {
00261       if(wPlaneChargePE[ipln]>0.0)
00262         {
00263           wPlaneMeanTPos[ipln]/=wPlaneChargePE[ipln];
00264         }
00265 
00266 //       if(ipln>0 && ipln<(cSM1StartPlane+cEndPlaneVeto)                 /*>0 && <1+5*/) fSM1begPlaneQ+=wPlaneChargePE[ipln];
00267 //       if(ipln>(cSM1StopPlane-cEndPlaneVeto) && ipln<cSM2StartPlane/* >248-5 && <250*/) fSM1endPlaneQ+=wPlaneChargePE[ipln];
00268 //       if(ipln>cSM1StopPlane && ipln<(cSM2StartPlane+cEndPlaneVeto)/* >248 && <250+5*/) fSM2begPlaneQ+=wPlaneChargePE[ipln];
00269 //       if(ipln>(cSM2StopPlane-cEndPlaneVeto) && ipln<486            /*>485-5 && <486*/) fSM2endPlaneQ+=wPlaneChargePE[ipln];
00270       if(wPlaneChargePE[ipln]>=cMinPlaneCharge) fNPlane++;
00271     }
00272   
00273   //============================================================
00274   //Guess 3d hit postion of each hit
00275   //
00276   cdhItr.Reset();
00277   //  Double_t fAbsTime = cdlh->GetAbsTime()*1.0E9; //upward muons
00278   while (CandDigitHandle *cdh = (cdhItr())) 
00279     {
00280        if(cdh->GetPlexSEIdAltL().size()==0)continue;
00281        
00282        const PlexSEIdAltLItem psealtli = cdh->GetPlexSEIdAltL().GetBestItem();
00283        const PlexSEIdAltL psealt = cdh->GetPlexSEIdAltL();
00284        const PlexStripEndId  pseid = cdh->GetPlexSEIdAltL().GetBestSEId();    
00285        const int plane = pseid.GetPlane();
00286        if(plane<0 || plane>=500 || !pseid.IsValid())continue; //skip dark channels and shield
00287        UgliStripHandle ush = ugh.GetStripHandle(pseid);
00288        MSG("FCPCFilter", Msg::kDebug)<< "Trying Q2: "<<plane<<":"<<psealt.GetDemuxVetoFlag()<<" Veto:"<< psealt.IsVetoShield()<<endl;
00289        const double  pecharge =cdh->GetCharge(CalDigitType::kPE);  
00290        MSG("FCPCFilter", Msg::kDebug)<< "Done Q2: "<<plane<<endl;
00291        const double  tpos = ush.GetTPos();
00292        if(psealtli.GetWeight()<cMinDmxWeight ||
00293           pecharge            <cMinDigitCharge ) continue;
00294        
00295        const CandDeMuxDigitHandle* demux_digit = dynamic_cast<const CandDeMuxDigitHandle*>(cdh);
00296        if(demux_digit)
00297          {
00298            if((demux_digit->GetDeMuxDigitFlagWord()&CandDeMuxDigit::kXTalk)!=0)continue;
00299          }
00300 //        Double_t time =  cdh->GetTime()*1.0E9; //upward muons
00301 //        float adc  = psealtli.GetSigLin(); //upward muons
00302 //        time = time - fAbsTime+10.-TimeWalk(adc); //upward muons
00303 //        float clear = ush.ClearFiber(psealt.GetEnd()); //upward muons
00304 //        float wls = ush.WlsPigtail(psealt.GetEnd()); //upward muons
00305 //        time = time - (clear + wls)*5.9; //upward muons
00306        double nz[2]={-1.0, -1.0};
00307        double ntpos[2]={-1.0, -1.0};
00308        double stp_guess_z=wPlaneZPos[plane];
00309        double stp_guess_tpos=-9999.0;
00310        if(plane>0 && plane<500 && wPlaneChargePE[plane-1]>cMinDigitCharge && wPlaneChargePE[plane+1]>cMinDigitCharge)
00311          {
00312            nz[0] = wPlaneZPos[plane-1];
00313            nz[1] = wPlaneZPos[plane+1];
00314            ntpos[0] = wPlaneMeanTPos[plane-1];
00315            ntpos[1] = wPlaneMeanTPos[plane+1];
00316          }
00317        else if(plane>2  && plane<500 && wPlaneChargePE[plane-1]>cMinDigitCharge && wPlaneChargePE[plane-3]>cMinDigitCharge)
00318          {
00319            nz[0] = wPlaneZPos[plane-1];
00320            nz[1] = wPlaneZPos[plane-3];
00321            ntpos[0] = wPlaneMeanTPos[plane-1];
00322            ntpos[1] = wPlaneMeanTPos[plane-3];
00323          }
00324        else if(plane<497  && plane>0 && wPlaneChargePE[plane+1]>cMinDigitCharge && wPlaneChargePE[plane+3]>cMinDigitCharge)
00325          {
00326            nz[0] = wPlaneZPos[plane+1];
00327            nz[1] = wPlaneZPos[plane+3];
00328            ntpos[0] = wPlaneMeanTPos[plane+1];
00329            ntpos[1] = wPlaneMeanTPos[plane+3];      
00330          }
00331        else if(plane>1 && wPlaneChargePE[plane-1]>cMinDigitCharge)
00332          {stp_guess_tpos = wPlaneMeanTPos[plane-1];}
00333        else if(plane<497 && wPlaneChargePE[plane+1]>cMinDigitCharge)
00334          {stp_guess_tpos = wPlaneMeanTPos[plane+1];}
00335        else{/*lone hit*/}
00336        
00337        
00338        //====================================
00339        //Extrapolate position
00340        //
00341        if(stp_guess_tpos<-1000.0 &&  nz[0]>0.0)
00342          {
00343            stp_guess_tpos =((ntpos[0]-ntpos[1])/(nz[0]-nz[1]))*stp_guess_z + 
00344              ((ntpos[1]*nz[0] - ntpos[0]*nz[1]) / (nz[0]-nz[1]));
00345          }
00346        
00347        //Find min fid side
00348        //
00349        if(stp_guess_tpos>-1000.0)
00350          {
00351 //         float green(0.0); //green fibre correction
00352 //         float dist(0.0);
00353 //         if(pseid.GetPlaneView()==PlaneView::kU) dist -= stp_guess_tpos;
00354 //         if(pseid.GetPlaneView()==PlaneView::kV) dist += stp_guess_tpos;
00355 //         if(psealt.GetEnd()==StripEnd::kNegative) green = ush.GetHalfLength()+dist;//green fibre correction
00356 //         if(psealt.GetEnd()==StripEnd::kPositive) green = ush.GetHalfLength()-dist;//green fibre correction
00357 //         time-=green*5.9; //green fibre correction to hit time
00358            double u = (pseid.GetPlaneView()==PlaneView::kU)?tpos:stp_guess_tpos;
00359            double v = (pseid.GetPlaneView()==PlaneView::kV)?tpos:stp_guess_tpos;
00360            double y = (u+v)*0.707106781;
00361            double x = (u-v)*0.707106781;
00362            if(x> 4.0)  x = 4.0;if(x<-4.0) x = -4.0;
00363            if(y> 4.0)  y = 4.0;if(y<-4.0) y = -4.0;
00364            if(u> 4.0)  u = 4.0;if(u<-4.0) u = -4.0;
00365            if(v> 4.0)  v = 4.0;if(v<-4.0) v = -4.0;
00366            
00367            int minFidSide = 0;
00368            double minFid  = fabs(4.0 - y);
00369            if(fabs( 4.0 - u)<minFid) {minFidSide = 1; minFid = fabs( 4.0 - u);}
00370            if(fabs( 4.0 - x)<minFid) {minFidSide = 2; minFid = fabs( 4.0 - x);}
00371            if(fabs(-4.0 - v)<minFid) {minFidSide = 3; minFid = fabs(-4.0 - v);}
00372            if(fabs(-4.0 - y)<minFid) {minFidSide = 4; minFid = fabs(-4.0 - y);}
00373            if(fabs(-4.0 - u)<minFid) {minFidSide = 5; minFid = fabs(-4.0 - u);} 
00374            if(fabs(-4.0 - x)<minFid) {minFidSide = 6; minFid = fabs(-4.0 - x);}
00375            if(fabs( 4.0 - v)<minFid) {minFidSide = 7; minFid = fabs( 4.0 - v);}
00376            if( minFid<cFidRadius && ( (plane>=(cSM1StartPlane+cEndPlaneVeto) && plane<=(cSM1StopPlane-cEndPlaneVeto) ) || (plane>=(cSM2StartPlane+cEndPlaneVeto) && plane<=(cSM2StopPlane-cEndPlaneVeto) ) ) )
00377              {
00378                fFidCharge[minFidSide] += pecharge;
00379                fSumFidCharge          += pecharge;
00380                //              wMeanTime[minFidSide]+=time*pecharge; //upward going mus
00381                fFidU[minFidSide] += u*pecharge;//use to calc mean hit position in a region
00382                fFidV[minFidSide] += v*pecharge;//use to calc mean hit position in a region
00383                fFidZ[minFidSide] += wPlaneZPos[plane]*pecharge;//use to calc mean hit position in a region
00384              }
00385            else if(plane>0 && plane<(cSM1StartPlane+cEndPlaneVeto)                 /*>0 && <1+5*/)
00386              {
00387                minFidSide=8;
00388                fSM1begPlaneQ     += pecharge; //new 
00389                //              wMeanTime[minFidSide]+=time*pecharge; //upward going mus
00390                fFidU[minFidSide] += u*pecharge;//use to calc mean hit position in a region
00391                fFidV[minFidSide] += v*pecharge;//use to calc mean hit position in a region
00392                fFidZ[minFidSide] += wPlaneZPos[plane]*pecharge;//use to calc mean hit position in a region
00393              }
00394            else if(plane>(cSM1StopPlane-cEndPlaneVeto) && plane<cSM2StartPlane/* >248-5 && <250*/) 
00395              {
00396                minFidSide=9;
00397                fSM1endPlaneQ     += pecharge; //new 
00398                //              wMeanTime[minFidSide]+=time*pecharge; //upward going mus
00399                fFidU[minFidSide] += u*pecharge;//use to calc mean hit position in a region
00400                fFidV[minFidSide] += v*pecharge;//use to calc mean hit position in a region
00401                fFidZ[minFidSide] += wPlaneZPos[plane]*pecharge;//use to calc mean hit position in a region
00402              }
00403            else if(plane>cSM1StopPlane && plane<(cSM2StartPlane+cEndPlaneVeto)/* >248 && <250+5*/) 
00404              {
00405                minFidSide=10;
00406                fSM2begPlaneQ     += pecharge; //new 
00407                //              wMeanTime[minFidSide]+=time*pecharge; //upward going mus
00408                fFidU[minFidSide] += u*pecharge;//use to calc mean hit position in a region
00409                fFidV[minFidSide] += v*pecharge;//use to calc mean hit position in a region
00410                fFidZ[minFidSide] += wPlaneZPos[plane]*pecharge;//use to calc mean hit position in a region
00411              }
00412            else if( plane>(cSM2StopPlane-cEndPlaneVeto) && plane<486            /*>485-5 && <486*/)
00413              {
00414                minFidSide=11;
00415                fSM2endPlaneQ     += pecharge; //new 
00416                //              wMeanTime[minFidSide]+=time*pecharge; //upward going mus
00417                fFidU[minFidSide] += u*pecharge;//use to calc mean hit position in a region
00418                fFidV[minFidSide] += v*pecharge;//use to calc mean hit position in a region
00419                fFidZ[minFidSide] += wPlaneZPos[plane]*pecharge;//use to calc mean hit position in a region
00420              }
00421            else if(minFid>=cFidRadius) // look at PH in each SM
00422              {
00423                if( plane>=(cSM1StartPlane+cEndPlaneVeto) && plane<=(cSM1StopPlane-cEndPlaneVeto) )  fChargeInSM1 += pecharge; //new 
00424                if( plane>=(cSM2StartPlane+cEndPlaneVeto) && plane<=(cSM2StopPlane-cEndPlaneVeto) )  fChargeInSM2 += pecharge; //new 
00425              }
00426          }      
00427      }
00428 
00429    //calculate the mean positions for each region
00430    for(int i=0; i<12; ++i)
00431      {
00432        Double_t regionQ=0.;
00433        if(i<8) regionQ = fFidCharge[i];
00434        if(i==8) regionQ = fSM1begPlaneQ;
00435        if(i==9) regionQ = fSM1endPlaneQ;
00436        if(i==10) regionQ = fSM2begPlaneQ;
00437        if(i==11) regionQ = fSM2endPlaneQ;
00438        if(regionQ>0)
00439          {
00440            wMeanTime[i]/=regionQ; //upward going mus
00441            fFidU[i]/=regionQ;
00442            fFidV[i]/=regionQ;
00443            fFidZ[i]/=regionQ;
00444          }
00445      }
00446    //======================================================
00447    //Determine the event type
00448    //
00449    int chargeside[2] ={-1, -1}; 
00450    int nside = 0;
00451    if(fSumFidCharge<cMinFidCharge) fEventClass = 2; //fSumFidCharge only includes contributions from the sides not the ends
00452    else
00453      {
00454        //deal with the sides of the detector
00455        for(unsigned int iside = 0; iside<=7; ++iside)
00456          {
00457            if(fFidCharge[iside]>cMinFidCharge)
00458              {
00459                if(nside<2)
00460                  {
00461                    chargeside[nside] = iside;
00462                  }
00463                else chargeside[1] =(chargeside[0]!=0 || iside!=7)?iside:chargeside[1]; //bit hacky
00464                nside++;
00465              }
00466          }
00467        if(nside<=2)
00468          {
00469            if(nside==1) fEventClass = 1;
00470            else
00471              {//went through 2 sides - are they adjacent?
00472                const int d = abs(chargeside[1] - chargeside[0]);
00473                if(d==7 || d==1) 
00474                  {//yes - check more carefully
00475                    //check separation of mean positions in each side
00476                    Double_t sep = (fFidU[chargeside[0]]-fFidU[chargeside[1]])*(fFidU[chargeside[0]]-fFidU[chargeside[1]]);
00477                    sep += (fFidV[chargeside[0]]-fFidV[chargeside[1]])*(fFidV[chargeside[0]]-fFidV[chargeside[1]]);
00478                    sep += (fFidZ[chargeside[0]]-fFidZ[chargeside[1]])*(fFidZ[chargeside[0]]-fFidZ[chargeside[1]]);
00479                    sep = pow(sep,0.5);
00480                    if(sep<cMaxRegionSep) fEventClass = 1; //hits are grouped close enough together to be from a single source
00481                    else fEventClass = 0; 
00482                  }
00483                else fEventClass = 0; //no - therefore uncontained
00484                
00485              }
00486          }
00487        else
00488          {//went through more than 2 sides therefore uncontained
00489            fEventClass=0;
00490          }
00491      }       
00492    int fEndEventClass(10); //set to a high number if it isn't going to be used
00493    if( fEventClass!=0 )
00494      {//event hasn't already been classified uncontained so check ends
00495        int ucb1 = (fSM1begPlaneQ>cEndPlaneVetoCharge)?1:0;
00496        int uce1 = (fSM1endPlaneQ>cEndPlaneVetoCharge)?1:0;
00497        int ucb2 = (fSM2begPlaneQ>cEndPlaneVetoCharge)?1:0;
00498        int uce2 = (fSM2endPlaneQ>cEndPlaneVetoCharge)?1:0;
00499        fEndEventClass = wPlaneContainment[ucb1][uce1][ucb2][uce2];
00500        if(fChargeInSM1<cMinPlaneCharge && fChargeInSM2<cMinPlaneCharge) fEndEventClass=0; //new
00501        if(fEndEventClass==2 && uce1==1 && ucb2==1)//new
00502          {//new
00503            if(fChargeInSM1<cMinPlaneCharge || fChargeInSM2<cMinPlaneCharge) fEndEventClass=-2; //new -2 is a special flag 
00504          }//new
00505       if(fEndEventClass ==-1)//new
00506          {//new
00507            if(fChargeInSM1<cMinPlaneCharge && uce2==1) { fEndEventClass=0; chargeside[1]=11; } //new
00508            if(fChargeInSM2<cMinPlaneCharge && ucb1==1) { fEndEventClass=0; chargeside[1]=8;  } //new
00509          }//new
00510        if(abs(fEndEventClass)==1 && fEventClass==1)
00511          {//if it appears PC from both sides and ends check how close the two regions are
00512            int SM = (fFidZ[chargeside[0]]<15.0)?1:2;
00513            int check(-1);
00514            int end(0);
00515            //work out which end if any needs to be checked against the side
00516            if(SM==1) //side hits are in SM1
00517              {
00518                if(fEndEventClass==1)
00519                  {
00520                    if(ucb1==1)
00521                      {
00522                        end=8;
00523                        check=1;
00524                      }
00525                    else if(uce1==1)
00526                      {
00527                        end=9;
00528                        check=1;
00529                      }
00530                    else fEndEventClass=0;
00531                  }
00532                else if(fEndEventClass==-1)
00533                  {
00534                    if(ucb1==1)fEndEventClass=1;
00535                    if(uce2==1){ fEndEventClass=0; chargeside[1]=11; } //new
00536                  }
00537              }
00538            if(SM==2) //side hits are in SM2
00539              {
00540                if(fEndEventClass==1)
00541                  {
00542                    if(ucb2==1)
00543                      {
00544                        end=10;
00545                        check=1;
00546                      }
00547                    else if(uce2==1)
00548                      {
00549                        end=11;
00550                        check=1;
00551                      }
00552                    else fEndEventClass=0;
00553                  }
00554                else if(fEndEventClass==-1)
00555                  {
00556                    if(ucb1==1) { fEndEventClass=0; chargeside[1]=8;  } //new
00557                    if(uce2==1)fEndEventClass=1;
00558                  }
00559              }
00560            if(check==1)
00561              {//compare positions if necessary
00562                Double_t sep = (fFidU[chargeside[0]]-fFidU[end])*(fFidU[chargeside[0]]-fFidU[end]);
00563                sep += (fFidV[chargeside[0]]-fFidV[end])*(fFidV[chargeside[0]]-fFidV[end]);
00564                sep += (fFidZ[chargeside[0]]-fFidZ[end])*(fFidZ[chargeside[0]]-fFidZ[end]);
00565                sep = pow(sep,0.5);
00566                if(sep<cMaxRegionSep) fEndEventClass = 1; //hits are grouped close enough together to be from a single source
00567                else 
00568                  {
00569                    fEndEventClass = 0; 
00570                    chargeside[1]=end;
00571                  }
00572              }
00573         }
00574        if(fEndEventClass==-2) fEndEventClass=1; //new deal with this new class of events
00575        //set the containment to the worst one
00576        fEventClass=(fEndEventClass<fEventClass)? abs(fEndEventClass):fEventClass;
00577      }
00578 
00579 
00580    //======================================================
00581    //Apply the cuts
00582    //
00583   
00584    //Total Plane cut
00585    if(fNPlane<cMinPlaneCut) 
00586      {
00587        MSG("FCPCFilter", Msg::kInfo) <<"FCPCFilter::Ana Event "<<fRun<< ":"<<fSnarl<< " Failed Min Plane Cut:"<< fNPlane<<"<"<<cMinPlaneCut<<endl;
00588        result.SetFailed();
00589      }
00590    if(fEventClass ==0 )
00591      {
00592        MSG("FCPCFilter", Msg::kInfo) <<"FCPCFilter::Ana Passed through regions "<<chargeside[0]<< " & "<<chargeside[1]<< "."<<endl;
00593 //        double dy(-1),dt(1);//upward going mus
00594 //        if(chargeside[1]!=-1 && chargeside[0]!=-1)//upward going mus
00595 //       {//upward going mus
00596 //         dy = fFidV[chargeside[1]]+fFidU[chargeside[1]]-fFidV[chargeside[0]]-fFidU[chargeside[0]]; //don't worry about rt2
00597 //         dt = wMeanTime[chargeside[1]]-wMeanTime[chargeside[0]];//upward going mus
00598 //       }
00599 //        MSG("FCPCFilter", Msg::kInfo) <<"FCPCFilter::Ana dy = "<<dy<< " & dt = "<<dt<<"."<<endl;//upward going mus
00600 //        if(dt!=0.0) //upward going mus
00601 //       {//upward going mus
00602 //         if( (dy/dt)>0)//upward going mus
00603 //           {//upward going mus
00604 //             MSG("FCPCFilter", Msg::kInfo) <<"FCPCFilter::Ana Event "<<fRun<< ":"<<fSnarl<< " Upward-going Muon Event"<<endl;//upward going mus
00605 //             if(!cDoUp) result.SetFailed();//upward going mus
00606 //           }//upward going mus
00607 //         else//upward going mus
00608 //           {   //upward going mus
00609 //             MSG("FCPCFilter", Msg::kInfo) <<"FCPCFilter::Ana Event "<<fRun<< ":"<<fSnarl<< " NON-contained Event"<<endl;
00610 //             if(!cDoUC) result.SetFailed();
00611 //           }//upward going mus
00612 //       }
00613 //        else//upward going mus
00614 //       {       //upward going mus
00615        MSG("FCPCFilter", Msg::kInfo) <<"FCPCFilter::Ana Event "<<fRun<< ":"<<fSnarl<< " NON-contained Event"<<endl;
00616        if(!cDoUC) result.SetFailed();
00617            //    }//upward going mus
00618      }
00619    else if(fEventClass ==2 )
00620      {
00621        MSG("FCPCFilter", Msg::kInfo) <<"FCPCFilter::Ana Event "<<fRun<< ":"<<fSnarl<< " FC Event"<<endl;
00622        if(!cDoFC){
00623          result.SetFailed();
00624        }
00625      }
00626    else if(fEventClass ==1 )
00627      {
00628        MSG("FCPCFilter", Msg::kInfo) <<"FCPCFilter::Ana Event "<<fRun<< ":"<<fSnarl<< " PC Event"<<endl;
00629        if(!cDoPC){
00630          result.SetFailed();
00631        }
00632      }
00633    else
00634      {
00635        MSG("FCPCFilter", Msg::kError) <<"FCPCFilter::Ana Event "<<fRun<< ":"<<fSnarl<< " Unknown Event Class"<<endl;
00636        result.SetFailed().SetFatal();
00637      }
00638    
00639    if(fCharge>0.0 && cMaxFidPctChargeCut >0.0 && cMaxFidPctChargeCut < (fSumFidCharge / fCharge))
00640      {
00641        MSG("FCPCFilter", Msg::kInfo) <<"FCPCFilter::Ana Event "<<fRun<< ":"<<fSnarl<< " Failed, Fraction "<< (fSumFidCharge / fCharge)<< " of charge outside fid volume"<<endl;
00642        result.SetFailed();
00643      }
00644      
00645    if(result.Passed())
00646      {
00647        MSG("FCPCFilter", Msg::kInfo) <<"FCPCFilter::Ana Event "<<fRun<< ":"<<fSnarl<< " Passed all cuts."<<endl;
00648      }
00649 
00650 //    string command;
00651 //    cin>>command;
00652    
00653    
00654    MSG("FCPCFilter", Msg::kDebug)<< "Exit FCPCFilterModule::Ana"<<endl;
00655    if(cMakeTree==1) fTree->Fill();
00656    return result; // kNoDecision, kFailed, etc.
00657 }
00658 
00659 //......................................................................
00660 Float_t FCPCFilterModule::TimeWalk(float adc)
00661 {
00662   
00663   // Time walk correction based on DATA parameterisation
00664   // from Andy Blake (new fit using form suggested by Roy Lee)
00665 
00666   Float_t logQ  = log(adc/2.3);
00667   Float_t logQ2 = logQ*logQ;
00668   Float_t logQ3 = logQ2*logQ;
00669   Float_t timeW = 5.7 - 0.76*logQ-0.038*logQ2+0.00541*logQ3; 
00670   
00671   return  timeW/0.3;
00672   
00673 }
00674 
00675 //......................................................................
00676 
00677 const Registry& FCPCFilterModule::DefaultConfig() const
00678 {
00679 //======================================================================
00680 // Supply the default configuration for the module
00681 //======================================================================
00682   static Registry r; // Default configuration for module
00683 
00684   // Set name of config
00685   std::string name = this->GetName();
00686   name += ".config.default";
00687   r.SetName(name.c_str());
00688 
00689   // Set values in configuration
00690   r.UnLockValues();
00691   r.Set("PC",cDoPC);
00692   r.Set("FC",cDoFC);
00693   r.Set("UC",cDoUC);
00694   r.Set("UpMu",cDoUp);
00695   r.Set("NameListIn","demuxdigitlist");
00696   r.Set("MinDmxWeight",cMinDmxWeight);
00697   r.Set("MinDigitCharge",cMinDigitCharge);
00698   r.Set("MinPlaneCharge",cMinPlaneCharge);
00699   r.Set("FidRadius",cFidRadius);
00700   r.Set("MinPlaneCut",cMinPlaneCut);
00701   r.Set("MinPlaneCharge",cMinPlaneCharge);
00702   r.Set("EndPlaneVeto",cEndPlaneVeto);
00703   r.Set("MinFidCharge",cMinFidCharge);
00704   r.Set("MaxFidPctChargeCut",cMaxFidPctChargeCut);
00705   r.Set("MaxRegionSeparation",cMaxRegionSep);
00706   r.Set("SMEndVetoCharge",cEndPlaneVetoCharge);
00707   r.Set("SM1StartPlane",cSM1StartPlane);
00708   r.Set("SM1StopPlane",cSM1StopPlane);
00709   r.Set("SM2StartPlane",cSM2StartPlane);
00710   r.Set("SM2StopPlane",cSM2StopPlane);
00711   r.Set("MakeTree",cMakeTree);
00712   r.Set("OuputFile",fFileName);
00713   r.LockValues();
00714 
00715 
00716   return r;
00717 }
00718 
00719 //......................................................................
00720 
00721 void FCPCFilterModule::Config(const Registry& r)
00722 {
00723 //======================================================================
00724 // Configure the module given the Registry r
00725 //======================================================================
00726   int    tmpi;
00727   const char* tmps;
00728   double tmpd;
00729 
00730   if (r.Get("PC",tmpi)) { cDoPC = tmpi; }
00731   if (r.Get("FC",tmpi)) { cDoFC = tmpi; }
00732   if (r.Get("UC",tmpi)) { cDoUC = tmpi; }
00733   if (r.Get("UpMu",tmpi)) { cDoUp = tmpi; }
00734   if (r.Get("NameListIn",tmps)) { cListIn = tmps; }
00735   if (r.Get("MinDmxWeight", tmpd)){ cMinDmxWeight = tmpd;}
00736   if (r.Get("MinDigitCharge", tmpd)){ cMinDigitCharge = tmpd;}
00737   if (r.Get("MinPlaneCharge", tmpd)){ cMinPlaneCharge = tmpd;}
00738   if (r.Get("FidRadius", tmpd)){ cFidRadius = tmpd;}
00739   if (r.Get("MinPlaneCut", tmpi)){ cMinPlaneCut = tmpi;}
00740   if (r.Get("MinPlaneCharge", tmpd)){ cMinPlaneCharge = tmpd;}
00741   if (r.Get("EndPlaneVeto", tmpi)){ cEndPlaneVeto = tmpi;}
00742   if (r.Get("MinFidCharge", tmpd)){ cMinFidCharge = tmpd;}
00743   if (r.Get("MaxFidPctChargeCut", tmpd)){ cMaxFidPctChargeCut = tmpd;}
00744   if (r.Get("MaxRegionSeparation",tmpd)){ cMaxRegionSep = tmpd;}
00745   if (r.Get("SMEndVetoCharge", tmpd)){ cEndPlaneVetoCharge = tmpd;}
00746   if (r.Get("SM1StartPlane", tmpi)){ cSM1StartPlane = tmpi;}
00747   if (r.Get("SM1StopPlane", tmpi)){  cSM1StopPlane = tmpi;}
00748   if (r.Get("SM2StartPlane", tmpi)){ cSM2StartPlane = tmpi;}
00749   if (r.Get("SM2StopPlane", tmpi)){  cSM2StopPlane = tmpi;}
00750   if (r.Get("MakeTree", tmpi)){ cMakeTree = tmpi;}
00751   if (r.Get("OuputFile", tmps)){ fFileName = tmps;}
00752   
00753 }
00754