UberPlotsModule.cxx

Go to the documentation of this file.

// $Id: UberPlotsModule.cxx,v 1.5 2007/03/01 17:37:57 rhatcher Exp $
//
// Fills histograms with data from UberRecords
//
// kordosky@hep.utexas.edu
#include "CalDetDST/UberPlotsModule.h"

#include "TString.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TProfile.h"
#include "TProfile2D.h"

#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro

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

TFile* UberPlotsModule::fOut=0;

JOBMODULE(UberPlotsModule, "UberPlotsModule",
          "Make some plots from UberRecords");
CVSID("$Id: UberPlotsModule.cxx,v 1.5 2007/03/01 17:37:57 rhatcher Exp $");
//......................................................................

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

}

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

UberPlotsModule::~UberPlotsModule()
{
//======================================================================
// does nothing
//======================================================================

/*
  // seg faulted -- why?
     if(fOut) {
          delete fOut;
          fOut=0;
     }
*/
}

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

void UberPlotsModule::BeginJob()
{
//======================================================================
// opens a file to write histograms to
// creates some histograms
//======================================================================

MSG("CalDetDST", Msg::kDebug)<<"UberPlotsModule::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);

fOut->cd();
TDirectory* d = new TDirectory(fDirName.c_str(), fDirName.c_str());
d->cd();

// define histograms
h_mips=new TH1F("mips", "MIPs", 1000, 0, 2000);
h_mips24=new TH1F("mips24", "MIPs (first 24 planes)", 1000, 0, 2000);

h_ph = new TH1F("ph", "# Planes Hit", 60, -0.5, 59.5);

for(int i=0; i<60; i++){
    TString n1="strip_occ_neg_p";n1+=i;
    TString t1="Strip Occupancy (-) side Plane "; t1+=i;
    TString n2="strip_occ_pos_p";n2+=i;
    TString t2="Strip Occupancy (+) side Plane "; t2+=i;
    fStripOcc1[i] = new TH1F(n1.Data(),t1.Data(),24,-0.5,23.5);
    fStripOcc2[i] = new TH1F(n2.Data(),t2.Data(),24,-0.5,23.5);
}

h_pse1 = new TH2F("pse1", "Plane vs. Strip: Even Planes (-) side", 
                  60, -0.5, 59.5, 24, -0.5, 23.5);
h_pse2 = new TH2F("pse2", "Plane vs. Strip: Even Planes (+) side", 
                  60, -0.5, 59.5, 24, -0.5, 23.5);
h_pso1 = new TH2F("pso1", "Plane vs. Strip: Odd Planes (-) side", 
                  60, -0.5, 59.5, 24, -0.5, 23.5);
h_pso2 = new TH2F("pso2", "Plane vs. Strip: Odd Planes (+) side", 
                  60, -0.5, 59.5, 24, -0.5, 23.5);

h_wpse1 = new TH2F("wpse1", 
                   "Plane vs. Strip: Weighted by MIPs Even Planes (-) side", 
                   60, -0.5, 59.5, 24, -0.5, 23.5);
h_wpse2 = new TH2F("wpse2", 
                   "Plane vs. Strip: Weighted by MIPs Even Planes (+) side", 
                   60, -0.5, 59.5, 24, -0.5, 23.5);
h_wpso1 = new TH2F("wpso1", 
                   "Plane vs. Strip: Weighted by MIPs Odd Planes (-) side", 
                   60, -0.5, 59.5, 24, -0.5, 23.5);
h_wpso2 = new TH2F("wpso2", 
                   "Plane vs. Strip: Weighted by MIPs Odd Planes (+) side", 
                   60, -0.5, 59.5, 24, -0.5, 23.5);

h_ptrig = new TH1F("ptrig", "Plane Trigger", 7, -0.5, 6.5);
h_ptrig->GetXaxis()->SetTitle("n for n/n+1");

h_edep = new TProfile("edep", "MIPs vs. Plane", 60, -0.5, 59.5);

// Last thing:
// switch back to old add dir status so as not to hose user
TH1::AddDirectory(old_dir_stat);
}


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

