TestClusterWithHit Class Reference

#include <TestClusterWithHit.h>

Inheritance diagram for TestClusterWithHit:

Public Member Functions
	TestClusterWithHit (const std::string &name, ISvcLocator *pSvcLocator)
	~TestClusterWithHit ()
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
void	reset ()
void	reset_hit ()
void	reset_cluster_1d ()
void	reset_cluster_2d ()
void	reset_map ()

Detailed Description

Definition at line 24 of file TestClusterWithHit.h.

Constructor & Destructor Documentation

TestClusterWithHit()

TestClusterWithHit::TestClusterWithHit	(	const std::string &	name,
		ISvcLocator *	pSvcLocator )

Definition at line 44 of file TestClusterWithHit.cxx.

                                                                                     :
    Algorithm(name,pSvcLocator){   
 
  //  declareProperty("CosmicRayDataSetID",    CosmicRayDataSetID = "CR201909");
  declareProperty("selGoodDigi",m_selGoodDigi=1);
  declareProperty("minDigiTime",m_minDigiTime=-8875);
  declareProperty("maxDigiTime",m_maxDigiTime=-8562);
  declareProperty("minQDigi",myQMin=0.0);
  declareProperty("LUTfile", lutfile = "/bes3fs/cgemCosmic/data/CGEM_cosmic_look_up_table_from_17_to_17.root");
 
 
  
}

~TestClusterWithHit()

TestClusterWithHit::~TestClusterWithHit

(

)

Definition at line 58 of file TestClusterWithHit.cxx.

                                       {
  // delete m_SvcCgem;
  // delete f;
  // delete Tdigi;  
  // delete anode;
  // delete tree;
  // delete output;
  delete lutreader;  
}

Member Function Documentation

execute()

StatusCode TestClusterWithHit::execute

(

)

read map std::map<int, std::vector < int > >::iterator it2; for(it2 = map_L1_S1_stripx_to_hit.begin(); it2 != map_L1_S1_stripx_to_hit.end(); it2++) { int istrip = it2->first; std::vector< int > hitlist = it2->second; cout << "L1/S1/x....reading strip " << istrip << " hit " << hitlist[0] << endl;
}

cout << "1D_CLUSTER index " << clusterid << " or " << ncluster << " layer " << layerid << " sheet " << sheetid << " view " << cluster_1d_view[ncluster] << endl; cout << " from " << flagB << " to " << flagE << endl; cout << "hits: ";

cout << "2D_CLUSTER index " << clusterid << " or " << ncluster << " layer " << layerid << " sheet " << sheetid << " view " << cluster_2d_view[ncluster] << endl; cout << " idx " << flagB << " idv " << flagE << endl;

cout << "************************" << endl; for(int iclu=0; iclu<ncluster; iclu++) { cout << "cluster 1d view " << cluster_1d_view[iclu] << " size " << cluster_1d_size[iclu] << " from " << cluster_1d_strip1[iclu] << " to " << cluster_1d_strip2[iclu] << endl;

for(int ihit = 0; ihit < cluster_1d_size[iclu] ; ihit++) { cout << "hit " << cluster_1d_hitindex[iclu][ihit] << endl; } }

