CalDet2003PlotsModule.cxx

Go to the documentation of this file.

// $Id: CalDet2003PlotsModule.cxx,v 1.13 2008/06/12 19:23:53 rhatcher Exp $
//
// Fills histograms with data from UberRecords
//
// kordosky@hep.utexas.edu
#include <cstdio>
#include <cmath>
#include "CalDetDST/CalDet2003PlotsModule.h"
#include "TROOT.h"
#include "TString.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TProfile2D.h"
#include "TDirectory.h"
#include "TEventList.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro

#include "CalDetDST/UberRecord.h"
#include "CalDetDST/UberHit.h"

int CalDet2003PlotsModule::fcount = 0;
int CalDet2003PlotsModule::fNumObj = 0;
std::string CalDet2003PlotsModule::fFilePath = "";
std::string CalDet2003PlotsModule::fFileName = "CDPlots-default.root";
TFile* CalDet2003PlotsModule::fOut=0;

JOBMODULE(CalDet2003PlotsModule, "CalDet2003PlotsModule",
          "Make some plots from UberRecord");
CVSID("$Id: CalDet2003PlotsModule.cxx,v 1.13 2008/06/12 19:23:53 rhatcher Exp $");
//......................................................................

CalDet2003PlotsModule::CalDet2003PlotsModule()
{
//======================================================================
// does nothing
//======================================================================

   fNumObj++;

}

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

CalDet2003PlotsModule::~CalDet2003PlotsModule()
{
//======================================================================
// does nothing
//======================================================================
   MSG("CalDetDST",Msg::kDebug)<<"in CalDet2003PlotsModule() Destructor"<<endl
                               <<"dir name "<<fDirName<<endl;

   fNumObj--;
   if(fNumObj==0){
      if(fOut){
         if(fOut->IsOpen()){
//       fOut->cd();
            MSG("CalDetDST",Msg::kDebug)<<"cd'd to "<<fFileName<<endl;
//       fOut->ls();
            fOut->Write();
            MSG("CalDetDST",Msg::kDebug)<<"wrote file "<<fFileName<<endl;
            fOut->Close();
            MSG("CalDetDST",Msg::kDebug)<<"closed file "<<fFileName<<endl;
            delete fOut;
            MSG("CalDetDST",Msg::kDebug)<<"deleted file "<<fFileName<<endl;
            fOut=0; 
         }
      }
   }       
   //  MSG("CalDetDST",Msg::kDebug)<<"~CalDet2003PlotsModule: printing fDir: "<<fDir<<endl;
   //fDir->Print();
   delete fDir;
   fDir=0;


   delete el_snarls;
}

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

