CMS 3D CMS Logo

Public Member Functions | Private Member Functions | Private Attributes

ZdcTestAnalysis Class Reference

#include <ZdcTestAnalysis.h>

Inheritance diagram for ZdcTestAnalysis:
SimWatcher Observer< const BeginOfJob * > Observer< const BeginOfRun * > Observer< const EndOfRun * > Observer< const BeginOfEvent * > Observer< const EndOfEvent * > Observer< const G4Step * >

List of all members.

Public Member Functions

 ZdcTestAnalysis (const edm::ParameterSet &p)
virtual ~ZdcTestAnalysis ()

Private Member Functions

void finish ()
void update (const EndOfRun *run)
 This routine will be called when the appropriate signal arrives.
void update (const G4Step *step)
 This routine will be called when the appropriate signal arrives.
void update (const EndOfEvent *evt)
 This routine will be called when the appropriate signal arrives.
void update (const BeginOfEvent *evt)
 This routine will be called when the appropriate signal arrives.
void update (const BeginOfJob *run)
 This routine will be called when the appropriate signal arrives.
void update (const BeginOfRun *run)
 This routine will be called when the appropriate signal arrives.

Private Attributes

int doNTzdcevent
int doNTzdcstep
int eventIndex
std::string eventNtFileName
int stepIndex
std::string stepNtFileName
int verbosity
Float_t zdceventarray [16]
TNtuple * zdceventntuple
TFile * zdcOutputEventFile
TFile * zdcOutputStepFile
Float_t zdcsteparray [18]
TNtuple * zdcstepntuple

Detailed Description

Definition at line 65 of file ZdcTestAnalysis.h.


Constructor & Destructor Documentation

ZdcTestAnalysis::ZdcTestAnalysis ( const edm::ParameterSet p)

Definition at line 25 of file ZdcTestAnalysis.cc.

References gather_cfg::cout, doNTzdcevent, doNTzdcstep, eventNtFileName, edm::ParameterSet::getParameter(), stepNtFileName, verbosity, zdceventntuple, and zdcstepntuple.

                                                        {
  //constructor
  edm::ParameterSet m_Anal = p.getParameter<edm::ParameterSet>("ZdcTestAnalysis");
  verbosity    = m_Anal.getParameter<int>("Verbosity");
  doNTzdcstep  = m_Anal.getParameter<int>("StepNtupleFlag");
  doNTzdcevent  = m_Anal.getParameter<int>("EventNtupleFlag");
  stepNtFileName = m_Anal.getParameter<std::string>("StepNtupleFileName");
  eventNtFileName = m_Anal.getParameter<std::string>("EventNtupleFileName");
   
   if (verbosity > 0)
   std::cout<<std::endl;
   std::cout<<"============================================================================"<<std::endl;
   std::cout << "ZdcTestAnalysis:: Initialized as observer"<< std::endl;
   if (doNTzdcstep  > 0){
     std::cout <<" Step Ntuple will be created"<< std::endl;
     std::cout <<" Step Ntuple file: "<<stepNtFileName<<std::endl;
   }
   if (doNTzdcevent > 0){
     std::cout <<" Event Ntuple will be created"<< std::endl;
     std::cout <<" Step Ntuple file: "<<stepNtFileName<<std::endl;
   }
   std::cout<<"============================================================================"<<std::endl;
   std::cout<<std::endl;

   if (doNTzdcstep  > 0)
     zdcstepntuple = 
       new TNtuple("NTzdcstep","NTzdcstep",
                   "evt:trackid:charge:pdgcode:x:y:z:stepl:stepe:eta:phi:vpx:vpy:vpz:idx:idl:pvtype:ncherphot");

   if (doNTzdcevent  > 0)
     zdceventntuple = 
       new TNtuple("NTzdcevent","NTzdcevent",
                   "evt:ihit:fiberid:zside:subdet:layer:fiber:channel:enem:enhad:hitenergy:x:y:z:time:etot");

   //theZdcSD = new ZdcSD("ZDCHITSB", new ZdcNumberingScheme());
}
ZdcTestAnalysis::~ZdcTestAnalysis ( ) [virtual]

