Merge pull request #364 from slava77/from-pr363/bs-standalone

write/read beam spot in a standalone setup

Event.cc

@@ -56,6 +56,7 @@ void Event::Reset(int evtID)
   fitTracksExtra_.clear();
   cmsswTracks_.clear();
   cmsswTracksExtra_.clear();
+  beamSpot_ = {};
 
   validation_.resetValidationMaps(); // need to reset maps for every event.
 }
@@ -172,6 +173,11 @@ void Event::write_out(DataFile &data_file)
     evsize += write_tracks(fp, cmsswTracks_);
   }
 
+  if (data_file.HasBeamSpot()) {
+    fwrite(&beamSpot_, sizeof(BeamSpot), 1, fp);
+    evsize += sizeof(BeamSpot);
+  }
+
   fseek(fp, start, SEEK_SET);
   fwrite(&evsize, sizeof(int), 1, fp);
   fseek(fp, 0, SEEK_END);
@@ -367,6 +373,11 @@ void Event::read_in(DataFile &data_file, FILE *in_fp)
   }
 #endif
 
+  if (data_file.HasBeamSpot())
+  {
+    fread(&beamSpot_, sizeof(BeamSpot), 1, fp);
+  }
+
   if (Config::kludgeCmsHitErrors)
   {
     kludge_cms_hit_errors();
@@ -413,7 +424,7 @@ int Event::read_tracks(FILE *fp, TrackVec& tracks, bool skip_reading)
 
   if (skip_reading)
   {
-    fseek(fp, data_size, SEEK_CUR);
+    fseek(fp, data_size-2*sizeof(int), SEEK_CUR);// -2 because data_size counts itself and n_tracks too
     n_tracks = -n_tracks;
   }
   else
@@ -881,7 +892,7 @@ void Event::fill_hitmask_bool_vectors(std::vector<int> &track_algo_vec,
 int DataFile::OpenRead(const std::string& fname, bool set_n_layers)
 {
   constexpr int min_ver = 4;
-  constexpr int max_ver = 5;
+  constexpr int max_ver = 6;
 
   f_fp = fopen(fname.c_str(), "r");
   assert (f_fp != 0 && "Opening of input file failed.");

Event.h

@@ -74,7 +74,7 @@ typedef std::vector<Event> EventVec;
 struct DataFileHeader
 {
   int f_magic          = 0xBEEF;
-  int f_format_version = 5;
+  int f_format_version = 6;
   int f_sizeof_track   = sizeof(Track);
   int f_sizeof_hit     = sizeof(Hit);
   int f_sizeof_hot     = sizeof(HitOnTrack);
@@ -96,7 +96,8 @@ struct DataFile
     ES_SimTrackStates = 0x1,
     ES_Seeds          = 0x2,
     ES_CmsswTracks    = 0x4,
-    ES_HitIterMasks   = 0x8
+    ES_HitIterMasks   = 0x8,
+    ES_BeamSpot       = 0x10
   };
 
   FILE *f_fp  =  0;
@@ -112,6 +113,7 @@ struct DataFile
   bool HasSeeds()          const { return f_header.f_extra_sections & ES_Seeds; }
   bool HasCmsswTracks()    const { return f_header.f_extra_sections & ES_CmsswTracks; }
   bool HasHitIterMasks()   const { return f_header.f_extra_sections & ES_HitIterMasks; }
+  bool HasBeamSpot()       const { return f_header.f_extra_sections & ES_BeamSpot; }
 
   int  OpenRead (const std::string& fname, bool set_n_layers = false);
   void OpenWrite(const std::string& fname, int nev, int extra_sections=0);

tkNtuple/WriteMemoryFile.cc

@@ -504,13 +504,29 @@ int main(int argc, char *argv[])
   vector<vector<float> >*    str_chargeFraction = 0;
   t->SetBranchAddress("str_chargeFraction", &str_chargeFraction);
 
+  // beam spot
+  float bsp_x;
+  float bsp_y;
+  float bsp_z;
+  float bsp_sigmax;
+  float bsp_sigmay;
+  float bsp_sigmaz;
+  t->SetBranchAddress("bsp_x", &bsp_x);
+  t->SetBranchAddress("bsp_y", &bsp_y);
+  t->SetBranchAddress("bsp_z", &bsp_z);
+  t->SetBranchAddress("bsp_sigmax", &bsp_sigmax);
+  t->SetBranchAddress("bsp_sigmay", &bsp_sigmay);
+  t->SetBranchAddress("bsp_sigmaz", &bsp_sigmaz);
+
   long long totentries = t->GetEntries();
   long long savedEvents = 0;
 
   DataFile data_file;
   int outOptions = DataFile::ES_Seeds;
   if (writeRecTracks) outOptions |= DataFile::ES_CmsswTracks;
   if (writeHitIterMasks) outOptions |= DataFile::ES_HitIterMasks;
+  outOptions |= DataFile::ES_BeamSpot;
+
   if (maxevt < 0) maxevt = totentries;
   data_file.OpenWrite(outputFileName, std::min(maxevt, totentries), outOptions);
 
@@ -530,6 +546,15 @@ int main(int argc, char *argv[])
 
     cout << "edm event=" << event << endl;
 
+    auto &bs = EE.beamSpot_;
+    bs.x = bsp_x;
+    bs.y = bsp_y;
+    bs.z = bsp_z;
+    bs.sigmaZ = bsp_sigmaz;
+    bs.beamWidthX = bsp_sigmax;
+    bs.beamWidthY = bsp_sigmay;
+    //dxdz and dydz are not in the trackingNtuple at the moment
+
     for (unsigned int istr = 0; istr < str_lay->size(); ++istr) {
       if(str_chargePerCM->at(istr) < cutValueCCC) numFailCCC++;
       numTotalStr++;
@@ -1119,17 +1144,18 @@ int main(int argc, char *argv[])
 	}
 
 	for (int i=0;i<ns;++i) {
-	  printf("seed id=%i label=%i q=%2i pT=%6.3f p=(%6.3f, %6.3f, %6.3f) x=(%6.3f, %6.3f, %6.3f)\n",i,
-		 seedTracks_[i].label(),seedTracks_[i].charge(),seedTracks_[i].pT(),seedTracks_[i].px(),seedTracks_[i].py(),seedTracks_[i].pz(),seedTracks_[i].x(),seedTracks_[i].y(),seedTracks_[i].z());
+	  printf("seed id=%i label=%i algo=%i q=%2i pT=%6.3f p=(%6.3f, %6.3f, %6.3f) x=(%6.3f, %6.3f, %6.3f)\n",i,
+		 seedTracks_[i].label(),seedTracks_[i].algorithm(),seedTracks_[i].charge(),
+                 seedTracks_[i].pT(),seedTracks_[i].px(),seedTracks_[i].py(),seedTracks_[i].pz(),seedTracks_[i].x(),seedTracks_[i].y(),seedTracks_[i].z());
 	  int nh = seedTracks_[i].nTotalHits();
 	  for (int ih=0;ih<nh;++ih) printf("seed #%i hit #%i idx=%i\n",i,ih,seedTracks_[i].getHitIdx(ih));
 	}
 
 	if (writeRecTracks){
 	  for (int i=0;i<nr;++i) {
 	    float spt = sqrt(pow(cmsswTracks_[i].px(),2)+pow(cmsswTracks_[i].py(),2));
-	    printf("rec track id=%i label%i chi2=%6.3f q=%2i p=(%6.3f, %6.3f, %6.3f) x=(%6.3f, %6.3f, %6.3f) pT=%7.4f nTotal=%i nFound=%i \n",
-		   i, cmsswTracks_[i].label(), cmsswTracks_[i].chi2(),
+	    printf("rec track id=%i label=%i algo=%i chi2=%6.3f q=%2i p=(%6.3f, %6.3f, %6.3f) x=(%6.3f, %6.3f, %6.3f) pT=%7.4f nTotal=%i nFound=%i \n",
+		   i, cmsswTracks_[i].label(), cmsswTracks_[i].algorithm(), cmsswTracks_[i].chi2(),
 		   cmsswTracks_[i].charge(),cmsswTracks_[i].px(),cmsswTracks_[i].py(),cmsswTracks_[i].pz(),cmsswTracks_[i].x(),cmsswTracks_[i].y(),cmsswTracks_[i].z(),spt,
 		   cmsswTracks_[i].nTotalHits(),cmsswTracks_[i].nFoundHits());
 	    int nh = cmsswTracks_[i].nTotalHits();