MuCalFitterModule.cxx

Go to the documentation of this file.

// $Id: MuCalFitterModule.cxx,v 1.12 2006/12/01 19:15:58 rhatcher Exp $
//
// Histograms and fits muon energy spectra for each strip end.
// Also returns truncated mean 
//
// cbs@hep.ucl.ac.uk
#include "MuCalFitterModule.h"
#include "SpectraFit.C"
#include "MuCalTrunc.h"
#include <cstdio>
#include <fstream>
#include <cmath>
// MINOS includes
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "CandData/CandRecord.h"
#include "CandDigit/CandDigitListHandle.h"
#include "CandDigit/CandDigitHandle.h"
#include "RecoBase/CandClusterHandle.h"
#include "RecoBase/CandShowerHandle.h"
#include "CandTrackSR/CandTrackSRHandle.h"
#include "CandTrackSR/CandTrackSRListHandle.h"
#include "RawData/RawDigit.h"
#include "RawData/RawHeader.h"
#include "RawData/RawRecord.h"
#include "RawData/RawChannelId.h"
#include "RawData/RawDaqSnarlHeader.h"
#include "RawData/RawDaqHeaderBlock.h"
#include "RawData/RawDigitDataBlock.h"
#include "RecoBase/CandShowerListHandle.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "UgliGeometry/UgliStripHandle.h"
#include "Validity/VldContext.h"
#include "Conventions/Munits.h"
#include "JobControl/JobCommand.h"
#include "JobControl/JobCModuleRegistry.h"
#include "Calibrator/CalADCtoPE.h"
#include "DatabaseInterface/DbiResultPtr.h"
#include "Calibrator/Calibrator.h"

CVSID("$Id: MuCalFitterModule.cxx,v 1.12 2006/12/01 19:15:58 rhatcher Exp $");
JOBMODULE(MuCalFitterModule,"MuCalFitterModule","Histogram and Fit Muon Energy Spectra for each Strip End");

//......................................................................

MuCalFitterModule::MuCalFitterModule() :
 
fWriteHistos(true),fwhichfit(0)

{
  
  ffile = new TFile("muonfits.root","RECREATE");
  
  char histname[256];
  for(int i=0;i<150;i++){
    for(int j=0;j<192;j++){
      for(int k=0;k<2;k++){
        sprintf(histname,"muhist%i_%i_%i",i,j,k);
        fHistoArray[i][j][k] = new TH1F(histname,histname,120,-0.5,19.5);
        fgain[i][j][k]=-1;
        fspewidth[i][j][k]=-1;
      }
    }
  }  


  fhistoraw0 = new TH1F("historaw0","historaw0",1400,-0.5,139.5);
  fhistoraw1 = new TH1F("historaw1","historaw1",1400,-0.5,139.5);
  fhistoang0 = new TH1F("histoang0","histoang0",1400,-0.5,139.5);
  fhistoang1 = new TH1F("histoang1","histoang1",1400,-0.5,139.5);
  fhistoatt0 = new TH1F("histoatt0","histoatt0",1400,-0.5,139.5);  
  fhistoatt1 = new TH1F("histoatt1","histoatt1",1400,-0.5,139.5);
  
}

//......................................................................

void MuCalFitterModule::HandleCommand(JobCommand* cmd)
{
//======================================================================
// Handle a command send from the job controller
// 
// Inputs: cmd - The parsed job command
//======================================================================
  if (cmd->HaveCmd()) {                      // If we have a command...
    string sc = cmd->PopCmd();               // Get the command
    if (sc == "WriteHistos") {
      fWriteHistos = true;
      return;
    }
    
    else if (sc == "NoFit") {
      fwhichfit = 0;
      return;
    }
    
    else if (sc == "SimpleFit") {
      fwhichfit = 1;
      return;
    }

    else if (sc == "ComplexFit") {
      fwhichfit = 2;
      return;
    }

    // Don't understand command
    MSG("Demo",Msg::kWarning) << "Don't understand '" << sc.c_str() << "'\n";
  }
}