Definition at line 62 of file ZdcTestAnalysis.cc.

References gather_cfg::cout, finish(), and verbosity.

                                  {
  // destructor
  finish();
  if (verbosity > 0) {
    std::cout << std::endl << "ZdcTestAnalysis Dextructor  -------->  End of ZdcTestAnalysis : "
      << std::endl << std::endl; 
  }

  //if (doNTzdcstep  > 0)delete zdcstepntuple;
  //if (doNTzdcevent > 0)delete zdceventntuple;

  std::cout<<"ZdcTestAnalysis: End of process"<<std::endl;
}

Member Function Documentation

void ZdcTestAnalysis::finish ( ) [private]

Definition at line 389 of file ZdcTestAnalysis.cc.

References gather_cfg::cout, doNTzdcevent, doNTzdcstep, eventIndex, zdceventntuple, zdcOutputEventFile, zdcOutputStepFile, and zdcstepntuple.

Referenced by ~ZdcTestAnalysis().

                            {
  if (doNTzdcstep) {
    zdcOutputStepFile->cd();
    zdcstepntuple->Write();
    std::cout << "ZdcTestAnalysis: Ntuple step  written for event: "<<eventIndex<<std::endl;
    zdcOutputStepFile->Close();
    std::cout << "ZdcTestAnalysis: Step file closed" << std::endl;
  }
  

 if (doNTzdcevent) {
   zdcOutputEventFile->cd();
   zdceventntuple->Write("",TObject::kOverwrite);
   std::cout << "ZdcTestAnalysis: Ntuple event written for event: "<<eventIndex<<std::endl;   
   zdcOutputEventFile->Close();
   std::cout << "ZdcTestAnalysis: Event file closed" << std::endl;
 }
 
}
void ZdcTestAnalysis::update ( const EndOfRun ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfRun * >.

Definition at line 387 of file ZdcTestAnalysis.cc.

{;}
void ZdcTestAnalysis::update ( const G4Step *  ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const G4Step * >.

Definition at line 115 of file ZdcTestAnalysis.cc.

References gather_cfg::cout, doNTzdcstep, eta(), eventIndex, UserOptions_cff::idx, findQualityFiles::jj, funct::log(), ntzdcs_charge, ntzdcs_eta, ntzdcs_evt, ntzdcs_idl, ntzdcs_idx, ntzdcs_ncherphot, ntzdcs_pdgcode, ntzdcs_phi, ntzdcs_pvtype, ntzdcs_stepe, ntzdcs_stepl, ntzdcs_trackid, ntzdcs_vpx, ntzdcs_vpy, ntzdcs_vpz, ntzdcs_x, ntzdcs_y, ntzdcs_z, phi, stepIndex, verbosity, zdcsteparray, and zdcstepntuple.

                                                 {
  //step;
  stepIndex++;

  if (doNTzdcstep) {
    G4StepPoint * preStepPoint = aStep->GetPreStepPoint();
    // G4StepPoint * postStepPoint= aStep->GetPostStepPoint();
    G4double stepL = aStep->GetStepLength();
    G4double stepE = aStep->GetTotalEnergyDeposit();
    
    if (verbosity >= 2)
      std::cout << "Step " << stepL << "," << stepE <<std::endl;
    
    G4Track * theTrack    = aStep->GetTrack();
    G4int theTrackID      = theTrack->GetTrackID();
    G4double theCharge    = theTrack->GetDynamicParticle()->GetCharge();
    G4String particleType = theTrack->GetDefinition()->GetParticleName();
    G4int pdgcode         = theTrack->GetDefinition()->GetPDGEncoding();
    
    G4ThreeVector vert_mom = theTrack->GetVertexMomentumDirection();
    G4double vpx = vert_mom.x();
    G4double vpy = vert_mom.y();
    G4double vpz = vert_mom.z();
    double eta = 0.5 * log( (1.+vpz) / (1.-vpz) );
    double phi = atan2(vpy,vpx);
    
    G4ThreeVector hitPoint = preStepPoint->GetPosition();
    G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()->
    GetTopTransform().TransformPoint(hitPoint);
    
    const G4VTouchable* touch = aStep->GetPreStepPoint()->GetTouchable();
    int idx = touch->GetReplicaNumber(0);
    int idLayer = -1;
    int thePVtype = -1;
    
    int historyDepth = touch->GetHistoryDepth();

    if (historyDepth > 0) {
      std::vector<int>                theReplicaNumbers(historyDepth);
      std::vector<G4VPhysicalVolume*> thePhysicalVolumes(historyDepth);
      std::vector<G4String>           thePVnames(historyDepth);
      std::vector<G4LogicalVolume*>   theLogicalVolumes(historyDepth);
      std::vector<G4String>           theLVnames(historyDepth);
      std::vector<G4Material*>        theMaterials(historyDepth);
      std::vector<G4String>           theMaterialNames(historyDepth);
      
      for (int jj = 0; jj < historyDepth; jj++) {
        theReplicaNumbers[jj] = touch->GetReplicaNumber(jj);
        thePhysicalVolumes[jj] = touch->GetVolume(jj);
        thePVnames[jj] = thePhysicalVolumes[jj]->GetName();
        theLogicalVolumes[jj] = thePhysicalVolumes[jj]->GetLogicalVolume();
        theLVnames[jj] = theLogicalVolumes[jj]->GetName();
        theMaterials[jj] = theLogicalVolumes[jj]->GetMaterial();
        theMaterialNames[jj] = theMaterials[jj]->GetName();
        if (verbosity >= 2)
          std::cout << " GHD " << jj << ": " << theReplicaNumbers[jj] << ","
                       << thePVnames[jj] << "," << theLVnames[jj] << ","
                    << theMaterialNames[jj]  << std::endl;
      }

      idLayer = theReplicaNumbers[1];
      if (thePVnames[0] == "ZDC_EMLayer")
        thePVtype = 1;
      else if (thePVnames[0] == "ZDC_EMAbsorber")
        thePVtype = 2;
      else if (thePVnames[0] == "ZDC_EMFiber")
        thePVtype = 3;
      else if (thePVnames[0] == "ZDC_HadLayer")
        thePVtype = 7;
      else if (thePVnames[0] == "ZDC_HadAbsorber")
        thePVtype = 8;
      else if (thePVnames[0] == "ZDC_HadFiber")
        thePVtype = 9;
      else if (thePVnames[0] == "ZDC_PhobosLayer")
        thePVtype = 11;
      else if (thePVnames[0] == "ZDC_PhobosAbsorber")
        thePVtype = 12;
      else if (thePVnames[0] == "ZDC_PhobosFiber")
      thePVtype = 13;
      else {
        thePVtype = 0;
        if (verbosity >= 2)
          std::cout << " pvtype=0 hd=" << historyDepth << " " << theReplicaNumbers[0] << ","
                    << thePVnames[0] << "," << theLVnames[0] << "," << theMaterialNames[0] << std::endl;
      }
    }    
    else if (historyDepth == 0) { 
      int theReplicaNumber = touch->GetReplicaNumber(0);
      G4VPhysicalVolume* thePhysicalVolume = touch->GetVolume(0);
      G4String thePVname = thePhysicalVolume->GetName();
      G4LogicalVolume * theLogicalVolume = thePhysicalVolume->GetLogicalVolume();
      G4String theLVname = theLogicalVolume->GetName();
      G4Material * theMaterial = theLogicalVolume->GetMaterial();
      G4String theMaterialName = theMaterial->GetName();
      if (verbosity >= 2)
        std::cout << " hd=0 " << theReplicaNumber << "," 
                  << thePVname << "," << theLVname << "," 
                  << theMaterialName << std::endl;
    }
    else {
      std::cout << " hd<0:  hd=" << historyDepth << std::endl;
    }

    
    double NCherPhot = -1;
    zdcsteparray[ntzdcs_evt] = (float)eventIndex;
    zdcsteparray[ntzdcs_trackid] = (float)theTrackID;
    zdcsteparray[ntzdcs_charge] = theCharge;
    zdcsteparray[ntzdcs_pdgcode] = (float)pdgcode;
    zdcsteparray[ntzdcs_x] = hitPoint.x();
    zdcsteparray[ntzdcs_y] = hitPoint.y();
    zdcsteparray[ntzdcs_z] = hitPoint.z();
    zdcsteparray[ntzdcs_stepl] = stepL;
    zdcsteparray[ntzdcs_stepe] = stepE;
    zdcsteparray[ntzdcs_eta] = eta;
    zdcsteparray[ntzdcs_phi] = phi;
    zdcsteparray[ntzdcs_vpx] = vpx;
    zdcsteparray[ntzdcs_vpy] = vpy;
    zdcsteparray[ntzdcs_vpz] = vpz;
    zdcsteparray[ntzdcs_idx] = (float)idx;
    zdcsteparray[ntzdcs_idl] = (float)idLayer;
    zdcsteparray[ntzdcs_pvtype] = thePVtype;
    zdcsteparray[ntzdcs_ncherphot] = NCherPhot;
    zdcstepntuple->Fill(zdcsteparray);

  }
}
void ZdcTestAnalysis::update ( const EndOfEvent ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const EndOfEvent * >.

Definition at line 243 of file ZdcTestAnalysis.cc.

References gather_cfg::cout, doNTzdcevent, eventIndex, CaloG4Hit::getEM(), CaloG4Hit::getEnergyDeposit(), CaloG4Hit::getHadr(), CaloG4Hit::getPosition(), CaloG4Hit::getTimeSliceID(), CaloG4Hit::getTrackID(), CaloG4Hit::getUnitID(), i, Association::map, npart, ntzdce_channel, ntzdce_enem, ntzdce_enhad, ntzdce_etot, ntzdce_evt, ntzdce_fiber, ntzdce_fiberid, ntzdce_hitenergy, ntzdce_ihit, ntzdce_layer, ntzdce_subdet, ntzdce_time, ntzdce_x, ntzdce_y, ntzdce_z, ntzdce_zside, funct::pow(), mathSSE::sqrt(), stepIndex, cond::rpcobgas::time, ZdcNumberingScheme::unpackZdcIndex(), verbosity, zdceventarray, and zdceventntuple.

                                                   {
  //end of event 
   
  // Look for the Hit Collection
  std::cout  << std::endl << "ZdcTest::upDate(const EndOfEvent * evt) - event #" << (*evt)()->GetEventID()
             << std::endl << "  # of aSteps followed in event = " << stepIndex << std::endl;
  
  // access to the G4 hit collections
  G4HCofThisEvent* allHC = (*evt)()->GetHCofThisEvent();
  std::cout << "  accessed all HC";
  
  int theZDCHCid = G4SDManager::GetSDMpointer()->GetCollectionID("ZDCHITS");
  std::cout << " - theZDCHCid = " << theZDCHCid;
  
  CaloG4HitCollection* theZDCHC = (CaloG4HitCollection*) allHC->GetHC(theZDCHCid);
  std::cout << " - theZDCHC = " << theZDCHC << std::endl;
  
  ZdcNumberingScheme * theZdcNumScheme = new ZdcNumberingScheme(1);
  
  float ETot=0., SEnergy=0.;
  int maxTime=0;
  int fiberID=0;
  unsigned int unsignedfiberID=0;
  std::map<int,float,std::less<int> > energyInFibers;
  std::map<int,float,std::less<int> > primaries;
  float totalEnergy = 0;
  int nentries = theZDCHC->entries();
  std::cout << "  theZDCHC has " << nentries << " entries" << std::endl;

  if (doNTzdcevent) {
    if (nentries > 0) {
      for (int ihit = 0; ihit < nentries; ihit++) {
        CaloG4Hit* caloHit = (*theZDCHC)[ihit];
        totalEnergy += caloHit->getEnergyDeposit();
      }

      for (int ihit = 0; ihit < nentries; ihit++) {
        CaloG4Hit* aHit = (*theZDCHC)[ihit];
            fiberID = aHit->getUnitID();
            unsignedfiberID = aHit->getUnitID();
            double enEm = aHit->getEM();
            double enHad = aHit->getHadr();
            math::XYZPoint hitPoint = aHit->getPosition();
            double hitEnergy = aHit->getEnergyDeposit();
            if (verbosity >= 1)
              std::cout << " entry #" << ihit << ": fiberID=0x" << std::hex 
                        << fiberID << std::dec << "; enEm=" << enEm 
                        << "; enHad=" << enHad << "; hitEnergy=" << hitEnergy
                        << "z=" << hitPoint.z() << std::endl;
            energyInFibers[fiberID]+= enEm + enHad;
            primaries[aHit->getTrackID()]+= enEm + enHad;
            float time = aHit->getTimeSliceID();
            if (time > maxTime) maxTime=(int)time;
            
            int thesubdet, thelayer, thefiber, thechannel, thez;
            theZdcNumScheme->unpackZdcIndex(fiberID, thesubdet, thelayer, thefiber, thechannel, thez);
            int unsignedsubdet, unsignedlayer, unsignedfiber, unsignedchannel, unsignedz;
            theZdcNumScheme->unpackZdcIndex(unsignedfiberID, unsignedsubdet, 
                                            unsignedlayer, unsignedfiber, unsignedchannel, unsignedz);
            
            // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
            // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
            // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);
            // unsigned int packidx1 = packZdcIndex(thesubdet, thelayer, thefiber, thechannel, thez);

            zdceventarray[ntzdce_evt] = (float)eventIndex;
            zdceventarray[ntzdce_ihit] = (float)ihit;
            zdceventarray[ntzdce_fiberid] = (float)fiberID;
            zdceventarray[ntzdce_zside] = (float)thez;
            zdceventarray[ntzdce_subdet] = (float)thesubdet;
            zdceventarray[ntzdce_layer] = (float)thelayer;
            zdceventarray[ntzdce_fiber] = (float)thefiber;
            zdceventarray[ntzdce_channel] = (float)thechannel;
            zdceventarray[ntzdce_enem] = enEm;
            zdceventarray[ntzdce_enhad] = enHad;
            zdceventarray[ntzdce_hitenergy] = hitEnergy;
            zdceventarray[ntzdce_x] = hitPoint.x();
            zdceventarray[ntzdce_y] = hitPoint.y();
            zdceventarray[ntzdce_z] = hitPoint.z();
            zdceventarray[ntzdce_time] = time;
            zdceventarray[ntzdce_etot] = totalEnergy;
            zdceventntuple->Fill(zdceventarray);
      }
      
      for (std::map<int,float,std::less<int> >::iterator is = energyInFibers.begin();
           is!= energyInFibers.end(); is++)
        {
          ETot = (*is).second;
          SEnergy += ETot;
        }

      // Find Primary info:
      int trackID = 0;
      G4PrimaryParticle* thePrim=0;
      G4int nvertex = (*evt)()->GetNumberOfPrimaryVertex();
      std::cout << "Event has " << nvertex << " vertex" << std::endl;
      if (nvertex==0)
        std::cout << "ZdcTest End Of Event  ERROR: no vertex" << std::endl;
      
      for (int i = 0 ; i<nvertex; i++) {
        
        G4PrimaryVertex* avertex = (*evt)()->GetPrimaryVertex(i);
        if (avertex == 0)
          std::cout << "ZdcTest End Of Event ERR: pointer to vertex = 0"
                       << std::endl;
        std::cout << "Vertex number :" <<i << std::endl;
        int npart = avertex->GetNumberOfParticle();
        if (npart ==0)
          std::cout << "ZdcTest End Of Event ERR: no primary!" << std::endl;
        if (thePrim==0) thePrim=avertex->GetPrimary(trackID);
      }
      
      double px=0.,py=0.,pz=0.;
      double pInit = 0.;
      
      if (thePrim != 0) {
        px = thePrim->GetPx();
        py = thePrim->GetPy();
        pz = thePrim->GetPz();
        pInit = sqrt(pow(px,2.)+pow(py,2.)+pow(pz,2.));
        if (pInit==0) {
          std::cout << "ZdcTest End Of Event  ERR: primary has p=0 " << std::endl;
        }
      } else {
        std::cout << "ZdcTest End Of Event ERR: could not find primary "
                     << std::endl;
      }
      
    } // nentries > 0
    
 
 } // createNTzdcevent
  
  int iEvt = (*evt)()->GetEventID();
  if (iEvt < 10)
    std::cout << " ZdcTest Event " << iEvt << std::endl;
  else if ((iEvt < 100) && (iEvt%10 == 0))
    std::cout << " ZdcTest Event " << iEvt << std::endl;
  else if ((iEvt < 1000) && (iEvt%100 == 0))
    std::cout << " ZdcTest Event " << iEvt << std::endl;
  else if ((iEvt < 10000) && (iEvt%1000 == 0))
    std::cout << " ZdcTest Event " << iEvt << std::endl;
}
void ZdcTestAnalysis::update ( const BeginOfEvent ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfEvent * >.

Definition at line 107 of file ZdcTestAnalysis.cc.

References gather_cfg::cout, eventIndex, and stepIndex.

                                                     {
  //event
  std::cout << "ZdcTest: Processing Event Number: "<<eventIndex<< std::endl;
  eventIndex++;
  stepIndex = 0;
}
void ZdcTestAnalysis::update ( const BeginOfJob ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfJob * >.

Definition at line 77 of file ZdcTestAnalysis.cc.

References gather_cfg::cout.

                                                   {
  //job
  std::cout<<"beggining of job"<<std::endl;;
}
void ZdcTestAnalysis::update ( const BeginOfRun ) [private, virtual]

This routine will be called when the appropriate signal arrives.

Implements Observer< const BeginOfRun * >.

Definition at line 84 of file ZdcTestAnalysis.cc.

References gather_cfg::cout, doNTzdcevent, doNTzdcstep, eventIndex, eventNtFileName, stepNtFileName, zdcOutputEventFile, and zdcOutputStepFile.

                                                   {
  //run

 std::cout << std::endl << "ZdcTestAnalysis: Begining of Run"<< std::endl; 
  if (doNTzdcstep) { 
    std::cout << "ZDCTestAnalysis: output step file created"<< std::endl;
    TString stepfilename = stepNtFileName;
    zdcOutputStepFile = new TFile(stepfilename,"RECREATE");

  }
  
  if (doNTzdcevent) {
    std::cout << "ZDCTestAnalysis: output event file created"<< std::endl;
    TString stepfilename = eventNtFileName;
    zdcOutputEventFile = new TFile(stepfilename,"RECREATE");
  }

  eventIndex = 0;
}

Member Data Documentation

Definition at line 92 of file ZdcTestAnalysis.h.

Referenced by finish(), update(), and ZdcTestAnalysis().

Definition at line 91 of file ZdcTestAnalysis.h.

Referenced by finish(), update(), and ZdcTestAnalysis().

Definition at line 102 of file ZdcTestAnalysis.h.

Referenced by finish(), and update().

std::string ZdcTestAnalysis::eventNtFileName [private]

Definition at line 94 of file ZdcTestAnalysis.h.

Referenced by update(), and ZdcTestAnalysis().

Definition at line 103 of file ZdcTestAnalysis.h.

Referenced by update().

std::string ZdcTestAnalysis::stepNtFileName [private]

Definition at line 93 of file ZdcTestAnalysis.h.

Referenced by update(), and ZdcTestAnalysis().

Definition at line 90 of file ZdcTestAnalysis.h.

Referenced by update(), ZdcTestAnalysis(), and ~ZdcTestAnalysis().

Float_t ZdcTestAnalysis::zdceventarray[16] [private]

Definition at line 106 of file ZdcTestAnalysis.h.

Referenced by update().

TNtuple* ZdcTestAnalysis::zdceventntuple [private]

Definition at line 100 of file ZdcTestAnalysis.h.

Referenced by finish(), update(), and ZdcTestAnalysis().

Definition at line 96 of file ZdcTestAnalysis.h.

Referenced by finish(), and update().

Definition at line 97 of file ZdcTestAnalysis.h.

Referenced by finish(), and update().

Float_t ZdcTestAnalysis::zdcsteparray[18] [private]

Definition at line 105 of file ZdcTestAnalysis.h.

Referenced by update().

TNtuple* ZdcTestAnalysis::zdcstepntuple [private]

Definition at line 99 of file ZdcTestAnalysis.h.

Referenced by finish(), update(), and ZdcTestAnalysis().