void UberPlotsModule::Config(const Registry& r)
{
//======================================================================
// Configure the module given the Registry r
//======================================================================
  const char*    t;

  if (r.Get("DirName",t)) { fDirName = t; }

}

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

const Registry& UberPlotsModule::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);
  
  r.LockValues();

  return r;
}

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

void UberPlotsModule::EndJob()
{
//======================================================================
// write and close file
//======================================================================

     if(fOut->IsOpen()){
          fOut->Write();
          fOut->Close();
     }
}

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

JobCResult UberPlotsModule::Ana(const MomNavigator* mom)
{
//======================================================================
// does the plotting
//======================================================================
const UberRecord* ur = 
    dynamic_cast<const UberRecord*>(mom->GetFragment("UberRecord"));

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

// require it to be a pion, nooverlap

//if(!( (ur->cpid.pid&0x04)&&(ur->cpid.nov==1)) ) return JobCResult::kAOK;

// loop over hits the first time to fill array of mips in each plane
float mip_by_plane[60]={0.0};
const TClonesArray* hitlist = ur->GetHitList();
for(int i=0; i<hitlist->GetSize(); i++){
    const UberHit* hit= static_cast<const UberHit*>(hitlist->At(i));
    if(!hit) continue;
    
    const int plane=hit->GetPlane();
    if((plane<0)||(plane>59)) continue;
    const int strip=hit->GetStrip();
    if((strip<0)||(strip>23)) continue;

    mip_by_plane[plane]+=(hit->GetPosMIP()+hit->GetNegMIP());

}

map<double,int> pmap;
for(int i=0; i<60; i++){pmap.insert(pair<double,int>(mip_by_plane[i],i));}


// ok, now apply a plane trigger!

const int trig = PlaneTrigger(mip_by_plane, 60);
h_ptrig->Fill(trig);


// cut on the trigger

const int lowtrig=2; // n in n/n+1 ... should make configurable

// not failing it, just don't want to histogram it
if(trig<lowtrig) return JobCResult::kAOK; 

h_mips->Fill(ur->totmip);
h_ph->Fill(ur->nhitplanes);

for(int i=0; i<60; i++) h_edep->Fill(i, mip_by_plane[i]);

// compute # of mips in first 24 planes
double mips24=0.0;
for(int i=1; i<24; i++) mips24+=mip_by_plane[i];
h_mips24->Fill(mips24);

// loop over hits again
for(int i=0; i<hitlist->GetSize(); i++){
    const UberHit* hit= static_cast<const UberHit*>(hitlist->At(i));
    if(!hit) continue;
    
    const int plane=hit->GetPlane();
    if((plane<0)||(plane>59)) continue;
    const int strip=hit->GetStrip();
    if((strip<0)||(strip>23)) continue;
    
    // fill strip occupancy histograms
    if(hit->GetNegNPE()>0.0) {
        fStripOcc1[plane]->Fill(strip);
        if(plane%2==0) {
            h_pse1->Fill(plane,strip);
            h_wpse1->Fill(plane,strip,hit->GetNegMIP()/ur->totmip);
            }
        else {
            h_pso1->Fill(plane,strip);
            h_wpso1->Fill(plane,strip,hit->GetNegMIP()/ur->totmip);
            }
        }
    if(hit->GetPosNPE()>0.0) {
        fStripOcc2[plane]->Fill(strip);
        if(plane%2==0) {
            h_pse2->Fill(plane,strip);
            h_wpse2->Fill(plane,strip,hit->GetPosMIP()/ur->totmip);
            }
        else {
            h_pso2->Fill(plane,strip);
            h_wpso2->Fill(plane,strip,hit->GetPosMIP()/ur->totmip);
            }
        }
    }


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

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

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



int UberPlotsModule::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;
    }

---