online_coinc.cpp 12.7 KB
Newer Older
1
// nptool headers
gypaos's avatar
gypaos committed
2 3 4
#include "NPOptionManager.h"
#include "RootOutput.h"
#include "NPDetectorManager.h"
5 6
#include "TComptonTelescopeData.h"
#include "TComptonTelescopePhysics.h"
7 8

// root headers
LAVIRON Adrien's avatar
LAVIRON Adrien committed
9
#include "TCutG.h"
10
//#include "TH2.h"
LAVIRON Adrien's avatar
LAVIRON Adrien committed
11 12
#include "TFile.h"
#include "TCanvas.h"
13 14

// custom headers
15
#include "DecodeR.h"
16
#include "DecodeD.h"
17 18
#include "DecodeT.h"

19 20 21 22
#define __RUN__ 3

#define __1DET__

23
#define __TEST_ZONE__
24
#undef __TEST_ZONE__
25 26 27

#define __CIRCULAR_TREE__
#undef __CIRCULAR_TREE__
LAVIRON Adrien's avatar
LAVIRON Adrien committed
28

LAVIRON Adrien's avatar
LAVIRON Adrien committed
29
#define __USE_CUTG__
30
#undef __USE_CUTG__
31 32 33 34

// C++ headers
#include <iostream>
#include <fstream>
gypaos's avatar
gypaos committed
35
#include <string>
36 37
using namespace std;

38 39

//--//--// One-line setter for DSSSD(s) //--//--//
40 41 42 43 44 45 46 47 48 49 50 51 52
void setCTTracker(TComptonTelescopeData* ccamData, DecodeD* DD)
{
  for (int i = 0; i < DD->getEventSize(); i++) {
    if (DD -> getFaceType(i) == 0) { // front
      ccamData->SetFrontE(1, DD->getDetNbr(i)+1, DD->getStripNbr(i), DD->getEnergy(i));
      ccamData->SetFrontT(1, DD->getDetNbr(i)+1, 33, DD->getTime());
    }
    else if (DD -> getFaceType(i) == 1) { // back
      ccamData->SetBackE(1, DD->getDetNbr(i)+1, DD->getStripNbr(i), DD->getEnergy(i));
      ccamData->SetBackT(1, DD->getDetNbr(i)+1, 33, DD->getTime());
    }
  }
}
53 54