Definition at line 303 of file TestClusterWithHit.cxx.

                                      {
  
  MsgStream log(msgSvc(), name());
  if(event%1000==0) cout << "->TestClusterWithHit::execute " << event << endl;
 
  //interface to event data service
  ISvcLocator* svcLocator = Gaudi::svcLocator();
  StatusCode sc=svcLocator->service("EventDataSvc", m_evtSvc);
  if (sc.isFailure())
    cout<<"Could not accesss EventDataSvc!"<<endl;
 
  // reset
  reset();
  // -------------------------
 
  //retrieve CGEM hits from TDS
  SmartDataPtr<CgemDigiCol> aDigiCol(m_evtSvc,"/Event/Digi/CgemDigiCol");
  if(!aDigiCol)
    cout<<"Could not retrieve CGEM digi collection"<<endl;
  nhit = aDigiCol->size();
  
  // loop on hits
  int ihit=0;
  
  CgemDigiCol::iterator iDigiCol;
  for(iDigiCol=aDigiCol->begin(); iDigiCol!=aDigiCol->end(); iDigiCol++)
    {
      const Identifier ident = (*iDigiCol)->identify();
      int strip = CgemID::strip(ident);
      int layer = CgemID::layer(ident);
      int sheet = CgemID::sheet(ident);
      int view = 1;
      bool is_xstrip = CgemID::is_xstrip(ident);
      if(is_xstrip == true) view = 0;
      double charge = (*iDigiCol)->getCharge_fc();
      double time = (*iDigiCol)->getTime_ns();
 
      hit_strip[ihit] = strip;
      hit_layer[ihit] = layer;
      hit_sheet[ihit] = sheet;
      hit_view[ihit] = view;
 
      if(selDigi(iDigiCol, m_selGoodDigi)==true) {
 
    //  cout << "HIT " << ihit << " strip " << strip << " layer " << layer << " sheet " << sheet << " view " << view << endl;
    // map strips to hits (maybe more hits for the same strip)
    std::map<int, std::vector < int > > *map_strip_to_hit = GetStripToHitMap(layer, sheet, view);
    std::map<int, std::vector < int > >::iterator it1 = map_strip_to_hit->find(strip);
    if (it1 != map_strip_to_hit->end()) {
      it1->second.push_back(ihit); 
    }
    else {
      std::vector< int > hitlist;
      hitlist.push_back(ihit);
      std::pair<int, std::vector < int> > pair(strip, hitlist);
      map_strip_to_hit->insert(pair);
    }
      }
      // --------------------------------------------------------
 
      ihit++;
    }
  /**
  //read map
  std::map<int, std::vector < int > >::iterator it2;
  for(it2 = map_L1_S1_stripx_to_hit.begin(); it2 != map_L1_S1_stripx_to_hit.end(); it2++)
    {
      int istrip = it2->first;
      std::vector< int > hitlist = it2->second;
      cout << "L1/S1/x....reading strip " << istrip << " hit " << hitlist[0] << endl;  
  }
    **/
 
  // -------------------------------------------------------------------  
  //retrieve CGEM clusters from TDS
  SmartDataPtr<RecCgemClusterCol> aClusterCol(m_evtSvc,"/Event/Recon/RecCgemClusterCol");
  if(!aClusterCol)
    cout<<"Could not retrieve CGEM cluster collection"<<endl;
  
  int nclu = aClusterCol->size();
  if(nclu > MAXNOFCLUSTERS) {
    std::cout << "nclu " << nclu << std::endl;
    nclu = MAXNOFCLUSTERS; // CHECK HACK
  }
  
  int nclusterx = 0;
  int nclusterv = 0;
  int nclusterxv = 0;
  
  anode = m_SvcCgem->getReadoutPlane(0, 0);
  anode_radius_L1_x = anode->getRX();
  anode_radius_L1_v = anode->getRV();
  anode_mid_gap_L1 = anode->getMidRAtGap();
  anode = m_SvcCgem->getReadoutPlane(1, 0);
  anode_radius_L2_x = anode->getRX();
  anode_radius_L2_v = anode->getRV();
  anode_mid_gap_L2 = anode->getMidRAtGap();
 
  int tmp_cluster_L1_S1 = -1;
  int tmp_cluster_L2_S1 = -1;
  int tmp_cluster_L2_S2 = -1;
  double tmp_charge_L1_S1 = 0.;
  double tmp_charge_L2_S1 = 0.;
  double tmp_charge_L2_S2 = 0.;
 
  ncluster = 0;
  RecCgemClusterCol::iterator iClusterCol;
  for(iClusterCol=aClusterCol->begin(); iClusterCol!=aClusterCol->end(); iClusterCol++)
     {
 
    if(ncluster==MAXNOFCLUSTERS) break;
    
        int flag = (*iClusterCol)->getflag(); // 0=x 1=v 2=xv
    int clusterid = (*iClusterCol)->getclusterid();
    int trkid = (*iClusterCol)->getTrkId();
    int layerid = (*iClusterCol)->getlayerid();
    int sheetid = (*iClusterCol)->getsheetid();
    double edep = (*iClusterCol)->getenergydeposit();
    double phi = (*iClusterCol)->getrecphi();
    double v = (*iClusterCol)->getrecv();
    double z = (*iClusterCol)->getRecZ();
    double phi_cc = (*iClusterCol)->getrecphi_CC();
    double v_cc = (*iClusterCol)->getrecv_CC();
    double z_cc = (*iClusterCol)->getRecZ_CC();
    double phi_tpc = (*iClusterCol)->getrecphi_mTPC();
    double v_tpc = (*iClusterCol)->getrecv_mTPC();
    double z_tpc = (*iClusterCol)->getRecZ_mTPC();
        double a_tpc = 0; //(*iClusterCol)->getparamA_mTPC();
        double b_tpc = 0; //(*iClusterCol)->getparamB_mTPC();
 
    if(sheetid == -1)  {
      ncluster++;
      continue; // CHECK HACK FOR NOW
    }   
     
    // anode mid gap radius
    anode = m_SvcCgem->getReadoutPlane(layerid, sheetid);
    double anode_mid_gap = anode->getMidRAtGap();
 
        // beginning and ending flags 
        // if 2D cluster --> they are the x & v 1D clusters
        // if 1D cluster --> they are the first and last contigous strips 
        int flagB = (*iClusterCol)->getclusterflagb();
        int flagE = (*iClusterCol)->getclusterflage();
 
        // ------------------------------------------------------------------- 1D cluster
    // positions CC & uTPC
    if(flag==0 || flag == 1) {
          cluster_1d_ID[ncluster] = clusterid;
      // time & charge
      cluster_1d_t[ncluster] = 0; // CHECK 
      cluster_1d_q[ncluster] = edep;
      // layer/sheet/view
      cluster_1d_layerid[ncluster] = layerid;
      cluster_1d_sheetid[ncluster] = sheetid;
      cluster_1d_view[ncluster] = flag;
      cluster_1d_r[ncluster] = anode_mid_gap;
      cluster_1d_phi[ncluster] = phi;
      cluster_1d_phi_cc[ncluster] = phi_cc;
      cluster_1d_phi_tpc[ncluster] = phi_tpc;
      cluster_1d_v[ncluster] = v;
      cluster_1d_v_cc[ncluster] = v_cc;
          cluster_1d_v_tpc[ncluster] = v_tpc;
          cluster_1d_a_tpc[ncluster] = a_tpc;
          cluster_1d_b_tpc[ncluster] = b_tpc;
      // from strip1 to strip2
          cluster_1d_strip1[ncluster] = flagB;
          cluster_1d_strip2[ncluster] = flagE;
      int size = flagE - flagB + 1;
          cluster_1d_size[ncluster] = size;
      
      //      cout << "getting the hits from strip1 "<< flagB << " to strip2 " << flagE << " size " << size << endl;
          // get the hits
      std::map<int, std::vector < int > > *read_map_strip_to_hit = GetStripToHitMap(layerid, sheetid, flag);
      std::map<int, std::vector < int > >::iterator it2;
      /**
         cout << "1D_CLUSTER index " << clusterid << " or " << ncluster
         << " layer " << layerid << " sheet " << sheetid
         << " view " << cluster_1d_view[ncluster] << endl;
         cout << " from " << flagB << " to " << flagE << endl;
         cout << "hits: ";
      **/        
 
      //      if(size >= 50)      cout << "size " << size << endl;      
      for(int istrip = 0; istrip < size; istrip++) {
            int stripid = flagB + istrip;
            it2 = read_map_strip_to_hit->find(stripid);
        std::vector< int > hitlist = it2->second;
        int size_hitlist = hitlist.size();
        // if(size_hitlist>1)  cout << event << " cluster " << ncluster << " hitlist > 1 " << size_hitlist << endl;
        cluster_1d_hitindex[ncluster][istrip] = hitlist[size_hitlist-1]; // take only the first hit CHECK THIS HAS TO BE CHANGED!
        // cout << " " << hitlist[0];
      }
      //      cout << endl;
      // -----------
      ncluster_1d++;
    }
    else if(flag==2 || flag==3) { // ------------------------------------- 2D - xv view
      cluster_2d_ID[ncluster] = clusterid;
      // time & charge 
          cluster_2d_t[ncluster] = 0; // CHECK 
          cluster_2d_q[ncluster] = edep;
          // layer/sheet/view 
          cluster_2d_layerid[ncluster] = layerid;
      cluster_2d_sheetid[ncluster] = sheetid;
          cluster_2d_view[ncluster] = flag;
      cluster_2d_r[ncluster] = anode_mid_gap;
      cluster_2d_z[ncluster] = z;
      cluster_2d_z_cc[ncluster] = z_cc;
          cluster_2d_z_tpc[ncluster] = z_tpc;
      cluster_2d_phi[ncluster] = phi;
      cluster_2d_phi_cc[ncluster] = phi_cc;
          cluster_2d_phi_tpc[ncluster] = phi_tpc;
      // composing 1D_clusters
      cluster_2d_idx[ncluster] = flagB;
          cluster_2d_idv[ncluster] = flagE;
              
      /**     cout << "2D_CLUSTER index " << clusterid << " or " << ncluster
               << " layer " << layerid << " sheet " << sheetid
               << " view " << cluster_2d_view[ncluster] << endl;
          cout << " idx " << flagB << " idv " << flagE << endl;
       
      **/
      // find the highest charge cluster
      if(flag==2) { // CHECK THIS
        if(layerid==0) {
          if(edep > tmp_charge_L1_S1) {
        tmp_charge_L1_S1 = edep; tmp_cluster_L1_S1 = ncluster; 
          }
        }
        else if(layerid==1 && sheetid==0) {
          if(edep > tmp_charge_L2_S1){
        tmp_charge_L2_S1 = edep; tmp_cluster_L2_S1 = ncluster; 
          }
        }
        else if(layerid==1 && sheetid==1) {
          if(edep > tmp_charge_L2_S2){
        tmp_charge_L2_S2 = edep; tmp_cluster_L2_S2 = ncluster; 
          }
        }
          }
      // -----------
      ncluster_2d++;
    }
 
        // cluster counters
    ncluster++;
    if(flag == 0) {
      nclusterx++;
      if(layerid == 0) ncluster_1d_L1_S1_x++;
      else if(layerid == 1) {
        if(sheetid == 0) ncluster_1d_L2_S1_x++;
        else             ncluster_1d_L2_S2_x++;
      }
    }
    else if(flag == 1) {
          nclusterv++;
      if(layerid == 0) ncluster_1d_L1_S1_v++;
      else if(layerid == 1) {
        if(sheetid == 0) ncluster_1d_L2_S1_v++;
        else             ncluster_1d_L2_S2_v++;
      }
    }
        else if(flag == 2) {
      nclusterxv++;
      if(layerid == 0) ncluster_2d_L1_S1++;
          else if(layerid == 1) {
            if(sheetid == 0) ncluster_2d_L2_S1++;
            else             ncluster_2d_L2_S2++;
          }
    }
    // --------
     }
 
   //    std::cout << "cluster " << tmp_cluster_L1_S1 << " " << tmp_cluster_L2_S1 << " " << tmp_cluster_L2_S2 << std::endl;
   //    std::cout << "charge " << tmp_charge_L1_S1 << " " << tmp_charge_L2_S1 << " " << tmp_charge_L2_S2 << std::endl;
 
  // set the highest charge 2D - cluster
  if(tmp_cluster_L1_S1!=-1)      cluster_2d_highest[tmp_cluster_L1_S1]=1;
  if(tmp_cluster_L2_S1!=-1)      cluster_2d_highest[tmp_cluster_L2_S1]=1;
  if(tmp_cluster_L2_S2!=-1)      cluster_2d_highest[tmp_cluster_L2_S2]=1;
  
  /**
  cout << "************************" << endl;
  for(int iclu=0; iclu<ncluster; iclu++) {
    cout << "cluster 1d view " << cluster_1d_view[iclu] << " size " << cluster_1d_size[iclu] << " from " << cluster_1d_strip1[iclu] << " to " << cluster_1d_strip2[iclu] << endl;
 
    for(int ihit = 0; ihit <  cluster_1d_size[iclu] ; ihit++) {
      cout << "hit " << cluster_1d_hitindex[iclu][ihit] << endl;
    }
  }
  **/
 
  tree->Fill();
  event++;
 
  return StatusCode::SUCCESS;
}