void CalDet2003PlotsModule::BeginJob()
{
//======================================================================
// opens a file to write histograms to
// creates some histograms
//======================================================================
   
   MSG("CalDetDST", Msg::kDebug)<<"CalDet2003PlotsModule::BeginJob()"<<endl;
/*   
   if(!fOut){
      fOut = new TFile("upm.root", "recreate");
   }
*/ 
// force histograms into the file
   Bool_t old_dir_stat = TH1::AddDirectoryStatus();
   TH1::AddDirectory(kTRUE);
   
   gROOT->cd();
   fDir = new TDirectory(fDirName.c_str(), fDirName.c_str());
   fDir->cd();
   
// define histograms
//tof histograms
   h_tdc0 = new TH1F("h_tdc0","TDC from TOF 0",4097,-0.5,4096.5);
   h_tdc1 = new TH1F("h_tdc1","TDC from TOF 1",4097,-0.5,4096.5);
   h_tdc2 = new TH1F("h_tdc2","TDC from TOF 2",4097,-0.5,4096.5);
   h_tdcdiff = new TH1F("h_tdcdiff","TDC2-TDC0",8193,-4096.5,4096.5);
   
//cerenkov histograms
   h_cera0 = new TH1F("h_cera0","CER0 ADC",500,0,5000);
   h_cera1 = new TH1F("h_cera1","CER1 ADC",500,0,5000);
   h_cera2 = new TH1F("h_cera2","CER2 ADC",500,0,5000);
   
   h_cert0 = new TH1F("h_cert0","CER0 Time",500,-250,250);
   h_cert1 = new TH1F("h_cert1","CER1 Time",500,-250,250);
   h_cert2 = new TH1F("h_cert2","CER2 Time",500,-250,250);
   
//hittime histograms
   h_hittime = new TH1F("h_hittime","Hit Times",801,-400.5,400.5);
   h_hittimenear = new TH1F("h_hittimenear","Hit Times, Near",801,-400.5,400.5);
   h_hittimefar = new TH1F("h_hittimefar","Hit Times, Far",801,-400.5,400.5);

   h_pid = new TH1F("h_pid","Particle Type",17,-0.5,16.5);
   h_nov = new TH1F("h_nov","Overlappers",2,0,2);
   h_inct = new TH1F("h_inct","Cerenkov in time",2,0,2);
   h_olchi2 = new TH1F("h_olchi2","Overlap Chi2",500,0,20);
   h_beamp = new TH1F("h_beamp","Beam Momentum",100,0,10);
   
//number of hits histograms
   h_nhits = new TH1F("h_nhits","Total Number of Hits",2881,-.5,2880.5);
   h_nhitsfarnear = new TH2F("h_nhitsfarnear","Total Number of Hits Near vs. Far",
                             2881,-.5,2880.5,2881,-.5,2880.5);
   h_nhitsfar = 0; //comes from projection of above
   h_nhitsnear = 0; //comes from projection of above
      
//number of hit strips histograms
   h_nhitstrips = new TH1F("h_nhitstrips","Total Number of Hit Strips",
                           201,-.5,200.5);
   h_nhitstripsfarnear = new TH2F("h_nhitstripsfarnear",
                                  "Total Number of Hit Strips Near vs. Far",
                                  201,-.5,200.5,201,-.5,200.5);
   h_nhitstripsfar = 0; //comes from projection of above
   h_nhitstripsnear = 0; //comes from projection of above

//number of hit planes histograms
   h_nhitplanes = new TH1F("h_nhitplanes","Total Number of Hit Planes",
                           61,-.5,60.5);
   h_nhitplanesfarnear = new TH2F("h_nhitplanesfarnear",
                                  "Total Number of Hit Planes Near vs. Far",
                                  61,-.5,60.5,61,-.5,60.5);
   h_nhitplanesfar = 0; //comes from projection of above
   h_nhitplanesnear = 0; //comes from projection of above
      
// strip vs. plane histograms
   h_svp = new TH2F("h_svp","Strip vs. Plane, both sides",60,0,60,24,0,24);
   h_svpfar = new TH2F("h_svpfar","Strip vs. Plane, Far",60,0,60,24,0,24);
   h_svpnear = new TH2F("h_svpnear","Strip vs. Plane, Near",60,0,60,24,0,24);

// weighted strip vs. plane histograms
   h_wsvp = new TH2F("h_wsvp","Weighted Strip vs. Plane, both sides",
                     60,0,60,24,0,24);
   h_wsvpfar = new TH2F("h_wsvpfar","Weighted Strip vs. Plane, Far",
                        60,0,60,24,0,24);
   h_wsvpnear = new TH2F("h_wsvpnear","Weighted Strip vs. Plane, Near",
                         60,0,60,24,0,24);

// total mip histograms
   h_tmip = new TH1F("h_tmip","Total MIP",140,0,700);
   h_tmipfarnear = new TH2F("h_tmipfarnear","Total MIP, Near vs. Far",
                            175,0,350,175,0,350);
   h_tmipfar = 0; //projection of above
   h_tmipnear = 0; //projection of above

//mip vs. plane distributions
   h_mipvplane = new TH2F("h_mipvplane","MIP vs. Plane, both sides",
                          60,0,60,50,0,100);
   h_mipvplanefar = new TH2F("h_mipvplanefar","MIP vs. Plane, Far",
                          60,0,60,50,0,100);
   h_mipvplanenear = new TH2F("h_mipvplanenear","MIP vs. Plane, Near",
                          60,0,60,50,0,100);

   h_mipvplane_pfx = 0; //profile of above
   h_mipvplanefar_pfx = 0; //profile of above
   h_mipvplanenear_pfx = 0; //profile of above
   
   
   h_mipvplanediff = new TH2F("h_mipvplanediff","MIP (Far-Near) vs. Plane",
                              60,0,60,500,-50,50);
   h_mipvplaneratio = new TH2F("h_mipvplaneratio","Ratio MIP (Near/Far) vs. Plane",
                               60,0,60,100,-1,1);

   h_mipvplanediff_pfx = 0; //profile of above
   h_mipvplaneratio_pfx = 0; //profile of above

//hit adc histograms
   h_hitadc = new TH1F("h_hitadc","Hit ADC",170,0,17000);
   h_hitadcfarnear = new TH2F("h_hitadcfarnear","Hit ADC Near vs. Far",
                              170,0,17000,170,0,17000);
   h_hitadcfar = 0; //projection of above
   h_hitadcnear = 0; //projection of above
      
   h_fartime = new TH1F("h_fartime","Far hit time-toftime", 1000,-500,500);
   h_neartime = new TH1F("h_neartime","Near hit time-toftime", 1000,-500,500);

   h_triggertime = new TH1F("h_triggertime","Trigger time",
                               9000,0,1800);
   h_triggertimenofar = new TH1F("h_triggertimenofar","Trigger time No Far",
                               9000,0,1800);
   h_triggertimenonear = new TH1F("h_triggertimenonear","Trigger time No Near",
                               9000,0,1800);

   // a couple of event lists

   el_snarls=new TEventList("snarls", "Offline snarl number");
   el_snarls->SetDirectory(0);

//   MSG("CalDetDST",Msg::kDebug)<<"gROOT::ls BeginJob"<<endl;
//   gROOT->ls();
//   MSG("CalDetDST",Msg::kDebug)<<"fDir::ls BeginJob"<<endl;
//   fDir->ls();
// Last thing:
// switch back to old add dir status so as not to hose user
TH1::AddDirectory(old_dir_stat);
}


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