//--//--//--//--//--//--//--//--//--//--//--//--//--//--//--//
55
int main(int argc, char** argv)
56
{
LAVIRON Adrien's avatar
LAVIRON Adrien committed
57
  int resetCountSearch = 3;
58
  int timestampDiffSearch = 1000;
LAVIRON Adrien's avatar
LAVIRON Adrien committed
59
  int timestampNBins = 100;
60 61 62 63
  Option_t* tfoption;
  if (argc > 1 and std::string(argv[1]) == "-r") {tfoption = "recreate";}
  else {tfoption = "update";}
  auto fout = new TFile("timesearch.root", tfoption);
LAVIRON Adrien's avatar
LAVIRON Adrien committed
64
  auto bidim = new TH2F("bidim", "bidim", 
65 66 67 68 69
#ifdef __TEST_ZONE__
      150, 0, 150000,
#else
      2*resetCountSearch+1, -resetCountSearch-.5, resetCountSearch+.5, 
#endif
LAVIRON Adrien's avatar
LAVIRON Adrien committed
70
      timestampNBins, -timestampDiffSearch, timestampDiffSearch);
71 72
  auto deltaT = new TH1F("deltaT", "deltaT", timestampNBins, -timestampDiffSearch, timestampDiffSearch);

LAVIRON Adrien's avatar
LAVIRON Adrien committed
73
#ifdef __USE_CUTG__
74 75
  //TFile* fcut = new TFile("/disk/proto-data/data/CUT_Compton.root");
  TFile* fcut = new TFile("../data/CUT_Compton.root");
LAVIRON Adrien's avatar
LAVIRON Adrien committed
76 77 78 79 80
  TCutG* mcut = (TCutG*) fcut -> Get("CUT_Compton");
  fcut -> Close();
  cout << fcut << endl;
  cout << mcut << endl;
#endif
81

82 83
  ///////////////////////////////////////////////////////////////////////////
  // configure option manager
gypaos's avatar
gypaos committed
84
//   NPOptionManager::getInstance()->Destroy();
85
#ifdef __CIRCULAR_TREE__
LAVIRON Adrien's avatar
LAVIRON Adrien committed
86
  string arg = "-D ./ComptonCAM.detector -C calibrations.txt -GH -E ./10He.reaction --circular";
LAVIRON Adrien's avatar
LAVIRON Adrien committed
87 88 89
#else
  string arg = "-D ./ComptonCAM.detector -C calibrations.txt -GH -E ./10He.reaction";
#endif
90
  NPOptionManager::getInstance(arg);
91

gypaos's avatar
gypaos committed
92
  // open ROOT output file
93
  RootOutput::getInstance("OnlineTree.root", "OnlineTree");
gypaos's avatar
gypaos committed
94 95
  // get tree pointer
  auto m_OutputTree = RootOutput::getInstance()->GetTree();
gypaos's avatar
gypaos committed
96

97 98 99 100
  // configure detector manager
  string detectorfileName = NPOptionManager::getInstance()->GetDetectorFile();
  NPL::DetectorManager* m_NPDetectorManager = new NPL::DetectorManager();
  m_NPDetectorManager->ReadConfigurationFile(detectorfileName);
101
  m_NPDetectorManager->InitializeRootOutput();
102

103 104 105
  // configure spectra server
  m_NPDetectorManager->SetSpectraServer();

106
  // instantiate raw ComptonCAM data pointer
gypaos's avatar
gypaos committed
107
  auto ccamData = new TComptonTelescopeData();
108
  // connect raw CCAM data pointer to physics class
gypaos's avatar
gypaos committed
109 110
  auto ccamPhys = (TComptonTelescopePhysics*) m_NPDetectorManager->GetDetector("ComptonTelescope");
  ccamPhys->SetRawDataPointer(ccamData);
111

112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
#ifdef __1DET__
  ifstream is;
  is.open("/disk/proto-data/data/20210510_Bi207/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-05-10_07_41_50.raw");
  is.seekg(0, ios::end);
  int len = is.tellg();
  is.seekg(0, ios::beg);
  char* buff = new char[len];
  is.read(buff, len);
  is.close();
  DecodeR* DR = new DecodeR(false, buff);
  while(DR -> getCursor() < len) {
    DR -> decodeBlobMFM();
    m_NPDetectorManager->ClearEventPhysics();
    m_NPDetectorManager->ClearEventData();
    ccamData -> SetCTCalorimeter(1, 4, DR->getPixelNumber(), DR->getTime(), DR->getData(), 64);
    m_NPDetectorManager->BuildPhysicalEvent();
    m_OutputTree->Fill();
    m_NPDetectorManager->CheckSpectraServer();
  }
  delete DR;
  delete [] buff;
  m_NPDetectorManager -> WriteSpectra();
#else

136 137
  // read data file/flux and fill ccamData object
  std::cout << "Reading data\n";
138
  DecodeR* DR = new DecodeR(false); // Instantiates DecodeR object reading calorimeter data flux
139 140
  DecodeT* DT = new DecodeT(false); // Instantiates DecodeT object reading trigger data flux
  DecodeD* DD = new DecodeD(false); // Instantiates DecodeD object reading DSSSD(s) data flux
141
//  newframe_t* event;
142 143 144 145 146 147 148
  #if __RUN__ == 0
  DD -> setTree("../data/20210210_run1/bb7_3309-7_cs137-20210210_11h05_coinc_run1_conv.root");
  #elif __RUN__ == 1
  DD -> setTree("/disk/proto-data/data/20210210_run1/bb7_3309-7_cs137-20210210_11h05_coinc_run1_conv.root");
  #elif __RUN__ == 2
  DD -> setTree("/disk/proto-data/data/20210304_run2/bb7_3309-7_cs137_20210304_14h35_conv.root");
  #elif __RUN__ == 3
149
  DD -> setTree("/disk/proto-data/data/20210305_run3/bb7_3309-7_cs137_20210305_14h53_conv.root");
150 151 152
  #elif __RUN__ == 4
  DD -> setTree("/disk/proto-data/data/20210407_run4/bb7_3309-7_cs137_20210407_14h53_conv.root");
  #endif
153
  int dlen = DD -> getLength();
154

155 156
  int i = 0;// ROSMAP files loop counter
  int c = 0;// Event counter
LAVIRON Adrien's avatar
LAVIRON Adrien committed
157
  int cc = 0;
158
  // Set some constants
gypaos's avatar
gypaos committed
159
  const int pixelNumber = 64;
160 161 162

  // Open data files
  ifstream iros, itrig;
163
  cout << "Loading data files " << std::flush;
164

165 166 167 168 169 170
  #if __RUN__ == 0
  itrig.open("../data/20210210_run1/mfm_trigger_202102101104.raw", ios::binary);
  #elif __RUN__ == 1
  itrig.open("/disk/proto-data/data/20210210_run1/mfm_trigger_202102101104.raw", ios::binary);
  #elif __RUN__ == 2
  #elif __RUN__ == 3
171
  itrig.open("/disk/proto-data/data/20210305_run3/mfm_trigger_20210305_run3.raw", ios::binary);
172 173 174
  #elif __RUN__ == 4
  itrig.open("/disk/proto-data/data/20210407_run4/mfm_trigger_20210407_run4.raw", ios::binary);
  #endif
175 176 177 178 179 180
  itrig.seekg(0, ios::end);
  int tlen = itrig.tellg();
  itrig.seekg(0, ios::beg);
  char* tbuff = new char[tlen];
  itrig.read(tbuff, tlen);
  itrig.close();
181
  cout << "... " << std::flush;
182

183 184 185 186 187 188
  #if __RUN__ == 0
  iros.open("../data/20210210_run1/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-02-10_10_04_59.raw.0001", ios::binary);
  #elif __RUN__ == 1
  iros.open("/disk/proto-data/data/20210210_run1/mfm_rdd_rosmap_04_mfm_rosmap_04_2021-02-10_10_04_59.raw.0001", ios::binary);
  #elif __RUN__ == 2
  #elif __RUN__ == 3
189
  iros.open("/disk/proto-data/data/20210305_run3/mfm_rosmap_20210305_run3.raw", ios::binary);
190 191 192
  #elif __RUN__ == 4
  iros.open("/disk/proto-data/data/20210407_run4/mfm_rosmap_20210407_run4.raw", ios::binary);
  #endif
193 194 195 196 197 198 199 200
  iros.seekg(0, ios::end);
  int rlen = iros.tellg();
  iros.seekg(0, ios::beg);
  char* rbuff = new char[rlen];
  iros.read(rbuff, rlen);
  iros.close();
  cout << "Done" << endl;

LAVIRON Adrien's avatar
LAVIRON Adrien committed
201
  // Search for reset count in trigger data
202
  DT -> setRaw(tbuff);
LAVIRON Adrien's avatar
LAVIRON Adrien committed
203
  DT -> decodeBlobMFM();
LAVIRON Adrien's avatar
LAVIRON Adrien committed
204
  int resetCount = DT -> getResetCount();
LAVIRON Adrien's avatar
LAVIRON Adrien committed
205 206
  //while (not(DT->hasTrigged(2))) {
  while (not(DT->hasTrigged(0))) {
LAVIRON Adrien's avatar
LAVIRON Adrien committed
207 208
    DT -> decodeBlobMFM();
  }
LAVIRON Adrien's avatar
LAVIRON Adrien committed
209
  resetCount = DT->getResetCount() - resetCount;
LAVIRON Adrien's avatar
LAVIRON Adrien committed
210
  resetCount = -resetCount; // T B C !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
LAVIRON Adrien's avatar
LAVIRON Adrien committed
211
  cout << "Found reset count: " << resetCount << endl;
212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227

  if (resetCount == 0) {
    string ans;
    cout << "reset count is 0. Continue ? y/[n]";
    cin >> ans;
    if (ans != "y" and ans != "Y") {
      DT -> setRaw(tbuff);
      DT -> decodeBlobMFM();
      resetCount = DT -> getResetCount();
      while (not(DT->hasTrigged(2))) {
        DT -> decodeBlobMFM();
      }
      resetCount = DT->getResetCount() - resetCount;
      cout << "Found reset count: " << resetCount << endl;
    }
  }
228
  int cr, cd, tr, td;
LAVIRON Adrien's avatar
LAVIRON Adrien committed
229

230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  if (argc > 1 and std::string(argv[1]) == "-r") {
    cout << "Now looking for a potential error in reset number by filling the time bidim" << endl;
  
    // Fill control bidim
    for (int reset=-resetCountSearch; reset<resetCountSearch+1; reset++)
    {
      cr = resetCount+reset;
      cout << "Filling coincidence bidim with reset #" << cr << endl;
      DD -> rewind();
  
      DR -> setRaw(rbuff);
      DR -> decodeBlobMFM();
      DD -> decodeEvent();
  
      cd = 0;
      c = 0;
      i = 0;
      tr = DR -> getTime();
      td = DD -> getTime();
      while(c < 1000)
      {
  //    cout << DR -> getTime() << " " << DD -> getTime() << endl; 
        if (cr == cd) {
          if (abs(td-tr) < timestampDiffSearch) {
            bidim->Fill(reset, td-tr);
            c++;
LAVIRON Adrien's avatar
LAVIRON Adrien committed
257
          }
258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273
          if (td < tr) {
            DD -> decodeEvent();
          } else {
            DR -> decodeBlobMFM();
          }
        } else if (cr < cd) {
          DR -> decodeBlobMFM();
        } else {
          DD -> decodeEvent();
        }
        if (DR -> getTime() < tr) {
          cr++;
        }
        tr = DR -> getTime();
        if (DD -> getTime() < td) {
          cd ++;
LAVIRON Adrien's avatar
LAVIRON Adrien committed
274
        }
275 276 277 278 279 280 281 282
        td = DD -> getTime();
        i++; 
      } // End of event loop
    } // End of for loop
  
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    c = 0;
    cout << "Now filling deltaT histogram with all available events" << endl;
283

284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305
    DR -> setRaw(rbuff);
    DR -> decodeBlobMFM();
    DD -> rewind();
    DD -> decodeEvent();

    cr = resetCount; cd = 0;
    tr = DR -> getTime();
    td = DD -> getTime();

    while(DR -> getCursor() < rlen and DD -> getCursor() < dlen)
    {
      if (cr == cd) {
        if (abs(td-tr) < timestampDiffSearch) {
          deltaT -> Fill(td-tr);
          c++;
        }
        if (td < tr) {
          DD -> decodeEvent();
        } else {
          DR -> decodeBlobMFM();
        }
      } else if (cr < cd) {
306
        DR -> decodeBlobMFM();
307 308
      } else {
        DD -> decodeEvent();
309
      }
310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326
      if (DR -> getTime() < tr) {
        cr++;
      }
      tr = DR -> getTime();
      if (DD -> getTime() < td) {
        cd ++;
      }
      td = DD -> getTime();
      i++;
    } // End of event loop
    cout << "Out of " << i << " events treated, " << c << " coincidences were found and added to the histogram" << endl;


/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
  } else {// non reset search mode

    c = 0;
327
    
328 329
    DR -> setRaw(rbuff);
    DR -> decodeBlobMFM();
330

331
//#ifdef __TEST_ZONE__
332
//  cout << "Entering test zone." << endl;
333 334

//#else
335 336 337 338
    DD -> rewind();
    DD -> decodeEvent();

    cr = resetCount; cd = 0;
339 340
    tr = DR -> getTime();
    td = DD -> getTime();
341

342
    while(DR -> getCursor() < rlen and DD -> getCursor() < dlen and c < 10000)
343 344 345 346
    {
      if (cr == cd) {
  #ifndef __TEST_ZONE__
        if (td-tr > 20 and td-tr < 120) { // That one is the real one
meyanne's avatar
meyanne committed
347
        //if (td-tr > -1000 and td-tr < 1000) { // That one is the real one
348
  #else
349
        if (td-tr > -1000  and td-tr < 0) {
350 351 352 353 354 355 356
  #endif
          //DR -> Dump();
          //DD -> Dump();
          c++;
          // Clear raw and physics data
          m_NPDetectorManager->ClearEventPhysics();
          m_NPDetectorManager->ClearEventData();
357
  
358 359 360 361
          // Fill data
          ccamData -> SetCTCalorimeter(1, 4, DR->getPixelNumber(), DR->getTime(), DR->getData(), pixelNumber);
          setCTTracker(ccamData, DD);// -> getEvent(), &nb_asic, &chain, stripNumber);
          ccamData -> SetResetCount(cr);
362
  
363 364 365 366 367 368 369 370 371 372 373 374 375 376
          // Build physical event
          m_NPDetectorManager->BuildPhysicalEvent();
  
          // Fill object in output ROOT file
          if (ccamPhys->EventMultiplicity > 0) {
  #ifdef __USE_CUTG__
            if (mcut->IsInside(ccamPhys->Half_Energy[0], ccamPhys->Calor_E[0])) {
              cc++;
              m_OutputTree->Fill();
            }
  #else
            m_OutputTree->Fill();
            cc++;
  #endif
377
            cout << cc << " " << c /*<< "(" << cr << ", " << cd << ") : " << tr << " " << td*/ << " |\t";
378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405
          }
  
          // check spectra
          m_NPDetectorManager->CheckSpectraServer();
          //cout << "d" << endl;
        }
        if (td < tr) {
          DD -> decodeEvent();
        } else {
          DR -> decodeBlobMFM();
        }
      } else if (cr < cd) {
        DR -> decodeBlobMFM();
      } else {
        DD -> decodeEvent();
      }
      if (DR -> getTime() < tr) {
        cr++;
      }
      tr = DR -> getTime();
      if (DD -> getTime() < td) {
        cd ++;
      }
      td = DD -> getTime();
    
      i++;
    
    } // End of main loop
406
//#endif
407
  } // End of mode if
408 409 410 411 412 413 414 415 416 417

  delete DR;
  delete DT;
  delete DD;
  delete [] rbuff;
  delete [] tbuff;

  // fill spectra
  m_NPDetectorManager -> WriteSpectra();

418 419 420 421 422
  if (argc > 1 and std::string(argv[1]) == "-r") {
    fout->cd();
    bidim->Write();
    deltaT->Write();
  } 
LAVIRON Adrien's avatar
LAVIRON Adrien committed
423
  fout->Close();
424
#endif
425 426 427 428 429 430 431 432 433 434
  // Essential
#if __cplusplus > 199711L && NPMULTITHREADING
  m_NPDetectorManager->StopThread();
#endif
  RootOutput::Destroy();

  return 0;
}