finalize()

StatusCode TestClusterWithHit::finalize

(

)

Definition at line 602 of file TestClusterWithHit.cxx.

                                       {
    MsgStream log(msgSvc(),name());
    log << MSG::INFO << "TestClusterWithHit finalize()" << endreq;
    output->Write();
    output->Close();
    return StatusCode::SUCCESS;
}

initialize()

StatusCode TestClusterWithHit::initialize

(

)

Definition at line 68 of file TestClusterWithHit.cxx.

                                         {
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "TestClusterWithHit initialize()" << endreq;
 
 
  // CgemGeomSvc                                                                                                                       
  StatusCode sc;       
  sc = service("CgemGeomSvc", m_SvcCgem);                                                                                        
  if(sc != StatusCode::SUCCESS)  { 
    log << MSG::ERROR << "can not use CgemGeomSvc" << endreq;                                                                   
    return StatusCode::FAILURE;                                                                                                 
  }  
  
  // LUT
  lutreader = new CgemLUTReader(lutfile);
  lutreader->ReadLUT();
 
  output = new TFile("cluster.root", "RECREATE");
  tree = new TTree("tree","cluster info tree");
  DefineClusterTree();
  DefineHitTree();
 
  event=0;
 
  return StatusCode::SUCCESS;
 
}