//......................................................................

JobCResult MuCalFitterModule::Ana(const MomNavigator *mom) 
{

  JobCResult result(JobCResult::kPassed);

  int vetoplane[500];
  for(int i=0;i<500;i++) vetoplane[i]=0;

  CandRecord* candrec = dynamic_cast<CandRecord*>
    (mom->GetFragment("CandRecord", "PrimaryCandidateRecord"));
  if (candrec==0) {
    result.SetWarning().SetFailed();
    return result;
  }

  CandTrackListHandle *tracklist = dynamic_cast<CandTrackListHandle*>
    (candrec->FindCandHandle("CandTrackListHandle"));
  if (!tracklist) {
    return result;
  }

  CandShowerListHandle *showerlist = dynamic_cast<CandShowerListHandle*>
    (candrec->FindCandHandle("CandShowerListHandle"));

  //unused: CandDigitListHandle *digitlist = dynamic_cast<CandDigitListHandle*>
  //unused:   (candrec->FindCandHandle("CandDigitListHandle"));



  RawRecord *rr = dynamic_cast<RawRecord *>(mom->GetFragment("RawRecord"));
  if (rr == 0) {
    return result;
  }

  VldContext vc = rr->GetRawHeader()->GetVldContext();

  const RawDigitDataBlock *rddb = dynamic_cast<const RawDigitDataBlock *>
    (rr->FindRawBlock("RawDigitDataBlock"));
  if (rddb == 0) {
    return result;
  }
  

  if(showerlist!=0){

    TIter showerItr(showerlist->GetDaughterIterator());

    while (CandShowerHandle *shower = dynamic_cast<CandShowerHandle*>(showerItr())) {
      CandShowerHandle *showersr = 0;
      if(shower->InheritsFrom("CandShowerHandle")) {
        showersr = dynamic_cast<CandShowerHandle*>(shower);
      }

      for(int i=0;i<=shower->GetLastCluster();i++){
        const CandClusterHandle *cluster = shower->GetCluster(i);
        for(int j=cluster->GetBegPlane();j<=cluster->GetEndPlane();j++) vetoplane[j] = 1;
      }    
    }
  }

  
  TIter trackItr(tracklist->GetDaughterIterator());
  
  while (CandTrackHandle *track = dynamic_cast<CandTrackHandle*>(trackItr())) {
    CandTrackSRHandle *tracksr = 0;
    if (track->InheritsFrom("CandTrackSRHandle")) {
      tracksr = dynamic_cast<CandTrackSRHandle*>(track);
    }
    
    TIter stripItr(track->GetDaughterIterator());
    
    while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(stripItr())) {
      if(strip->GetNDigit(StripEnd::kWhole)==2 &&
         vetoplane[strip->GetPlane()]==0){
        
        TIter digitItr(strip->GetDaughterIterator());
        PlexStripEndId *plexstripendid=0;
        UgliGeomHandle *ugh = 0;
        PlexSEIdAltLItem *plexstripenditem=0;
        
        int cur_strip = strip->GetStrip();
        int cur_plane = strip->GetPlane();

        while (CandDigitHandle *digit = dynamic_cast<CandDigitHandle*>(digitItr())) {
          
          plexstripendid = 
            new PlexStripEndId(digit->GetPlexSEIdAltL().GetBestSEId());
          
          ugh = new UgliGeomHandle(*digit->GetVldContext());
          UgliStripHandle striphandle = ugh->GetStripHandle(*plexstripendid);

          plexstripenditem = 
            new PlexSEIdAltLItem(digit->GetPlexSEIdAltL().GetBestItem());

          DbiResultPtr<CalADCtoPE> resptr;
          resptr.NewQuery(vc,0);
          const CalADCtoPE* dpgc = resptr.GetRowByIndex(plexstripendid->BuildPlnStripEndKey());
          if(dpgc ==0) {
            MSG("MuCal",Msg::kWarning) 
              << "No database row for StripEnd " << *plexstripendid 
              << " (indexed as " << plexstripendid->BuildPlnStripEndKey() 
              << " )\n";
          }
          
          Calibrator::Instance().Reset(vc);

          float WLSLen = 0;  //green fibre length
          float ToZero = 0;  //green fibre length to middle of strip
          double plcor = 0;  //path length correction factor
          double attcor = 0;  //attenuation correction the the middle of strip

          switch (digit->GetPlexSEIdAltL().GetBestSEId().GetEnd()) {
          case StripEnd::kNegative: case StripEnd::kWhole: case StripEnd::kUnknown:

            if(fgain[cur_plane][cur_strip][0]==-1) {
              fgain[cur_plane][cur_strip][0]= dpgc->GetGain();
              fspewidth[cur_plane][cur_strip][0] = dpgc->GetSPEWidth();
            }

            // path length correction: 
            // <ds> = (ds/dz) * ( 4.050 / ((4.050/1.0) + |dT/dz|) ) 
            // dz/ds(plane) = tracksr->GetDirCosZ(plane)
            // dT/dz = du/dz or dv/dz depending on plane
            // du/ds(plane) =  tracksr->GetDirCosU(plane)
            // dv/ds(plane) =  tracksr->GetDirCosV(plane)
            
            WLSLen = striphandle.GetHalfLength() 
              + striphandle.WlsPigtail(StripEnd::kNegative);

            if (strip->GetPlaneView()==PlaneView::kU) {
              ToZero = -track->GetV(cur_plane);
              WLSLen += ToZero;
              attcor = TMath::Exp(ToZero/6.);
              plcor = 4.050 / (4.050 + fabs(tracksr->GetDirCosU(cur_plane)
                                            /tracksr->GetDirCosZ(cur_plane)));
            }

            else if (strip->GetPlaneView()==PlaneView::kV) {
              ToZero = track->GetU(cur_plane);
              WLSLen += ToZero;
              attcor = TMath::Exp(ToZero/6.);
              plcor = 4.050 / (4.050 + fabs(tracksr->GetDirCosV(cur_plane)
                                            /tracksr->GetDirCosZ(cur_plane)));
            }


            fhistoraw0->Fill(plexstripenditem->GetPE());

            fhistoang0->Fill(plexstripenditem->GetPE()
                             *tracksr->GetDirCosZ(cur_plane)/plcor);

            fhistoatt0->Fill(plexstripenditem->GetPE()
                             *attcor*tracksr->GetDirCosZ(cur_plane)/plcor);
            
            fHistoArray[cur_plane][cur_strip][0]
              ->Fill(plexstripenditem->GetPE()
                     *attcor*tracksr->GetDirCosZ(cur_plane)/plcor);
           
            break;
          

          case StripEnd::kPositive:

            if(fgain[cur_plane][cur_strip][1]==-1) {
              fgain[cur_plane][cur_strip][1]= dpgc->GetGain();
              fspewidth[cur_plane][cur_strip][1] = dpgc->GetSPEWidth();
            }

            WLSLen = striphandle.GetHalfLength() 
              + striphandle.WlsPigtail(StripEnd::kPositive);

            if (strip->GetPlaneView()==PlaneView::kU) {
              ToZero = track->GetV(cur_plane);
              WLSLen += ToZero;
              attcor = TMath::Exp(ToZero/6.);
              plcor = 4.050 / (4.050 + fabs(tracksr->GetDirCosU(cur_plane)
                                            /tracksr->GetDirCosZ(cur_plane)));
            }

            else if (strip->GetPlaneView()==PlaneView::kV) {
              ToZero = -track->GetU(cur_plane);
              WLSLen += ToZero;
              attcor = TMath::Exp(ToZero/6.);
              plcor = 4.050 / (4.050 + fabs(tracksr->GetDirCosV(cur_plane)
                                            /tracksr->GetDirCosZ(cur_plane)));
            }

            fhistoraw1->Fill(plexstripenditem->GetPE());

            fhistoang1->Fill(plexstripenditem->GetPE()
                             *tracksr->GetDirCosZ(cur_plane)/plcor);

            fhistoatt1->Fill(plexstripenditem->GetPE()
                             *attcor*tracksr->GetDirCosZ(cur_plane)/plcor);   

            fHistoArray[cur_plane][cur_strip][1]
              ->Fill(plexstripenditem->GetPE()
                     *attcor*tracksr->GetDirCosZ(cur_plane)/plcor);
            
            break;

          }
        }
      
      delete plexstripendid;
      delete ugh;

      }
    }

  }

  return result;

}