void CalDet2003PlotsModule::Config(const Registry& r)
{
//======================================================================
// Configure the module given the Registry r
//======================================================================
  const char*    t;
  if (r.Get("DirName",t)) { fDirName = t; }

  const char *f;
  if(r.Get("OutFilePath",f)) { fFilePath = f; }

  if(r.Get("FileName",f)){ fFileName = f;}
}

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

const Registry& CalDet2003PlotsModule::DefaultConfig() const
{
//======================================================================
// Supply the default configuration for the module
//======================================================================
  static Registry r; // Default configuration for module

  // Set name of config
  std::string name = this->JobCModule::GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  // Set values in configuration
  r.UnLockValues();

  const char* s = "all";
  r.Set("DirName", s);
  
  const char* f = "./";
  r.Set("OutFilePath", f);

  const char *n = "CDPlots-default.root";
  r.Set("FileName",n);
  r.LockValues();

  return r;
}

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

void CalDet2003PlotsModule::EndJob()
{
//======================================================================
// write and close file
//======================================================================
   MSG("CalDetDST",Msg::kDebug)<<"in CalDet2003PlotsModule() EndJob()"<<endl
                               <<"dir name "<<fDirName<<endl;

   if(fFileName==""){
      MSG("CalDetDST",Msg::kError)<<"FileName not set.  Can Not write file."<<endl;
      return;
   }
      


   if(!fOut){
      MSG("CalDetDST",Msg::kDebug)<<"Making output file "<<fFileName<<endl;
      string fullfn = fFilePath+fFileName;
      fOut = new TFile(fullfn.c_str(), "recreate");
   }


   fDir->cd();

   MSG("CalDetDST",Msg::kDebug)<<"EndJob() cd'd"<<endl;


//make projection and profile histograms
   h_nhitsfar = h_nhitsfarnear->ProjectionX(); 
   h_nhitsfar->SetNameTitle("h_nhitsfar","Number of Hits, Far");
   h_nhitsnear = h_nhitsfarnear->ProjectionY(); 
   h_nhitsnear->SetNameTitle("h_nhitsnear","Number of Hits, Near");

   MSG("CalDetDST",Msg::kDebug)<<"EndJob() made nhits proj"<<endl;

   h_nhitstripsfar = h_nhitstripsfarnear->ProjectionX();
   h_nhitstripsfar->SetNameTitle("h_nhitstripsfar","Number of Hit Strips, Far");

   h_nhitstripsnear = h_nhitstripsfarnear->ProjectionY(); 
   h_nhitstripsnear->SetNameTitle("h_nhitstripsnear","Number of Hit Strips, Near");

   MSG("CalDetDST",Msg::kDebug)<<"EndJob() made nhitstrips proj"<<endl;

   h_nhitplanesfar = h_nhitplanesfarnear->ProjectionX();
   h_nhitplanesfar->SetNameTitle("h_nhitplanesfar","Number of Hit Planes, Far");

   h_nhitplanesnear = h_nhitplanesfarnear->ProjectionY(); 
   h_nhitplanesnear->SetNameTitle("h_nhitplanesnear","Number of Hit Planes, Near");

   MSG("CalDetDST",Msg::kDebug)<<"EndJob() made nhitplanes proj"<<endl;

   h_tmipfar = h_tmipfarnear->ProjectionX();
   h_tmipfar->SetNameTitle("h_tmipfar","Total MIP, Far");

   h_tmipnear = h_tmipfarnear->ProjectionY();
   h_tmipnear->SetNameTitle("h_tmipnear","Total MIP, Near");

   MSG("CalDetDST",Msg::kDebug)<<"EndJob() made tmip proj"<<endl;

   h_mipvplane_pfx = h_mipvplane->ProfileX();
   h_mipvplane_pfx->SetNameTitle("h_mipvplane_pfx","MIP vs. Plane, Profile");

   h_mipvplanefar_pfx = h_mipvplanefar->ProfileX();
   h_mipvplanefar_pfx->SetNameTitle("h_mipvplanefar_pfx",
                                    "MIP vs. Plane Profile, Far");

   h_mipvplanenear_pfx = h_mipvplanenear->ProfileX();
   h_mipvplanenear_pfx->SetNameTitle("h_mipvplanenear_pfx",
                                    "MIP vs. Plane Profile, Near");
   
   h_mipvplanediff_pfx = h_mipvplanediff->ProfileX();
   h_mipvplanediff_pfx->SetNameTitle("h_mipvplanediff_pfx",
                                     "MIP vs. Plane diff Far-Near, Profile");

   h_mipvplaneratio_pfx = h_mipvplaneratio->ProfileX();
   h_mipvplaneratio_pfx->SetNameTitle("h_mipvplaneratio_pfx",
                                      "MIP vs. Plane Near/Far, Profile");

   MSG("CalDetDST",Msg::kDebug)<<"EndJob() made profs"<<endl;

   h_hitadcfar = h_hitadcfarnear->ProjectionX();
   h_hitadcfar->SetNameTitle("h_hitadcfar","Hit ADC, Far");

   h_hitadcnear = h_hitadcfarnear->ProjectionY();
   h_hitadcnear->SetNameTitle("h_hitadcnear","Hit ADC, Near");

   MSG("CalDetDST",Msg::kDebug)<<"EndJob() made adc proj"<<endl;

   fOut->cd();
   TDirectory *newdir=new TDirectory(fDirName.c_str(),fDirName.c_str());
   newdir->cd();
//   MSG("CalDetDST",Msg::kDebug)<<"gROOT::ls"<<endl;
//   gROOT->ls();
//   MSG("CalDetDST",Msg::kDebug)<<"fDir::ls"<<endl;
//   fDir->ls();
   fDir->GetList()->Write();

   // jam it for the stupid event list
   el_snarls->Write();


//   newdir->ls();
//   fDir->ls();
//   fOut->Append(fDir);
   
//   fDir->Write();
/*
   MSG("CalDetDST",Msg::kDebug)<<"Appeneded directory "
                               <<fDirName<<" to "<<fFileName<<endl;
   MSG("CalDetDST",Msg::kDebug)<<"Directory "<<fDirName
                               <<" IsWritable? "<<fDir->IsWritable()<<endl;
   fDir->SetWritable();
   MSG("CalDetDST",Msg::kDebug)<<"Directory "<<fDirName
                               <<" IsWritable? "<<fDir->IsWritable()<<endl;
*/
}

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