reset()

void TestClusterWithHit::reset

(

)

Definition at line 159 of file TestClusterWithHit.cxx.

                               {
  reset_map();
  reset_cluster_1d();
  reset_cluster_2d();
  reset_hit();
 
}

Referenced by execute().

reset_cluster_1d()

void TestClusterWithHit::reset_cluster_1d

(

)

Definition at line 178 of file TestClusterWithHit.cxx.

                                          {
 
  ncluster_1d = 0;
  ncluster_1d_L1_S1_x = 0;
  ncluster_1d_L2_S1_x = 0;
  ncluster_1d_L2_S2_x = 0;
  ncluster_1d_L1_S1_v = 0;
  ncluster_1d_L2_S1_v = 0;
  ncluster_1d_L2_S2_v = 0;
 
  for(int iclu = 0; iclu < MAXNOFCLUSTERS; iclu++)
    {
      cluster_1d_t[iclu] = 0.;
      cluster_1d_q[iclu] = 0.;
      cluster_1d_r[iclu] = 0.;
      cluster_1d_v[iclu] = 0.;
      cluster_1d_phi[iclu] = 0.;
      cluster_1d_v_cc[iclu] = 0.;
      cluster_1d_phi_cc[iclu] = 0.;
      cluster_1d_v_tpc[iclu] = 0.;
      cluster_1d_phi_tpc[iclu] = 0.;
      cluster_1d_a_tpc[iclu] = 0.;
      cluster_1d_b_tpc[iclu] = 0.;
      cluster_1d_layerid[iclu] = -1;
      cluster_1d_sheetid[iclu] = -1;
      cluster_1d_view[iclu] = -1;
      cluster_1d_size[iclu] = -1;
      cluster_1d_strip1[iclu] = -1;
      cluster_1d_strip2[iclu] = -1;
      for(int ihit = 0; ihit < MAXNOFHITS; ihit++) cluster_1d_hitindex[iclu][ihit] = -1;
    }
 
 
 
}