//......................................................................

void MuCalFitterModule::EndJob() 
{

  if(fwhichfit==1 || fwhichfit==2) 
    MSG("User",Msg::kInfo) << "Starting to Fit... " << endl;

  int numpar;
  if(fwhichfit==0) numpar = 0;
  else if(fwhichfit==1) numpar = 10;
  else if(fwhichfit==2) numpar = 12;
  else {
    numpar=0;
    MSG("User",Msg::kInfo) << "Not Fitting! No fit type specified" << endl;
  }

  ofstream mudat;
  mudat.open("muonfits.dat");

  TH1F *gains = new TH1F("gains","gains",250,-50,200);

  ffile->cd();

  for(int i=0;i<50;i++){
    for(int j=0;j<192;j++) {
      for(int k=0;k<2;k++) {
        
        gains->Fill(fgain[i][j][k]);

        if(fHistoArray[i][j][k]->GetEntries()>400){

          MuCalTrunc trunc(fHistoArray[i][j][k]);
          double *trunc_info = trunc.GetTruncMean();

          double *fit_info = 0;
          if(fwhichfit==1) fit_info = 
                             SpectraFit(fHistoArray[i][j][k],fgain[i][j][k]);
          else if(fwhichfit==2) 
            fit_info = BigSpectraFit(fHistoArray[i][j][k],
                                     fgain[i][j][k],fspewidth[i][j][k]);

          if (fWriteHistos) fHistoArray[i][j][k]->Write();
          
          mudat << i << " " << j << " " << k << " ";
          for(int i=0;i<numpar;i++) mudat << fit_info[i] << " ";
          mudat << trunc_info[0] << " " << trunc_info[1] << endl;
          
          delete fit_info;
          delete [] trunc_info;
        }

      }
    }
  }

  gains->Write();
  delete gains;

  for(int i=0;i<2;i++){

    double *fit_info = 0;

    MuCalTrunc trunc;

    if(i==0) trunc.SetHisto(fhistoang0);
    else trunc.SetHisto(fhistoang1);

    double *trunc_info = trunc.GetTruncMean();

    if(i==0) {
      if(fwhichfit==1) fit_info = SpectraFit(fhistoraw0,1);
      else if(fwhichfit==2) 
        fit_info = BigSpectraFit(fhistoraw0,1,0.5);
    }

    else {
      if(fwhichfit==1) fit_info = SpectraFit(fhistoraw1,1);
      else if(fwhichfit==2) 
        fit_info = BigSpectraFit(fhistoraw1,1,0.5);
    }

    mudat << "-1 " << "-1 " << i << " ";
    for(int i=0;i<numpar;i++) mudat << fit_info[i] << " ";
    mudat << trunc_info[0] << " " << trunc_info[1] << endl;

    delete trunc_info;
    delete fit_info;

  }
  
  if(fWriteHistos) {
    fhistoraw0->Write();
    fhistoraw1->Write();
    fhistoang0->Write();
    fhistoang1->Write();
    fhistoatt0->Write();
    fhistoatt1->Write();
  }

  mudat.close();  
  ffile->Close();
  
}

---