JobCResult CalDet2003PlotsModule::Ana(const MomNavigator* mom)
{
//======================================================================
// does the plotting
//======================================================================

if(fcount%1000==0){
   MSG("CalDetDST",Msg::kInfo)<<"On Event "<<fcount<<endl;
   // MAK: Feb 15, 2005
   // the following should be commented out in the CVS version
   //   system("date"); // list date
   //   system("tokens"); // list afs tokens
}
fcount++;

const UberRecord* ur = 
    dynamic_cast<const UberRecord*>(mom->GetFragment("UberRecord"));

if(!ur)   return JobCResult::kAOK; // don't care if there isn't a record!

/*
if(fFileName=="CDPlots-default"){
   MSG("CalDetDST",Msg::kDebug)<<"Trying to set filename"<<endl;
   const Int_t runnum = (ur->GetHeader()).GetRunNo();
   const Int_t subrunnum = (ur->GetHeader()).GetSubRunNo();
   char fn[100];
   sprintf(fn,"CDPlots-%d-%d.root",runnum,subrunnum);
   fFileName = fn;
   MSG("CalDetDST",Msg::kDebug)<<"Filename will be "<<fFileName<<endl;
}
*/

// fill event list
el_snarls->Enter(ur->snarlno);

//fill tdc hists
h_tdc0->Fill(ur->toftdc[0]);
h_tdc1->Fill(ur->toftdc[1]);
h_tdc2->Fill(ur->toftdc[2]);
h_tdcdiff->Fill(ur->toftdc[2]-ur->toftdc[0]);

//fill cer hists
h_cera0->Fill(ur->ceradc[0]);
h_cera1->Fill(ur->ceradc[1]);
h_cera2->Fill(ur->ceradc[2]);
if(ur->ceradc[0]>0){
   h_cert0->Fill(ur->certime[0]-ur->toftime);
}
if(ur->ceradc[1]>0){
   h_cert1->Fill(ur->certime[1]-ur->toftime);
}
if(ur->ceradc[2]>0){
   h_cert2->Fill(ur->certime[2]-ur->toftime);
}

h_pid->Fill(ur->cpid.pid);
h_nov->Fill(ur->cpid.nov);
h_inct->Fill(ur->cpid.inct);
h_olchi2->Fill(ur->cpid.olchi2);
h_beamp->Fill(ur->GetHeader().GetBeamMomentum());

const TClonesArray* hitlist = ur->GetHitList();
int prevplane = -1;
float planemip = 0.;
float planemipfar = 0.;
float planemipnear = 0.;
float tmipfar = 0.;
float tmipnear = 0.;
int nhitsfar = 0;
int nhitsnear = 0;
int nhitstripsnear = 0;
int nhitstripsfar = 0;
int nhitplanesnear = 0;
int nhitplanesfar = 0;
Bool_t planehitnear = kFALSE;
Bool_t planehitfar = kFALSE;
for(int i=0; i<ur->nhitstrips; i++){
    const UberHit* uh = static_cast<const UberHit*>(hitlist->At(i));
    if(!uh) continue;
    
    //should never fail these by construction of hitlist, but . . .
    const int plane=uh->GetPlane();
    if((plane<0)||(plane>59)) continue;
    const int strip=uh->GetStrip();
    if((strip<0)||(strip>23)) continue;

    if(prevplane!=uh->GetPlane()&&prevplane!=-1){
       h_mipvplane->Fill(prevplane,planemip);
       h_mipvplanefar->Fill(prevplane,planemipfar);
       h_mipvplanenear->Fill(prevplane,planemipnear);
       h_mipvplanediff->Fill(prevplane,(planemipfar-planemipnear));
       if(planemipfar+planemipnear!=0){
          h_mipvplaneratio->Fill(prevplane,
                                 (planemipfar-planemipnear)/((planemipnear+planemipfar)/2.));
       }     
       if(planehitnear){
          nhitplanesnear++;
       }
       if(planehitfar){
          nhitplanesfar++;
       }
       planemip=0.;
       planemipfar=0.;
       planemipnear=0.;
       planehitnear = kFALSE;
       planehitfar = kFALSE;
    }
    prevplane = uh->GetPlane();
    planemip+=(uh->GetPosMIP()+uh->GetNegMIP());
    planemipfar+=uh->GetNegMIP();
    planemipnear+=uh->GetPosMIP();
    if(uh->GetPlane()!=0){
       tmipfar+=uh->GetNegMIP();
       tmipnear+=uh->GetPosMIP();
    }
    h_hitadc->Fill(uh->GetNegADC());
    h_hitadc->Fill(uh->GetPosADC());
    h_hitadcfarnear->Fill(uh->GetNegADC(),uh->GetPosADC());
    h_svp->Fill(uh->GetPlane(),uh->GetStrip());
    h_wsvp->Fill(uh->GetPlane(),uh->GetStrip(),uh->GetPosMIP()+uh->GetNegMIP());

    if(uh->GetPosADC()>0){
       h_hittime->Fill(uh->GetPosTime()-ur->toftime);
       h_hittimenear->Fill(uh->GetPosTime()-ur->toftime);
       h_svpnear->Fill(uh->GetPlane(),uh->GetStrip());
       h_wsvpnear->Fill(uh->GetPlane(),uh->GetStrip(),uh->GetPosMIP());
       if(uh->GetPlane()!=0){
          h_neartime->Fill(uh->GetPosTime()-ur->toftime);
       }
       nhitsnear++;
       nhitstripsnear++;
       planehitnear=kTRUE;
    }
    if(uh->GetNegADC()>0){
       h_hittime->Fill(uh->GetNegTime()-ur->toftime);
       h_hittimefar->Fill(uh->GetNegTime()-ur->toftime);
       h_svpfar->Fill(uh->GetPlane(),uh->GetStrip());
       h_wsvpfar->Fill(uh->GetPlane(),uh->GetStrip(),uh->GetNegMIP());
       if(uh->GetPlane()!=0){
          h_fartime->Fill(uh->GetNegTime()-ur->toftime);
       }
       nhitsfar++;
       nhitstripsfar++;
       planehitfar=kTRUE;
    }
}
//fill nhits hists
h_nhits->Fill(ur->nhits);
h_nhitsfarnear->Fill(nhitsfar,nhitsnear);

//fill nhitstrips hists
h_nhitstrips->Fill(ur->nhitstrips);
h_nhitstripsfarnear->Fill(nhitstripsfar,nhitstripsnear);

//fill nhitplanes hists
h_nhitplanes->Fill(ur->nhitplanes);
h_nhitplanesfarnear->Fill(nhitplanesfar,nhitplanesnear);

//fill total mip hists
h_tmip->Fill(ur->totmip);
h_tmipfarnear->Fill(tmipfar,tmipnear);

h_triggertime->Fill(ur->triggertime);
if(tmipfar==0){
   h_triggertimenofar->Fill(ur->triggertime);
}
if(tmipnear==0){
   h_triggertimenonear->Fill(ur->triggertime);
}

  return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

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

void CalDet2003PlotsModule::Help()
{
//======================================================================
// I should print something here....
//======================================================================
}



int CalDet2003PlotsModule::PlaneTrigger(float* qvec, int len)
    {
    // qvec= array of charge,pes,mips by plane
    // len= the length of qvec
    const float thresh = 0.0; // threshold to count a plane
    int trig=0;
    int lown=0; // n-1 for the lowest n to try (eg lown=1 implies 2/3)
    for(int i=0; i<len; i++){
        // form a window with j+1 elements
        // j=n in n/n+1
        for(int j=6; j>lown; j--){
            int count=0;
            // loop over the window, careful not to go overboard
            for(int k=i; (k<=i+j && k<len); k++){
                if(qvec[k]>thresh) count++;
                }
            // don't assign if we already have a more restrictive trigger
            if(count>=j && count>trig) { 
                trig=j; // assign the trigger
                // bump the window so we only look for more restrictive triggers 
                lown=j;
                }
            }
        }
    return trig;
    }

---