Referenced by reset().

reset_cluster_2d()

void TestClusterWithHit::reset_cluster_2d

(

)

Definition at line 214 of file TestClusterWithHit.cxx.

                                          {
 
  ncluster_2d = 0;
  ncluster_2d_L1_S1 = 0;
  ncluster_2d_L2_S1 = 0;
  ncluster_2d_L2_S2 = 0;
 
  for(int iclu = 0; iclu < MAXNOFCLUSTERS; iclu++)
    {
      cluster_2d_t[iclu] = 0.;
      cluster_2d_q[iclu] = 0.;
      cluster_2d_r[iclu] = 0.;
      cluster_2d_z[iclu] = 0.;
      cluster_2d_phi[iclu] = 0.;
      cluster_2d_z_cc[iclu] = 0.;
      cluster_2d_phi_cc[iclu] = 0.;
      cluster_2d_z_tpc[iclu] = 0.;
      cluster_2d_phi_tpc[iclu] = 0.;
      cluster_2d_layerid[iclu] = -1;
      cluster_2d_sheetid[iclu] = -1;
      cluster_2d_view[iclu] = -1;
      cluster_2d_idx[iclu] = -1;
      cluster_2d_idv[iclu] = -1;
      cluster_2d_highest[iclu] = -1;
    }
}

Referenced by reset().

reset_hit()

void TestClusterWithHit::reset_hit

(

)

Definition at line 166 of file TestClusterWithHit.cxx.

                                   {
  // event = 0;
  nhit = 0;
  
  for(int ihit = 0; ihit < MAXNOFHITS; ihit++) {
    hit_strip[ihit] = -1;
    hit_view[ihit] = -1;
    hit_layer[ihit] = -1;
    hit_sheet[ihit] = -1;
  }
}

Referenced by reset().

reset_map()

void TestClusterWithHit::reset_map

(

)

Definition at line 243 of file TestClusterWithHit.cxx.

                                   {
  map_L1_S1_stripx_to_hit.clear();
  map_L2_S1_stripx_to_hit.clear();
  map_L2_S2_stripx_to_hit.clear();
  map_L1_S1_stripv_to_hit.clear();
  map_L2_S1_stripv_to_hit.clear();
  map_L2_S2_stripv_to_hit.clear();
}

Referenced by reset().

---

The documentation for this class was generated from the following files:

• TestClusterWithHit.h
• TestClusterWithHit.cxx