loon: AltDeMuxModule.cxx Source File

00001 
00002 //$Id: AltDeMuxModule.cxx,v 1.12 2006/12/01 20:39:57 rhatcher Exp $
00003 //
00004 //AltDeMuxModule.cxx
00005 //
00006 // An alternative multi-purpose demultiplexing algorithm module for 
00007 // demultiplexing events in the far detector
00008 //
00009 // Author:  M. Thosmson April 2003
00011 
00012 #include "AltDeMuxModule.h"
00013 #include "AltDeMuxDisplay.h"
00014 #include "AltDeMuxCalc.h"
00015 #include "AlgAltDeMux.h"
00016 #include "JobControl/JobCModuleRegistry.h"
00017 #include "JobControl/JobCommand.h"
00018 #include "MessageService/MsgService.h"
00019 #include "MinosObjectMap/MomNavigator.h"
00020 #include "Algorithm/AlgConfig.h"
00021 #include "CandDigit/CandDeMuxDigitListHandle.h"
00022 #include "CandDigit/CandDigitHandle.h"
00023 #include "CandDigit/CandDeMuxDigitHandle.h"
00024 #include "DataUtil/GetCandidate.h"
00025 #include "Algorithm/AlgFactory.h"
00026 #include "Algorithm/AlgHandle.h"
00027 #include "CandData/CandRecord.h"
00028 #include "CandData/CandHeader.h"
00029 #include "Candidate/CandContext.h"
00030 #include "Navigation/NavKey.h"
00031 #include "Navigation/NavSet.h"
00032 #include "Conventions/PlaneView.h"
00033 #include "CandTrackSR/CandTrackSR.h"
00034 #include "Navigation/XxxItr.h"
00035 #include "RecoBase/CandTrackHandle.h"
00036 #include "RecoBase/CandTrackListHandle.h"
00037 #include "TH1.h"
00038 #include "TH2.h"
00039 #include "TClonesArray.h"
00040 #include "RawData/RawRecord.h"
00041 #include "RawData/RawDaqSnarlHeader.h"
00042 #include "Validity/VldContext.h"
00043 #include "UgliGeometry/UgliGeomHandle.h"
00044 #include "TFolder.h"
00045 #include "TObjArray.h"
00046 #include "TROOT.h"
00047 #include <cassert>
00048 #include <cmath>
00049 
00050 ClassImp(AltDeMuxModule)
00051 
00052 CVSID("$Id: AltDeMuxModule.cxx,v 1.12 2006/12/01 20:39:57 rhatcher Exp $");
00053 
00054 // Declare this module to JobControl. Arguments are:
00055 //  (1) The class name 
00056 //  (2) The human-readable name 
00057 //  (3) A short, human-readable description of what the module does
00058 JOBMODULE(AltDeMuxModule, 
00059           "AltDeMuxModule",
00060           "An alternatice module used for demultiplexing far detector events");
00061 
00062 //......................................................................
00063 AltDeMuxModule::AltDeMuxModule() :
00064   fDraw(0),
00065   fDeMuxAlg("AlgAltDeMux"),
00066   fHistoFileName("altDeMux.root")
00067 {
00068  
00069   AlgFactory &af = AlgFactory::GetInstance();
00070 
00071   
00072   af.Register("AlgAltDeMux", "default", "libAltDeMux.so", "AlgConfig");
00073   AlgHandle ah = af.GetAlgHandle(fDeMuxAlg, "default");
00074 
00075   //get the AlgConfig class and set some variables for demuxing
00076   AlgConfig &acd = ah.GetAlgConfig();
00077   acd.UnLockValues();
00078    
00079   acd.Set("XTalkLowPEFraction", 0.1);
00080   acd.Set("XTalkFraction", 0.025);
00081   acd.Set("XTalkLowPECut", 1.5);
00082   
00083   acd.LockValues();
00084   acd.LockKeys();
00085 
00086   if(fDraw==0)fDiagnosticPlots = false;
00087   if(fDraw!=0)fDiagnosticPlots = true;
00088   fDiagnosticPlots = false;
00089 
00090   MSG("AltDeMuxModule", Msg::kDebug) << "AltDeMuxModule::Constructor" << endl;
00091   Float_t fScintillatorN = 1.79;  
00092   Float_t fWLSFibreN     = 1.92; 
00093   Float_t fClearFibreN   = 1.82;  
00094   // velocities of light in m/ns
00095   fClearFibreC   = 0.30/fClearFibreN;
00096   fWLSFibreC     = 0.30/fWLSFibreN;
00097   fScintillatorC = 0.30/fScintillatorN;
00098   //  
00099   fNumberOfStrips = 192;
00100   fDiagnosticPlots = true;
00101 
00102   pCalculator = new AltDeMuxCalc();
00103   pCalculator->SetClearFibreC(fClearFibreC);
00104   pCalculator->SetWLSFibreC(fWLSFibreC);
00105   pCalculator->SetScintillatorC(fScintillatorC);
00106   pCalculator->SetNumberOfStrips(fNumberOfStrips);
00107 
00108   fUZFitSSM1= NULL;
00109   fUZFitSSM2= NULL;
00110   fUZFitTSM1= NULL;
00111   fUZFitTSM2= NULL;
00112   fVZFitSSM1= NULL;
00113   fVZFitSSM2= NULL;
00114   fVZFitTSM1= NULL;
00115   fVZFitTSM2= NULL;
00116 
00117 
00118 
00119   this->BookHistos();
00120   
00121 }
00122 
00123 
00124 void AltDeMuxModule::BookHistos(){
00125 
00126   fQ1Q2Good=new TH2F("fQ1Q2Good","Q1 vs. Q2",200, 0.0, 20.0, 200, 0.0, 20.);
00127   fQ1Q2All=new TH2F("fQ1Q2All",  "Q1 vs. Q2",200, 0.0, 20.0, 200, 0.0, 20.);
00128   fUrat=new TH2F("fUrat","rat vs plane",96, 0., 192., 100, -1.0, 1.0);
00129   fVrat=new TH2F("fVrat","rat vs plane",96, 0., 192., 100, -1.0, 1.0);
00130   fUQdl=new TH2F("fUQdl","qE/qW vs le-lw",100, -10., 10., 100, 0.0, 10.0);
00131   fVQdl=new TH2F("fVQdl","qE/qW vs le-lw",100, -10., 10., 100, 0.0, 10.0);
00132   fDStripVsPlaneU=new TH2F("fDStripVsPlaneU","dstrip vs plane",500,0.0,500.0,       100,-50.,50.);
00133   fDStripVsPlaneV=new TH2F("fDStripVsPlaneV","dstrip vs plane",500,0.0,500.0,
00134     100,-50.,50.);
00135   fUdt=new TH2F("fUdt","dt vs plane",96, 0., 192., 100, -50.0, 50.0);
00136   fVdt=new TH2F("fVdt","dt vs plane",96, 0., 192., 100, -50.0, 50.0);
00137   fUwlsdt=new TH2F("fUwlsdt","dt vs plane",100, -2., 2., 100, -50.0, 50.0);
00138   fVwlsdt=new TH2F("fVwlsdt","dt vs plane",100, -2., 2., 100, -50.0, 50.0);
00139   fUcleardt=new TH2F("fUcleardt","dt vs plane",100, -2., 2., 100, -50.0, 50.0);
00140   fVcleardt=new TH2F("fVcleardt","dt vs plane",100, -2., 2., 100, -50.0, 50.0);
00141   fUwlsNdt=new TH2F("fUwlsNdt","dt vs plane",100, -2., 2., 100, -50.0, 50.0);
00142   fVwlsNdt=new TH2F("fVwlsNdt","dt vs plane",100, -2., 2., 100, -50.0, 50.0);
00143   fUclearNdt=new TH2F("fUclearNdt","Ndt vs plane",100, -2., 2., 100, -50.0, 50.0);
00144   fVclearNdt=new TH2F("fVclearNdt","Ndt vs plane",100, -2., 2., 100, -50.0, 50.0);
00145 
00146   fUG=new TH1F("fUG",  "sigmaQU ",100, -5.0, 5.0);
00147   fVG=new TH1F("fVG",  "sigmaQV ",100, -5.0, 5.0);
00148   fUsdt=new TH1F("fUsdt",  "sigmatV ",100, -50.0, 50.0);
00149   fVsdt=new TH1F("fVsdt",  "sigmatV ",100, -50.0, 50.0);
00150   fUsdtH=new TH1F("fUsdtH",  "sigmatV ",100, -50.0, 50.0);
00151   fVsdtH=new TH1F("fVsdtH",  "sigmatV ",100, -50.0, 50.0);
00152   fUsdq=new TH1F("fUsdq",  "sigma q U",100, -10.0, 10.0);
00153   fVsdq=new TH1F("fVsdq",  "sigma q V",100, -10.0, 10.0);
00154   fVsdqQ=new TH2F("fVsdqQ",  "sigma q V vs Q",20,0.,1000.,100, -10.0, 10.0);
00155   fVsdqV=new TH2F("fVsdqV",  "sigma q V vs V",20,0.,192.,100, -10.0, 10.0);
00156   fUsdq1=new TH1F("fUsdq1",  "sigma q U",100, -10.0, 10.0);
00157   fVsdq1=new TH1F("fVsdq1",  "sigma q V",100, -10.0, 10.0);
00158   fUsdq2=new TH1F("fUsdq2",  "sigma q U",100, -10.0, 10.0);
00159   fVsdq2=new TH1F("fVsdq2",  "sigma q V",100, -10.0, 10.0);
00160 
00161   fUsdq3=new TH1F("fUsdq3",  "sigma q U",100, -10.0, 10.0);
00162   fVsdq3=new TH1F("fVsdq3",  "sigma q V",100, -10.0, 10.0);
00163 
00164   fUsdtc=new TH1F("fUsdtc",  "sigmatUc ",100, -50.0, 50.0);
00165   fVsdtc=new TH1F("fVsdtc",  "sigmatVc ",100, -50.0, 50.0);
00166 
00167 
00168   return;
00169 }
00170 
00171 
00172 bool AltDeMuxModule::SelectCleanMuons()
00173 {
00174 
00175   Float_t NU=0;
00176   Float_t NV=0;
00177   Float_t NUold=0;
00178   Float_t NVold=0;
00179   Float_t xU=0;
00180   Float_t xV=0;
00181   Float_t xxU=0;
00182   Float_t xxV=0;
00183   Float_t yU=0;
00184   Float_t yV=0;
00185   Float_t yUT=0;
00186   Float_t yVT=0;
00187   Float_t xyU=0;
00188   Float_t xyV=0;
00189   Float_t xyUT=0;
00190   Float_t xyVT=0;
00191   Float_t detU;
00192   Float_t detV;
00193   Float_t mU = 0.;
00194   Float_t cU = 0.;
00195   Float_t mV = 0.;
00196   Float_t cV = 0.;
00197   Float_t mUT;
00198   Float_t cUT;
00199   Float_t mVT;
00200   Float_t cVT;
00201   Float_t xp[2];
00202   Float_t yp[2];
00203 
00204   Int_t lowestUplane=999;
00205   Int_t lowestVplane=999;
00206   Int_t highestUplane=-999;
00207   Int_t highestVplane=-999;
00208   cout << " IN ALTDEMUXMODULE::SELECTCLEANMUONS  " << endl;
00209 
00210 
00211   for(UInt_t ipair=0;ipair<fDeMuxedPairs.size();ipair++){
00212     Int_t iplane = fDeMuxedPairs[ipair].altListE->GetPlane();
00213     PlaneView::PlaneView_t kView = fDeMuxedPairs[ipair].altListE->GetPlaneView();
00214     if(kView==PlaneView::kU){
00215       if(iplane<lowestUplane)lowestUplane=iplane;
00216       if(iplane>highestUplane)highestUplane=iplane;
00217     }
00218     if(kView==PlaneView::kV){
00219       if(iplane<lowestVplane)lowestVplane=iplane;
00220       if(iplane>highestVplane)highestVplane=iplane;
00221     }
00222   }
00223   
00224   for(UInt_t ipair=0;ipair<fDeMuxedPairs.size();ipair++){
00225     if(fDeMuxedPairs[ipair].altListE->GetDemuxVetoFlag()==false &&
00226        fDeMuxedPairs[ipair].altListW->GetDemuxVetoFlag()==false){
00227       Int_t iplane = fDeMuxedPairs[ipair].altListE->GetPlane();
00228       PlaneView::PlaneView_t kView = fDeMuxedPairs[ipair].altListE->GetPlaneView();
00229       pCalculator->SetPlane(iplane);
00230       pCalculator->SetView(kView);
00231       Int_t istrip = fDeMuxedPairs[ipair].altListE->GetBestSEId().GetStrip();
00232       pCalculator->SetEast(fDeMuxedPairs[ipair].altListE,istrip);
00233       pCalculator->SetWest(fDeMuxedPairs[ipair].altListW,istrip);
00234       pCalculator->CalcBestEastWest();
00235       Int_t jplane = iplane;
00236       if(jplane>249)jplane+=20;
00237       if(kView==PlaneView::kV){
00238         if(iplane>=lowestUplane-1&&iplane<=highestUplane+1){
00239           NV  += 1.;
00240           xV  += jplane;
00241           xxV += jplane*jplane;
00242           yV  += istrip;
00243           yVT += pCalculator->StripAim();
00244           xyV += jplane*istrip;
00245           xyVT+= jplane*pCalculator->StripAim();
00246         }else{
00247           cout << "Ignoring hits in plane : " << iplane << endl;
00248         }
00249       }
00250       if(kView==PlaneView::kU){
00251         if(iplane>=lowestVplane-1&&iplane<=highestVplane+1){
00252           NU  += 1.;
00253           xU  += jplane;
00254           xxU += jplane*jplane;
00255           yU  += istrip;
00256           yUT += pCalculator->StripAim();
00257           xyU += jplane*istrip;
00258           xyUT+= jplane*pCalculator->StripAim();
00259         }else{
00260           cout << "Ignoring hits in plane : " << iplane << endl;
00261         }
00262       }
00263     }
00264   }
00265   NUold = NU;
00266   NVold = NV;
00267   
00268   detU = xxU*NU-xU*xU;
00269   detV = xxV*NV-xV*xV;
00270   if(detU>0){
00271     mU = (NU*xyU - xU*yU)/detU;
00272     cU = (-xU*xyU + xxU*yU)/detU;
00273     xp[0] = 0;
00274     xp[1] = 249.;
00275     yp[0] = cU;
00276     yp[1] = mU*249.+cU;
00277     if(fUZFitSSM1)delete fUZFitSSM1;
00278     fUZFitSSM1 = new TPolyLine(2,xp,yp);
00279     fUZFitSSM1->SetLineColor(4);
00280     fUZFitSSM1->SetLineStyle(1);
00281     fUZFitSSM1->SetLineWidth(2);
00282     xp[0] = 249;
00283     xp[1] = 500.;
00284     yp[0] = mU*269.+cU;
00285     yp[1] = mU*520.+cU;
00286     if(fUZFitSSM2)delete fUZFitSSM2;
00287     fUZFitSSM2 = new TPolyLine(2,xp,yp);
00288     fUZFitSSM2->SetLineColor(4);
00289     fUZFitSSM2->SetLineStyle(1);
00290     fUZFitSSM2->SetLineWidth(2);
00291 
00292     mVT = (NU*xyUT - xU*yUT)/detU;
00293     cVT = (-xU*xyUT + xxU*yUT)/detU;
00294     xp[0] = 0;
00295     xp[1] = 249.;
00296     yp[0] = cVT;
00297     yp[1] = mVT*249.+cVT;
00298     if(fVZFitTSM1)delete fVZFitTSM1;
00299     fVZFitTSM1 = new TPolyLine(2,xp,yp);
00300     fVZFitTSM1->SetLineColor(4);
00301     fVZFitTSM1->SetLineStyle(2);
00302     fVZFitTSM1->SetLineWidth(2);
00303     xp[0] = 249;
00304     xp[1] = 500.;
00305     yp[0] = mVT*269.+cVT;
00306     yp[1] = mVT*520.+cVT;
00307     if(fVZFitTSM2)delete fVZFitTSM2;
00308     fVZFitTSM2 = new TPolyLine(2,xp,yp);
00309     fVZFitTSM2->SetLineColor(4);
00310     fVZFitTSM2->SetLineStyle(2);
00311     fVZFitTSM2->SetLineWidth(2);
00312   }
00313 
00314   if(detV>0){
00315     mV = (NV*xyV - xV*yV)/detV;
00316     cV = (-xV*xyV + xxV*yV)/detV;
00317     xp[0] = 0;
00318     xp[1] = 249.;
00319     yp[0] = cV;
00320     yp[1] = mV*249.+cV;
00321     if(fVZFitSSM1)delete fVZFitSSM1;
00322     fVZFitSSM1 = new TPolyLine(2,xp,yp);
00323     fVZFitSSM1->SetLineColor(4);
00324     fVZFitSSM1->SetLineStyle(1);
00325     fVZFitSSM1->SetLineWidth(2);
00326     xp[0] = 249;
00327     xp[1] = 500.;
00328     yp[0] = mV*269.+cV;
00329     yp[1] = mV*520.+cV;
00330     if(fVZFitSSM2)delete fVZFitSSM2;
00331     fVZFitSSM2 = new TPolyLine(2,xp,yp);
00332     fVZFitSSM2->SetLineColor(4);
00333     fVZFitSSM2->SetLineStyle(1);
00334     fVZFitSSM2->SetLineWidth(2);
00335 
00336     mUT = (NV*xyVT - xV*yVT)/detV;
00337     cUT = (-xV*xyVT + xxV*yVT)/detV;
00338     xp[0] = 0;
00339     xp[1] = 249.;
00340     yp[0] = cUT;
00341     yp[1] = mUT*249.+cUT;
00342     if(fUZFitTSM1)delete fUZFitTSM1;
00343     fUZFitTSM1 = new TPolyLine(2,xp,yp);
00344     fUZFitTSM1->SetLineColor(4);
00345     fUZFitTSM1->SetLineStyle(2);
00346     fUZFitTSM1->SetLineWidth(2);
00347     xp[0] = 249;
00348     xp[1] = 500.;
00349     yp[0] = mUT*269.+cUT;
00350     yp[1] = mUT*520.+cUT;
00351     if(fUZFitTSM2)delete fUZFitTSM2;
00352     fUZFitTSM2 = new TPolyLine(2,xp,yp);
00353     fUZFitTSM2->SetLineColor(4);
00354     fUZFitTSM2->SetLineStyle(2);
00355     fUZFitTSM2->SetLineWidth(2);
00356   }
00357 
00358 
00359   Float_t chi2U = 0;
00360   Float_t chi2V = 0;
00361   Float_t deltaStrip = 0.;
00362   for(UInt_t ipair=0; ipair<fDeMuxedPairs.size(); ipair++){
00363     if(fDeMuxedPairs[ipair].altListE->GetDemuxVetoFlag()==false &&
00364        fDeMuxedPairs[ipair].altListW->GetDemuxVetoFlag()==false){
00365       Int_t iplane = fDeMuxedPairs[ipair].altListE->GetPlane();
00366       pCalculator->SetPlane(iplane);
00367       PlaneView::PlaneView_t kView = fDeMuxedPairs[ipair].altListE->GetPlaneView();
00368       pCalculator->SetView(kView);
00369       Int_t istrip = fDeMuxedPairs[ipair].altListE->GetBestSEId().GetStrip();
00370       pCalculator->SetEast(fDeMuxedPairs[ipair].altListE,istrip);
00371       pCalculator->SetWest(fDeMuxedPairs[ipair].altListW,istrip);
00372       pCalculator->CalcBestEastWest();
00373       //      Int_t isaim = pCalculator->StripAim();
00374       Int_t jplane = iplane;
00375       if(jplane>249)jplane+=20;
00376       
00377       if(kView==PlaneView::kU){
00378         if(iplane>=lowestUplane-1&&iplane<=highestUplane+1){
00379            deltaStrip = (istrip - mU*jplane-cU);
00380           chi2U += deltaStrip*deltaStrip;
00381         }
00382       }
00383       if(kView==PlaneView::kV){
00384         if(iplane>=lowestUplane-1&&iplane<=highestUplane+1){
00385           deltaStrip = (istrip - mV*jplane-cV);
00386           chi2V += deltaStrip*deltaStrip;
00387         }
00388       }
00389       if(fabs(deltaStrip)>15){
00390         fDeMuxedPairs[ipair].altListE->SetDemuxVetoFlag(true);
00391         fDeMuxedPairs[ipair].altListW->SetDemuxVetoFlag(true);
00392         if(kView==PlaneView::kV){
00393           NV  -= 1.;
00394           xV  -= jplane;
00395           xxV -= jplane*jplane;
00396           yV  -= istrip;
00397           yVT -= pCalculator->StripAim();
00398           xyV -= jplane*istrip;
00399           xyVT-= jplane*pCalculator->StripAim();
00400         }
00401         if(kView==PlaneView::kU){
00402           NU  -= 1.;
00403           xU  -= jplane;
00404           xxU -= jplane*jplane;
00405           yU  -= istrip;
00406           yUT -= pCalculator->StripAim();
00407           xyU -= jplane*istrip;
00408           xyUT-= jplane*pCalculator->StripAim();
00409         }
00410       }
00411     }
00412   }
00413 
00414   detU = xxU*NU-xU*xU;
00415   detV = xxV*NV-xV*xV;
00416   if(detU>0){
00417     mU = (NU*xyU - xU*yU)/detU;
00418     cU = (-xU*xyU + xxU*yU)/detU;
00419     xp[0] = 0;
00420     xp[1] = 249.;
00421     yp[0] = cU;
00422     yp[1] = mU*249.+cU;
00423     if(fUZFitSSM1)delete fUZFitSSM1;
00424     fUZFitSSM1 = new TPolyLine(2,xp,yp);
00425     fUZFitSSM1->SetLineColor(4);
00426     fUZFitSSM1->SetLineStyle(1);
00427     fUZFitSSM1->SetLineWidth(2);
00428     xp[0] = 249;
00429     xp[1] = 500.;
00430     yp[0] = mU*269.+cU;
00431     yp[1] = mU*520.+cU;
00432     if(fUZFitSSM2)delete fUZFitSSM2;
00433     fUZFitSSM2 = new TPolyLine(2,xp,yp);
00434     fUZFitSSM2->SetLineColor(4);
00435     fUZFitSSM2->SetLineStyle(1);
00436     fUZFitSSM2->SetLineWidth(2);
00437 
00438     mVT = (NU*xyUT - xU*yUT)/detU;
00439     cVT = (-xU*xyUT + xxU*yUT)/detU;
00440     xp[0] = 0;
00441     xp[1] = 249.;
00442     yp[0] = cVT;
00443     yp[1] = mVT*249.+cVT;
00444     if(fVZFitTSM1)delete fVZFitTSM1;
00445     fVZFitTSM1 = new TPolyLine(2,xp,yp);
00446     fVZFitTSM1->SetLineColor(4);
00447     fVZFitTSM1->SetLineStyle(2);
00448     fVZFitTSM1->SetLineWidth(2);
00449     xp[0] = 249;
00450     xp[1] = 500.;
00451     yp[0] = mVT*269.+cVT;
00452     yp[1] = mVT*520.+cVT;
00453     if(fVZFitTSM2)delete fVZFitTSM2;
00454     fVZFitTSM2 = new TPolyLine(2,xp,yp);
00455     fVZFitTSM2->SetLineColor(4);
00456     fVZFitTSM2->SetLineStyle(2);
00457     fVZFitTSM2->SetLineWidth(2);
00458   }
00459 
00460   if(detV>0){
00461     mV = (NV*xyV - xV*yV)/detV;
00462     cV = (-xV*xyV + xxV*yV)/detV;
00463     xp[0] = 0;
00464     xp[1] = 249.;
00465     yp[0] = cV;
00466     yp[1] = mV*249.+cV;
00467     if(fVZFitSSM1)delete fVZFitSSM1;
00468     fVZFitSSM1 = new TPolyLine(2,xp,yp);
00469     fVZFitSSM1->SetLineColor(4);
00470     fVZFitSSM1->SetLineStyle(1);
00471     fVZFitSSM1->SetLineWidth(2);
00472     xp[0] = 249;
00473     xp[1] = 500.;
00474     yp[0] = mV*269.+cV;
00475     yp[1] = mV*520.+cV;
00476     if(fVZFitSSM2)delete fVZFitSSM2;
00477     fVZFitSSM2 = new TPolyLine(2,xp,yp);
00478     fVZFitSSM2->SetLineColor(4);
00479     fVZFitSSM2->SetLineStyle(1);
00480     fVZFitSSM2->SetLineWidth(2);
00481 
00482     mUT = (NV*xyVT - xV*yVT)/detV;
00483     cUT = (-xV*xyVT + xxV*yVT)/detV;
00484     xp[0] = 0;
00485     xp[1] = 249.;
00486     yp[0] = cUT;
00487     yp[1] = mUT*249.+cUT;
00488     if(fUZFitTSM1)delete fUZFitTSM1;
00489     fUZFitTSM1 = new TPolyLine(2,xp,yp);
00490     fUZFitTSM1->SetLineColor(4);
00491     fUZFitTSM1->SetLineStyle(2);
00492     fUZFitTSM1->SetLineWidth(2);
00493     xp[0] = 249;
00494     xp[1] = 500.;
00495     yp[0] = mUT*269.+cUT;
00496     yp[1] = mUT*520.+cUT;
00497     if(fUZFitTSM2)delete fUZFitTSM2;
00498     fUZFitTSM2 = new TPolyLine(2,xp,yp);
00499     fUZFitTSM2->SetLineColor(4);
00500     fUZFitTSM2->SetLineStyle(2);
00501     fUZFitTSM2->SetLineWidth(2);
00502   }
00503 
00504   chi2U = 0;
00505   chi2V = 0;
00506   for(UInt_t ipair=0; ipair<fDeMuxedPairs.size(); ipair++){
00507     if(fDeMuxedPairs[ipair].altListE->GetDemuxVetoFlag()==false &&
00508        fDeMuxedPairs[ipair].altListW->GetDemuxVetoFlag()==false){
00509       Int_t iplane = fDeMuxedPairs[ipair].altListE->GetPlane();
00510       pCalculator->SetPlane(iplane);
00511       PlaneView::PlaneView_t kView = fDeMuxedPairs[ipair].altListE->GetPlaneView();
00512       pCalculator->SetView(kView);
00513       Int_t istrip = fDeMuxedPairs[ipair].altListE->GetBestSEId().GetStrip();
00514       pCalculator->SetEast(fDeMuxedPairs[ipair].altListE,istrip);
00515       pCalculator->SetWest(fDeMuxedPairs[ipair].altListW,istrip);
00516       pCalculator->CalcBestEastWest();
00517       //      Int_t isaim = pCalculator->StripAim();
00518       Int_t jplane = iplane;
00519       if(jplane>249)jplane+=20;
00520       
00521       if(kView==PlaneView::kU){
00522         if(iplane>=lowestUplane-1&&iplane<=highestUplane+1){
00523           deltaStrip = (istrip - mU*jplane-cU);
00524           chi2U += deltaStrip*deltaStrip;
00525         }
00526       }
00527       if(kView==PlaneView::kV){
00528         if(iplane>=lowestUplane-1&&iplane<=highestUplane+1){
00529           deltaStrip = (istrip - mV*jplane-cV);
00530           chi2V += deltaStrip*deltaStrip;
00531         }
00532       }
00533       if(deltaStrip>15){
00534         fDeMuxedPairs[ipair].altListE->SetDemuxVetoFlag(true);
00535         fDeMuxedPairs[ipair].altListW->SetDemuxVetoFlag(true);
00536         if(kView==PlaneView::kV){
00537           NV  -= 1.;
00538           xV  -= jplane;
00539           xxV -= jplane*jplane;
00540           yV  -= istrip;
00541           yVT -= pCalculator->StripAim();
00542           xyV -= jplane*istrip;
00543           xyVT-= jplane*pCalculator->StripAim();
00544         }
00545         if(kView==PlaneView::kU){
00546           NU  -= 1.;
00547           xU  -= jplane;
00548           xxU -= jplane*jplane;
00549           yU  -= istrip;
00550           yUT -= pCalculator->StripAim();
00551           xyU -= jplane*istrip;
00552           xyUT-= jplane*pCalculator->StripAim();
00553         }
00554       }
00555     }
00556   }
00557 
00558   cout << "AltDeMuxModule : " << "chi2U " << chi2U << " ndof = " << NV-2 << endl;
00559   cout << "AltDeMuxModule : " << "chi2V " << chi2V << " ndof = " << NU-2 << endl;
00560   bool goodMuon = false;
00561   if(NV>10&&NU>10&&chi2U<10.0*NV&&chi2V<10.0*NV)goodMuon = true;
00562   if(!goodMuon)return false;
00563 
00564   cout << "AltDeMuxModule : " << "Selected Good Muon" << endl;
00565 
00566   for(UInt_t ipair=0; ipair<fDeMuxedPairs.size(); ipair++){
00567     if(fDeMuxedPairs[ipair].altListE->GetDemuxVetoFlag()==false &&
00568        fDeMuxedPairs[ipair].altListW->GetDemuxVetoFlag()==false){
00569       bool onTrack = true;
00570       Int_t iplane = fDeMuxedPairs[ipair].altListE->GetPlane();
00571       pCalculator->SetPlane(iplane);
00572       PlaneView::PlaneView_t kView = fDeMuxedPairs[ipair].altListE->GetPlaneView();
00573       pCalculator->SetView(kView);
00574       Int_t istrip = fDeMuxedPairs[ipair].altListE->GetBestSEId().GetStrip();
00575       pCalculator->SetEast(fDeMuxedPairs[ipair].altListE,istrip);
00576       pCalculator->SetWest(fDeMuxedPairs[ipair].altListW,istrip);
00577       pCalculator->CalcBestEastWest();
00578       Int_t isaim = pCalculator->StripAim();
00579       Int_t jplane = iplane;
00580       if(jplane>249)jplane+=20;
00581             
00582       if(kView==PlaneView::kU){
00583         if(iplane>=lowestUplane-1&&iplane<=highestUplane+1){
00584           deltaStrip = (istrip - mU*jplane-cU);
00585           if(deltaStrip>2.0)onTrack=false;
00586         }
00587       }
00588       if(kView==PlaneView::kV){
00589         if(iplane>=lowestUplane-1&&iplane<=highestUplane+1){
00590           deltaStrip = (istrip - mV*jplane-cV);
00591           if(deltaStrip>2.0)onTrack=false;
00592         }
00593       }
00594       if(onTrack){
00595         Float_t wlsE = pCalculator->GetWlsE();
00596         Float_t wlsW = pCalculator->GetWlsW();
00597         Float_t clearE = pCalculator->GetClearE();
00598         Float_t clearW = pCalculator->GetClearW();
00599         pCalculator->CalcDt();
00600         Float_t dt   = pCalculator->DT();
00601         Double_t tE  = pCalculator->GetTE();
00602         Double_t tW  = pCalculator->GetTW();
00603 
00604         //      cout << "MOD : " << iplane << " : " << tE << " : " << tW << endl;
00605         //      cout << "MOD : " << clearE << "," << clearW << ":" <<  (clearW-clearE)/fWLSFibreC << endl;
00606 
00607         Float_t target = 0.;
00608         //      Float_t mydt = (tW-tE)-(clearW-clearE)/fClearFibreC-(wlsW-wlsE)/fWLSFibreC;
00609         //      cout << dt << ":" << ":" << mydt << endl;
00610 
00611         if(kView==PlaneView::kU)target = mV*jplane+cV;
00612         if(kView==PlaneView::kV)target = mU*jplane+cU;
00613         Float_t yclear = 2.0*fScintillatorN*(target-96.0)*0.04166-
00614           0.3*(tW-tE) + fWLSFibreN*(wlsW-wlsE);
00615         Float_t dtclear = (tW-tE) - (wlsW-wlsE)/fWLSFibreC - 
00616           2.0*(target-96.0)*0.04166/fScintillatorC
00617           - (clearW-clearE)/fClearFibreC;
00618         Float_t ywls = 2.0*fScintillatorN*(target-96.0)*0.04166-
00619           0.3*(tW-tE) + fClearFibreN*(clearW-clearE);
00620 
00621         Float_t dtwls;
00622         if(kView==PlaneView::kU){
00623           dtwls= (tW-tE) - (clearW-clearE)/fClearFibreC - 
00624             2.0*(target-96.0)*0.04166/fScintillatorC
00625             - (wlsW-wlsE)/fWLSFibreC;
00626         }else{
00627           dtwls= (tW-tE) - (clearW-clearE)/fClearFibreC + 
00628             2.0*(target-96.0)*0.04166/fScintillatorC
00629             - (wlsW-wlsE)/fWLSFibreC;
00630         }
00631 
00632 
00633         if(kView==PlaneView::kU){
00634           fDStripVsPlaneU->Fill(iplane,isaim-target,1.);
00635           fUdt->Fill(target,dt,1.);
00636           fUsdt->Fill(target-isaim,1.);
00637           fUsdtc->Fill(target-isaim,1.);
00638 
00639           if(iplane>250)fUsdtH->Fill(target-isaim,1.);
00640           fUclearNdt->Fill(clearW-clearE,yclear,1.);
00641           fUcleardt->Fill(clearW-clearE,dtclear,1.);
00642           fUwlsNdt->Fill(wlsW-wlsE,ywls,1.);
00643           fUwlsdt->Fill(wlsW-wlsE,dtwls,1.);
00644         }
00645         if(kView==PlaneView::kV){
00646           fVsdt->Fill(target-isaim,1.);
00647           if(iplane>250)fVsdtH->Fill(target-isaim,1.);
00648           fDStripVsPlaneV->Fill(iplane,isaim-target,1.);
00649           fVdt->Fill(target,dt,1.);
00650           fVclearNdt->Fill(clearW-clearE,yclear,1.);
00651           fVcleardt->Fill(clearW-clearE,dtclear,1.);
00652           fVwlsNdt->Fill(wlsW-wlsE,ywls,1.);
00653           fVwlsdt->Fill(wlsW-wlsE,dtwls,1.);
00654         }
00655       }
00656     }
00657   }
00658 
00659   return true;
00660 
00661 }
00662 
00663 //----------------------------------------------------------------------
00664 AltDeMuxModule::~AltDeMuxModule()
00665 {
00666   MSG("JobC", Msg::kDebug) << "AltDeMuxModule::Destructor" << endl;
00667   //  if(fDeMuxFile) delete fDeMuxFile;
00668 }
00669 
00670 //......................................................................
00671 
00672 void AltDeMuxModule::BeginJob()
00673 {  // Set Up Canvases
00674 
00675   MSG("JobC", Msg::kDebug) << "AltDeMuxModule::Begin" << endl;
00676  
00677   // if requested draw canvas - useful visual debugging tool
00678 
00679   //if(fDraw)fAltDeMuxDisplay = new AltDeMuxDisplay();
00680 
00681   return;
00682 }
00683 
00684 //......................................................................
00685 
00686 JobCResult AltDeMuxModule::Reco(MomNavigator *mom)
00687 {
00688 
00689   MSG("AltDeMuxModule", Msg::kDebug) << "AltDeMuxModule::Reco" << endl;
00690 
00691   JobCResult result(JobCResult::kPassed);
00692 
00693   // Check that mom exists.
00694   assert(mom);
00695 
00696   // Find PrimaryCandidateRecord fragment in MOM.
00697   CandRecord *candrec = dynamic_cast<CandRecord *>
00698     (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
00699   if (candrec == 0) {
00700     result.SetError().SetFailed();
00701     return result;
00702   }
00703 
00704   //  Find Raw Data CandDeMuxDigitList fragment in PrimaryCandidateRecord.
00705   MSG("AltDeMuxModule", Msg::kDebug) << "got cand record" << endl;
00706   //  Find Raw Data CandDeMuxDigitList fragment in PrimaryCandidateRecord.
00707 
00708   //  CandDeMuxDigitListHandle *crdlh = dynamic_cast<CandDeMuxDigitListHandle *>
00709   //   (candrec->FindCandHandle("CandDeMuxDigitListHandle", "canddigitlist"));
00710   CandDigitListHandle *crdlh = dynamic_cast<CandDigitListHandle *>
00711     (candrec->FindCandHandle("CandDigitListHandle", "canddigitlist"));
00712   if (crdlh == 0) {
00713     MSG("AltDeMuxModule", Msg::kDebug) << "crdlh = 0" << endl;
00714     result.SetError().SetFailed();
00715     return result;
00716    }
00717     
00718 
00719   cout << " ASKING FACTORY FOR : " << fDeMuxAlg << endl;
00720 
00721   AlgFactory &af = AlgFactory::GetInstance();
00722   AlgHandle ah = af.GetAlgHandle(fDeMuxAlg, "default");
00723     
00724   CandContext cx(this, mom);
00725   //cx.SetDataIn(status);
00726   
00727   MSG("AltDeMuxModule", Msg::kDebug) << "starting demux algorithm " << endl;
00728   ah.RunAlg(*crdlh, cx);
00729   MSG("AltDeMuxModule", Msg::kDebug) << "finished demux algorithm " << endl;
00730   
00731   
00732   MSG("AltDeMuxModule", Msg::kDebug) << "Event DeMuxed" << endl;
00733   
00734   return result;
00735 }
00736 
00737 //......................................................................
00738 //this method fills various histograms used to test the demuxing algorithm
00739 //efficiency
00740 
00741 JobCResult AltDeMuxModule::Ana(const MomNavigator *mom)
00742 {
00743   MSG("JobC", Msg::kDebug) << "AltDeMuxModule::Ana" << endl;
00744 
00745   // if the event displat has been asked for - fraw it.
00746   //if(fDiagnosticPlots)fAltDeMuxDisplay->Ana(mom);
00747 
00748   JobCResult result(JobCResult::kPassed);
00749   // Find PrimaryCandidateRecord fragment in MOM.
00750   CandRecord* candrecord = dynamic_cast<CandRecord*>(mom->GetFragment("CandRecord","PrimaryCandidateRecord"));
00751   
00752   if (candrecord == 0) {
00753     MSG("AltDeMuxModule", Msg::kWarning) << " AltDeMuxModule::Ana No PrimaryCandidateRecord in MOM." << endl;
00754     return JobCResult::kError;
00755   }
00756 
00757 
00758 
00759 
00760   fDeMuxedPairs.erase(fDeMuxedPairs.begin(),fDeMuxedPairs.end());
00761   for(UInt_t i=0; i<500;i++){
00762     fPlanesAltLists[i][ALG_EAST].erase(fPlanesAltLists[i][ALG_EAST].begin(),fPlanesAltLists[i][ALG_EAST].end());
00763     fPlanesAltLists[i][ALG_WEST].erase(fPlanesAltLists[i][ALG_WEST].begin(),fPlanesAltLists[i][ALG_WEST].end());
00764   }
00765 
00766   
00767   // Find PrimaryCandidateRecord fragment in MOM.
00768   const CandDeMuxDigitListHandle *canddigit = 
00769     DataUtil::GetCandidate<CandDeMuxDigitListHandle>(mom,
00770                                                      "CandDeMuxDigitListHandle", 
00771                                                      "altdemux");
00772   if (canddigit == 0) {
00773     MSG("MyAna",Msg::kWarning) << "Failed to get CandDeMuxDigitListHandle!\n";
00774     return JobCResult::kFailed;
00775   }
00776 
00777 
00778   // get handle to UgliGeometry
00779   const VldContext* vldc = candrecord->GetVldContext();
00780   UgliGeomHandle ugh(*vldc); 
00781   pUgh = &ugh;
00782   pCalculator->SetUgli(pUgh);
00783 
00784   CandDeMuxDigitHandleItr cdhItr(canddigit->GetDaughterIterator());
00785   
00786   PlexSEIdAltL* paltlist;
00787 
00788   Int_t ndigits = 0;
00789   //  Double_t fAbsTime = canddigit->GetAbsTime()*1.0E9;
00790 
00791   while ( CandDeMuxDigitHandle *cdh = cdhItr() ) {
00792     // for each digit
00793     //    Double_t rawTime = cdh->GetTime(CalTimeType::kNone)*1.0E9;
00794     // Double_t corTime = cdh->GetTime(CalTimeType::kT0)*1.0E9;
00795     paltlist = &(cdh->GetPlexSEIdAltLWritable());
00796     Float_t plane = paltlist->GetPlane();
00797     //    cout << plane << " : " << rawTime << " : " << corTime << " : " << corTime-rawTime << endl;;
00798     if(!paltlist->IsVetoShield()){
00799       if(cdh->GetDeMuxDigitFlagWord()==0){  
00800         ndigits++;
00801         if(paltlist->GetBestWeight()>0.99){
00802           Int_t  iplane = static_cast<int>(plane);
00803           if(cdh->GetDeMuxDigitFlagWord()==false){
00804             if(paltlist->GetEnd()==StripEnd::kEast){
00805               fPlanesAltLists[iplane][ALG_EAST].push_back(paltlist);
00806             }
00807             if(paltlist->GetEnd()==StripEnd::kWest){
00808               fPlanesAltLists[iplane][ALG_WEST].push_back(paltlist);
00809             }
00810           }
00811         }
00812       }
00813     }
00814   }
00815   
00816   cout << "GetFibreLengths" << endl;
00817   this->GetFibreLengths();
00818   cout << "GetFibreLengths DONE" << endl;
00819 
00820   cout << "ALTDEMUXMODULE NDIGITS : " << ndigits << endl;
00821 
00822   Int_t npairs = 0;
00823   Int_t   fUniqueDeMuxedGroupID = 0;
00824   vector <PlexSEIdAltL*>::iterator literE;
00825   vector <PlexSEIdAltL*>::iterator literW;
00826   for(Int_t iplane = 0; iplane <MAX_NUMBER_OF_PLANES; iplane++){
00827     if(fPlanesAltLists[iplane][ALG_EAST].size()&&fPlanesAltLists[iplane][ALG_WEST].size()){
00828       literE = fPlanesAltLists[iplane][ALG_EAST].begin();
00829       while (literE != fPlanesAltLists[iplane][ALG_EAST].end()){
00830         Int_t iStripE = (*literE)->GetBestSEId().GetStrip();
00831         literW = fPlanesAltLists[iplane][ALG_WEST].begin();
00832         while ( literW != fPlanesAltLists[iplane][ALG_WEST].end()){
00833           Int_t iStripW = (*literW)->GetBestSEId().GetStrip();
00834           if(iStripE==iStripW){
00835             npairs++;
00836             DeMuxedPair thisPair;
00837             thisPair.status   = true;
00838             thisPair.altListE = *literE;
00839             thisPair.altListW = *literW;
00840             PlaneView::PlaneView_t kView = (*literE)->GetPlaneView();
00841             pCalculator->SetPlane(iplane);
00842             pCalculator->SetView(kView);
00843             pCalculator->SetEast(*literE,iStripE);
00844             pCalculator->SetWest(*literW,iStripW);
00845             pCalculator->CalcBestEastWest();
00846             if(fabs(pCalculator->SigmaDQ())<3.0){
00847               thisPair.orthogonalStripFromTiming = pCalculator->StripAim();
00848               thisPair.uniqueGroupID = fUniqueDeMuxedGroupID;
00849               thisPair.weightQ = 1.0;
00850               thisPair.pairQCor = pCalculator->PairQCor();
00851               fDeMuxedPairs.push_back(thisPair);
00852               fUniqueDeMuxedGroupID++;
00853             }
00854           }
00855           literW++;
00856         }
00857         literE++;
00858       }
00859     }
00860   }
00861   cout << "ALTDEMUXMODULE NPAIRS : " << npairs << endl;
00862 
00863   this->SelectCleanMuons();
00864 
00865   //this->TimeAnalysis();
00866 
00867   return result;
00868 }
00869 
00870 //......................................................................
00871 
00872 void AltDeMuxModule::HandleCommand(JobCommand *command) 
00873 {
00874   //
00875   //  Purpose:  Method to interpret module commands.
00876   //
00877   //  Arguments:
00878   //    command   in    Command to interpret.
00879   //
00880   //  Return:   n/a
00881   //
00882   //  Commands implemented:
00883   //    Set TestToMake  Set the test to make to obtain a figure of merit - fraction of
00884   //                    digits within n strips from a track, fraction of mated signal
00885   //                    DigitsAlongTrack, and MatedSignalFraction respectively
00886   //    Set DataType    Set the type of data to test - MonteCarlo or FarDetector
00887   //    Set MatedSignalLimit  Set the lower limit for fraction of signal coming from opposite
00888   //                          sides of the same strip in a plane.
00889   //    Set DeltaStripsFromTrack  Set the number of strips a digit can be off from a track and
00890   //                              still be considered along that track.
00891   //
00892   //  see http://beaker.astro.indiana.edu/brebel/demux_notes/how_to_demux.html for details
00893 
00894    MSG("AltDeMux", Msg::kDebug) << "AltDeMuxModule::HandleCommand" << endl;
00895 
00896    TString cmd = command->PopCmd();
00897 
00898    if(cmd=="Set"){
00899      TString opt = command->PopOpt();
00900    }
00901 
00902    return;
00903 }
00904 
00905 
00906 
00907 //......................................................................
00908 
00909 const Registry& AltDeMuxModule::DefaultConfig() const
00910 {
00911   //  Commands implemented:
00912   //    Set TestToMake  Set the test to make to obtain a figure of merit - fraction of
00913   //                    digits within n strips from a track, fraction of mated signal
00914   //                    DigitsAlongTrack, and MatedSignalFraction respectively
00915   //    Set DataType    Set the type of data to test - MonteCarlo or FarDetector
00916   //    Set MatedSignalLimit  Set the lower limit for fraction of signal coming from opposite
00917   //                          sides of the same strip in a plane.
00918   //    Set DeltaStripsFromTrack  Set the number of strips a digit can be off from a track and
00919   //                              still be considered along that track.
00920   //
00921   //  see http://beaker.astro.indiana.edu/brebel/demux_notes/how_to_demux.html for details
00922 
00923   //   int itrue = 1;  // work around for lack of bool in registry
00924   // int ifalse = 0; // work around for lack of bool in registry
00925 
00926    MSG("JobC", Msg::kDebug) << "AltDeMuxModule::DefaultConfig" << endl;
00927 
00928    static Registry r;
00929  
00930    r.UnLockValues();
00931    
00932    r.Set("Draw","fDraw");
00933    r.Set("DeMuxAlg","fDeMuxAlg");
00934 
00935    r.LockValues();
00936       
00937    return r;
00938 }
00939 
00940 //......................................................................
00941 
00942 
00943 void AltDeMuxModule::GetFibreLengths()
00944 {
00945   // Store the fibre lengths for strips in this event.  
00946 
00947   UgliStripHandle ush;
00948 
00949   vector <PlexSEIdAltL*>::iterator literA;
00950   MSG("AlgAltDeMux", Msg::kVerbose) << "AlgAltDeMuxBase::GetFibreLengths " << endl;
00951  
00952   for(Int_t iplane=0;iplane<MAX_NUMBER_OF_PLANES;iplane++){
00953     for(Int_t iew = ALG_EAST; iew <= ALG_WEST; iew++){
00954       if(fPlanesAltLists[iplane][iew].size()>0){
00955         literA = fPlanesAltLists[iplane][iew].begin();
00956         while (literA != fPlanesAltLists[iplane][iew].end()){
00957           (*literA)->SetFirst();
00958           while( (*literA)->IsValid() ){
00959             ush = pUgh->GetStripHandle((*literA)->GetCurrentSEId());
00960             if(!ush.IsValid())MSG("AlgAltDeMux", Msg::kFatal) << "AlgAltDeMuxBase::GetFibreLengths => UgliStripHandle NOT VALID : plane " << iplane << " (Carry on regardless)" <<  endl;
00961             if(iew==ALG_EAST)pCalculator->SetFibreLengthE(iplane,*literA);
00962             if(iew==ALG_WEST)pCalculator->SetFibreLengthW(iplane,*literA);
00963             (*literA)->Next();
00964           }
00965           *literA++;
00966         }
00967       }
00968     }
00969   }
00970   return;
00971 
00972 }
00973 
00974 
00975 void AltDeMuxModule::Config(const Registry& r)
00976 {
00977   Int_t tmpint; 
00978   //Double_t tmpdub; 
00979   const char* tmpchar = 0;
00980 
00981   if(r.Get("Draw", tmpint)){    fDraw = tmpint;
00982     MSG("AltDeMux", Msg::kDebug) << "AltDeMuxModule::Config\n" << 
00983       "Drawing on/off " << fDraw << endl;
00984   }
00985 
00986   if(r.Get("DeMuxAlg", tmpchar)){
00987     fDeMuxAlg = tmpchar;
00988     MSG("AltDeMux", Msg::kDebug) << "AltDeMuxModule::Config\n" << 
00989       "DeMuxing Algorithm : " << fDeMuxAlg << endl;
00990   }
00991 
00992   return;
00993 }
00994 
00995 
00996 
00997 
00998 //----------------------------------------------------------------------------------
00999 void AltDeMuxModule::Help() 
01000 {
01001   MSG("JobC", Msg::kDebug) 
01002     << "AltDeMuxModule::Help\n"
01003     <<"AltDeMuxModule is an alternatice demultiplexer events "
01004     <<"for far detector events." << endl
01005     << "THIS IS A VERY PRELIMINARY VERSION"
01006     << endl;
01007 }
01008 
01009 void AltDeMuxModule::EndJob() 
01010 {
01011 
01012   TFile *hfile = new TFile("AltDeMuxHistos.root","recreate");
01013  
01014   fQ1Q2All->TH2F::Write();
01015   fQ1Q2Good->TH2F::Write();
01016   fUrat->TH2F::Write();
01017   fVrat->TH2F::Write();
01018   fUdt->TH2F::Write();
01019   fVdt->TH2F::Write();
01020   fUQdl->TH2F::Write();
01021   fVQdl->TH2F::Write();
01022   fUwlsdt->TH2F::Write();
01023   fVwlsdt->TH2F::Write();
01024   fUcleardt->TH2F::Write();
01025   fVcleardt->TH2F::Write();
01026   fUwlsNdt->TH2F::Write();
01027   fVwlsNdt->TH2F::Write();
01028   fUclearNdt->TH2F::Write();
01029   fVclearNdt->TH2F::Write();
01030   fUG->TH1F::Write();
01031   fVG->TH1F::Write();
01032   fUsdt->TH1F::Write();
01033   fVsdt->TH1F::Write();
01034   fUsdtH->TH1F::Write();
01035   fVsdtH->TH1F::Write();
01036   fUsdq->TH1F::Write();
01037   fVsdq->TH1F::Write();
01038   fVsdqQ->TH2F::Write();
01039   fVsdqV->TH2F::Write();
01040   fUsdq1->TH1F::Write();
01041   fVsdq1->TH1F::Write();
01042   fUsdq2->TH1F::Write();
01043   fVsdq2->TH1F::Write();
01044   fUsdq3->TH1F::Write();
01045   fVsdq3->TH1F::Write();
01046   fVsdtc->TH1F::Write();
01047   fUsdtc->TH1F::Write();
01048   fDStripVsPlaneU->TH2F::Write();
01049   fDStripVsPlaneV->TH2F::Write();
01050 
01051   hfile->Close();
01052   delete hfile;
01053 
01054   return;
01055 }
01056