Prototype2/hcalLab/hLabHelper.C

0001 hLabHelper:: hLabHelper(){
0002   //  whenever it is instanciated manually (not from detector helpers) it sets all dimensions to
0003   //  defaults one can change to get access to something in the data
0004   //  if request for Helper is from one of the project specific helpers then hLabHelper job is
0005   //  only to handle input/output related functions. The problem - all my detector specific
0006   //  helpers want to create detector with all presentation tools (and there are many) first
0007   //  and would like to have all that material available in the root file. Unfortunately it
0008   //  can't be created unique before the run number is known (it's a part of its name).
0009   //  I want all my presentation material to be in the TFile directory so they will go outside
0010   //   by default and so I create a temporary storage TFile with "host" dependent name for
0011   //  the period the program runs and use system tools to rename that file when program is
0012   //  ready for termination.
0013     runnumber      = 0;
0014     runKind        = "";     //   default
0015     eventseq       = 0;
0016     eventsread     = 0;
0017     fhcl           = NULL;
0018     thcl           = NULL;
0019     rootTempFileName = "";
0020     host = gSystem->HostName();
0021     host.Remove(host.First('.'));
0022
0023     if(ACTIVECHANNELS==0){
0024       CHTOTAL        = chinspected;
0025       ACTIVECHANNELS = CHTOTAL;
0026       detchannels    = CHTOTAL - triggerinsp;
0027       //  compute total number of HG &  LG channels in detector setup
0028       hgDetChannels = 0;
0029       lgDetChannels = 0;
0030       for(int ich = 0; ich<detchannels; ich++){
0031     if(feechinsp[ich]%2==HIGH) hgDetChannels ++;  else lgDetChannels ++;
0032       }
0033       cout<<"<hLabHelper:: hLabHelper >  Initializing as primary with "<<ACTIVECHANNELS<<" active channels"<<endl;
0034     } else {
0035       //  detectorHelper will taks care of all its data - it needs only file handling utilities from labHelper
0036       //  create and open buffer TFile so detector helper will be able to use it to place all its presentation
0037       //  data in its directory
0038       cout<<"<hLabHelper:: hLabHelper>  Initializing as secondary";
0039       rootTempFileName = (host.Contains("rcas"))?   RCFrootDir : HLABrootDir;
0040       TDatime T;
0041       TString fN = "hcalHelper_"; fN += T.GetDate(); fN += "_"; fN += T.GetTime(); fN += ".root";
0042       rootTempFileName += fN;
0043       cout<<"<hLabHelper:: hLabHelper> ROOT graphic objects will be written into TFile "<<rootTempFileName<<endl;
0044       fhcl = new TFile(rootTempFileName, "recreate");
0045       FileStat_t buf;
0046       if(gSystem->GetPathInfo(rootTempFileName, buf)) {
0047     cout<<"<runToRootFile> : File  "<<rootTempFileName<<"  still does not exist"<<endl;
0048       } else {cout<<"<hLabHelper:: hLabHelper> ROOT file "<<rootTempFileName<<"  created"<<endl;}
0049       fhcl->cd();
0050       return;
0051     }
0052     active     = new Int_t      [ACTIVECHANNELS];
0053     adc        = new Float_t  * [ACTIVECHANNELS];
0054     fitResults = new Double_t * [ACTIVECHANNELS];
0055     for (int ich=0; ich<ACTIVECHANNELS; ich++) {
0056       adc[ich]        = new Float_t  [NSAMPLES];
0057       fitResults[ich] = new Double_t [NPARAMETERS];
0058     }
0059     pedestal   = new Float_t  [ACTIVECHANNELS];   //     [ACTIVECHANNELS] */
0060     rawPeak    = new Float_t  [ACTIVECHANNELS];
0061     rawTime    = new Float_t  [ACTIVECHANNELS];
0062     fitPeak    = new Float_t  [ACTIVECHANNELS];
0063     fitTime    = new Float_t  [ACTIVECHANNELS];
0064     fitInt     = new Float_t  [ACTIVECHANNELS];
0065     fitChi2    = new Float_t  [ACTIVECHANNELS];
0066     fitPeakRMS = new Float_t  [ACTIVECHANNELS];
0067     fitTimeRMS = new Float_t  [ACTIVECHANNELS];
0068     rdata        = new TH2F("rdata", "Raw Peak Values for all active channels",ACTIVECHANNELS, 0, ACTIVECHANNELS, 1024, -2048., 2048.);
0069     shapes       = new TF1 *  [ACTIVECHANNELS];
0070     sigFit       = new TF1 *  [ACTIVECHANNELS];
0071     rvsc         = new TF1 *  [ACTIVECHANNELS];
0072     fmsc         = new TF1 *  [ACTIVECHANNELS];
0073     gpulse       = new TGraph * [ACTIVECHANNELS];
0074     ft           = new TH1 *  [ACTIVECHANNELS];
0075     fm           = new TH1 *  [ACTIVECHANNELS];
0076     fw           = new TH1 *  [ACTIVECHANNELS];
0077     fint         = new TH1 *  [ACTIVECHANNELS];
0078     fpd          = new TH1 *  [ACTIVECHANNELS];
0079     fchi2        = new TH1 *  [ACTIVECHANNELS];
0080     fitPar       = new TH1 ** [ACTIVECHANNELS];
0081     rpeak        = new TH1 *  [ACTIVECHANNELS];
0082     rtm          = new TH1 *  [ACTIVECHANNELS];
0083     for (int ich = 0; ich < ACTIVECHANNELS; ich++) {
0084       fitPar[ich] = new TH1 * [NPARAMETERS];
0085       shapes[ich] = NULL;
0086       gpulse[ich] = NULL;
0087     }
0088     for (int ig = 0; ig<2; ig++){
0089       evtadcsum[ig]  = new Double_t  [NSAMPLES];
0090       evtGraph[ig]   = new TGraph;
0091       TString sum    = "adcsum_"; sum += (ig==HIGH)? "HG"  :  "LG";
0092       TString sig    = "signal_"; sig += (ig==HIGH)? "HG"  :  "LG";
0093       evtShape[ig]   = new TF1(sum, &signalShape, 0., (Double_t)NSAMPLES, NPARAMETERS);
0094       evtSignal[ig]  = new TF1(sig, &signalShape, 0., (Double_t)NSAMPLES, NPARAMETERS);
0095     //    Double_t    times[ACTIVECHANNELS];
0096     //    Double_t    chi2 [ACTIVECHANNELS];
0097       evtfitpar[ig]  = new Double_t [NPARAMETERS];
0098     }
0099   }
0100 // **************************************************************************
0101
0102 bool hLabHelper::evDisplay(Int_t rn){
0103   if(rn<=0)    return false;
0104   runnumber    = rn;
0105   if(!runToPRDFFile(rn)) {
0106     cout<<"<setPRDFRun>  Failed to find the .prdf file for run number "<<rn<<endl;
0107     return false;
0108   }
0109   // Open up PRDFF
0110   int status;
0111   Eventiterator *it =  new fileEventiterator(prdfName, status);
0112   if (status)
0113   {
0114     cout << "Couldn't open input file " << prdfName << endl;
0115     delete(it);
0116     return false;
0117   }
0118   //  update channel map
0119   updateMap();
0120   int chperpanel = DISPLAYX*DISPLAYX;
0121   int panels     = (ACTIVECHANNELS%chperpanel==0?  ACTIVECHANNELS/chperpanel : ACTIVECHANNELS/chperpanel+1);
0122   int nx_c       = DISPLAYX;
0123   int ny_c(0);
0124   if(panels>1||ACTIVECHANNELS%chperpanel==0) {
0125     ny_c       = DISPLAYX;
0126   } else {
0127     ny_c       = ACTIVECHANNELS/DISPLAYX + (ACTIVECHANNELS%DISPLAYX==0? 0 : 1);
0128   }
0129
0130   //  cout<<"Display "<<panels<<"  "<<nx_c<<"/"<<ny_c<<endl;
0131
0132   gStyle->SetOptFit(111);
0133   gStyle->SetFitFormat("5.2g");
0134   //  creating Canvases
0135   TCanvas * evdisplay[panels];
0136   TString cnvsName  = "EvDisplay_";
0137   TString cnvsTitle = "sPhenix Channel display  for CHANNELS ";
0138   for (int iC = 0; iC < panels; iC++) {
0139     TString cN  = cnvsName;
0140     cN          += iC;
0141     TString cT  = cnvsTitle;
0142     cT          += (chperpanel*iC);
0143     cT  += " - ";
0144     cT  += (chperpanel*(iC+1)-1);
0145     evdisplay[iC]  = new TCanvas(cN, cT, 300*nx_c, 150*ny_c);
0146     evdisplay[iC]->SetEditable(kTRUE);
0147     evdisplay[iC]->Divide(nx_c,ny_c);
0148     //    if(iC==panels) {
0149     //    evdisplay[iC]->SetBit(TPad::kClearAfterCR);
0150     //    }
0151   }
0152   //   ====================================================================================
0153   //  now TGraph for every channel
0154   for (int ich = 0; ich<ACTIVECHANNELS;         ich++)
0155     {
0156       gpulse[ich] = new TGraph(NSAMPLES);
0157       TString name = "gch"; name += ich;
0158       gpulse[ich]->SetName(name);
0159       name = "Channel  "; name += ich; name += " / "; name += feechinsp[ich];  name += (feechinsp[ich]%2==HIGH? "  HG" : "  LG");
0160       gpulse[ich]->SetTitle(name);
0161       gpulse[ich]->SetMarkerStyle(20);
0162       gpulse[ich]->SetMarkerSize(0.4);
0163       //gpulse[ich]->SetMarkerColor(ich+1);
0164       //gpulse[ich]->SetLineColor(ich+1);
0165     }
0166   Event *evt;
0167
0168   // Loop over events
0169   while ( (evt = it->getNextEvent()) != 0 )
0170     {
0171       eventseq = evt->getEvtSequence();
0172
0173       if ( evt->getEvtType() != 1 ) {
0174     cout << "eventseq " << eventseq << " event type = " << evt->getEvtType() <<"  ignored "<< endl;
0175       } else {
0176     // Get sPHENIX Packet
0177     Packet_hbd_fpgashort *spacket = dynamic_cast<Packet_hbd_fpgashort*>( evt->getPacket(21101) );
0178     if ( spacket ) {
0179       spacket->setNumSamples( NSAMPLES );
0180       cout <<"PRDF FILE "<<prdfName<<"  EVENT "<<eventseq <<endl;
0181       for (int ich=0; ich<ACTIVECHANNELS; ich++) {
0182         for (int is=0; is<NSAMPLES; is++) {
0183           adc[ich][is] = spacket->iValue(active[ich],is);
0184           //          cout<<"<evDisplay  Event "<< eventseq<<"  ich/is/adc  "<<ich<<" - "<<is<<" - "<<adc[ich][is]<<endl;
0185         }
0186       }
0187       dofixes();
0188       getDetectorTiming();
0189       //      cout<<"<TileFit>  Peak  "<< evtsum.evtpeak<<"  Time "<<evtsum.evttime<<"  Chi2 "<<evtsum.tChi2 <<endl;
0190       displaySumEvent();
0191       //  we do generic display of square TCanvases DISPLAYX x DISPLAYX size
0192       //      for (int evD = 0; evD<panels; evD++){
0193       //        TString evDName  = cnvsName  + evD;
0194       //        TString evDTitle = cnvsTitle + evD;
0195       //        evdisplay->SetBit(TPad::kClearAfterCR);
0196       //        evdisplay->Divide(nx_c,ny_c);
0197
0198         for (int ich=0; ich<ACTIVECHANNELS; ich++) {
0199           for (int is=0; is<NSAMPLES; is++)  gpulse[ich]->SetPoint(is,(Double_t)is, adc[ich][is]);
0200           int iC   = ich/chperpanel;
0201           int iPad = (ich-chperpanel*iC);
0202           int iPx  = iPad%nx_c;
0203           int iPy  = ny_c - iPad/nx_c - 1;
0204           iPad     = iPy*nx_c + iPx + 1;
0205           evdisplay[iC] ->cd(iPad);
0206           //   cout<<iC<< "  "<<ich<<"   "<<iPx<<"  "<<iPy<<"  "<<iPad<<endl;
0207           gpulse[ich]->Draw("acp");
0208           // Double_t peakVal;
0209           // Int_t    peakPos;
0210           // fitShape(ich, peakVal, peakPos, (ich<detchannels? 1 : 0));
0211           // Double_t fMax  = shapes[ich]->GetMaximum(par0Min[1], par0Max[1]);
0212           // Double_t fPos  = shapes[ich]->GetMaximumX(par0Min[1], par0Max[1]);
0213           // Double_t fChi2 =  shapes[ich]->GetChisquare()/shapes[ich]->GetNDF();
0214           // cout<<"EVENT  "<<eventseq<<" SHAPE  "<<(char*)(shapes[ich]->GetName())<<"  peakVal  "<<peakVal<<"  peakPos  "<<peakPos<<"  fitMax  "<<fMax<<"  fitPos  "<<fPos<<"   fChi2  "<<fChi2<<endl;
0215           // shapes[ich]->Draw("same");
0216           //          evdisplay[iC]->Update();
0217         }
0218         for (int iC = 0; iC < panels; iC++) evdisplay[iC]->Update();
0219         char dummy;
0220         cin>>dummy;
0221         //    }
0222       delete spacket;
0223       //      for (int iC = 0; iC < panels; iC++) evdisplay[iC]->Clear();
0224     } else cout<<"<evDisplay> Packet hbd not found. EventSeq "<<eventseq<<" EventType "<< evt->getEvtType()<<endl;
0225       }
0226       delete evt;
0227     }
0228   cout<<"<evDisplay> EoF reached "<<endl;
0229
0230   return true;
0231 }
0232
0233 // **************************************************************************
0234 void hLabHelper::updateMap(){
0235   for (int ich = 0; ich<ACTIVECHANNELS;   ich++) active[ich] = feechinsp[ich]  ;
0236 }
0237
0238 // **************************************************************************
0239 void hLabHelper::setRunKind(TString kind){
0240   runKind = kind;
0241   cout<<"hLabHelper::setRunKind> will look for run data in "<<runKind<<endl;
0242 }
0243
0244
0245
0246 // **************************************************************************
0247 void hLabHelper::getFileLists(){
0248   //  cout<<"<getFileLists>  hostname  "<<host<<endl;
0249   //   ====================================================================================
0250   //  parsing .root file name from runnumber
0251   rootdirname  = (host.Contains("rcas"))?   RCFrootDir : HLABrootDir;
0252   prdfdirname  = (host.Contains("rcas"))?   RCFdataDir : HLABdataDir;
0253   //  in t1044 we have "junk", "cosmics", "led"  and "beam" kinds of events
0254   cout<<"<hLabHelper::getFileLists> "<<rootdirname<<endl;
0255   cout<<"<hLabHelper::getFileLists> "<<prdfdirname<<endl;
0256   if(!runKind.IsNull()) {
0257     prdfdirname += runKind; prdfdirname +="/";
0258     cout<<"hLabHelper::getFileLists> PRDF Directory name updated "<<prdfdirname<<endl;
0259   }
0260   TSystemDirectory rdir(rootdirname, rootdirname);
0261   TSystemDirectory pdir(prdfdirname, prdfdirname);
0262   TList * rootfiles = rdir.GetListOfFiles();
0263   TList * prdffiles = pdir.GetListOfFiles();
0264   Int_t   rn;
0265   TString rid;
0266   if(rootfiles){
0267     TSystemFile * file;
0268     TString       fname;
0269     TIter         next(rootfiles);
0270     while((file = (TSystemFile *) next())) {
0271       fname = file->GetName();
0272       if(!file->IsDirectory() && fname.EndsWith(".root")) {
0273     //  .root file found - get run number
0274     decoderun(fname, rn, rid);
0275     TString name = rootdirname;
0276     name += fname;
0277     rootfile rf(name, rn, rid);
0278     roots.push_back(rf);
0279     //  cout<<".root entry "<<name<<endl;
0280       } else {
0281     //  cout<<"<getFileList>  candidate root  "<<fname<<"  ignored "<<endl;
0282
0283       }
0284     }
0285   }
0286   if(prdffiles){
0287     TSystemFile * file;
0288     TString       fname;
0289     TIter         next(prdffiles);
0290     while((file = (TSystemFile *) next())) {
0291       fname = file->GetName();
0292       //      cout<<fname<<endl;
0293       if(!file->IsDirectory() && fname.EndsWith(".prdf")) {
0294     //  .prdf file found - get run number
0295     decoderun(fname, rn, rid);
0296     TString name = prdfdirname;
0297     name += fname;
0298     prdffile rf(name, rn, rid);
0299     prdfs.push_back(rf);
0300     //  cout<<".prdf entry "<<name<<endl;
0301       } else {
0302     //  cout<<"<getFileList>  candidate prdf  "<<fname<<"  ignored "<<endl;
0303
0304       }
0305     }
0306
0307   }
0308   cout<<"<getFileLists>  .root files found "<<roots.size()<<"   .prdf files found "<<prdfs.size()<<endl;
0309 }
0310
0311 // **************************************************************************
0312 void hLabHelper::decoderun(TString & fname, Int_t & rn, TString & rid){
0313   if(fname.Contains(".root") || fname.Contains(".prdf")) {
0314     //    cout<<"<decoderun> : HCAL file "<<fname<<"  confirmed"<<endl;
0315     rid          = fname;   // get run number from root file name
0316     rid.ReplaceAll(".root",""); rid.ReplaceAll(".prdf","");
0317     Ssiz_t lastslash = rid.Last('/');
0318     rid.Remove(0,lastslash+1);
0319     rid.ReplaceAll("HClab_","");   rid.ReplaceAll("rc-","");  rid.ReplaceAll("cosmics_",""); rid.ReplaceAll("beam_","");
0320
0321     TString id       = rid;
0322     id.Remove(rid.Last('-'));
0323     rn               = id.Atoi();
0324     //    cout<<"<decoderun> "<<fname<<"  runnumber "<<rn<<"  runId "<<rid<<endl;
0325   }
0326 }
0327
0328 // **************************************************************************
0329 TFile * hLabHelper::attachrootrun(Int_t rn){
0330   if(!roots.size()) {
0331     cout<<"<attachrootrun>  No .root files in the list: Collecting "<<endl;
0332     getFileLists();
0333   }
0334   //  find the rootfile rn in the local list of files kept by Helper
0335   std::list<rootfile>::iterator it = roots.begin();
0336
0337   while (it!=roots.end())  {
0338     //    it->print();
0339     if(it->runnumber==rn) {
0340       runnumber    = rn;
0341       runId        = it->runid;
0342       cout<<"<hLabHelper::attachrootrun> ROOT file with "<<runId<<" identifier found: "<<it->name<<endl;
0343       // this is dangerous
0344       return new TFile(it->name);
0345     } else it++;
0346   }
0347   cout<<"<hLabHelper::attachrootrun> Requested ROOT file for Run "<<rn<<" not found"<<endl;
0348   return NULL;
0349   // auto it = std::find_if( std::begin( roots ),
0350   //              std::end(   roots ),
0351   //                           [&]( const rootfile rf ){ return 1 == ( rf.runnumber==rn ); } );
0352   // if(roots.end() == it)  return NULL;
0353   //  return new TFile(it->name);
0354 }
0355
0356
0357
0358 //         const int pos = std::distance( myList.begin(), it ) + 1;
0359 //         std::cout << "item is found at position " << pos << std::endl;
0360
0361
0362 // **************************************************************************
0363 Bool_t hLabHelper::setPRDFRun(int run, Bool_t bookH){
0364   if(run<=0)   return false;
0365   runnumber    = run;
0366   if(!runToPRDFFile(run)) {
0367     cout<<"<setPRDFRun>  Failed to find the .prdf file for run number "<<run<<endl;
0368     return false;
0369   }
0370   initPRDFRun(bookH);
0371   return true;
0372 }
0373
0374
0375
0376 // **************************************************************************
0377 Bool_t hLabHelper::setROOTRun(int run){
0378   if(!runToRootFile(run)) {
0379     cout<<"<setROOTRun>  Failed to find the .root file for run number "<<run<<endl;
0380     return false;
0381   }
0382   return true;
0383 }
0384
0385
0386 // **************************************************************************
0387
0388 // The assumption hear is that the .root file is either already attached to this ROOT session or
0389 // must be looked for using its runnumber
0390 Int_t hLabHelper:: runToRootFile(Int_t run){
0391   if(run<=0) {
0392     cout<<"<runToRootFile> WARNING  :  Requested run number = "<<run<<endl;
0393     TSeqCollection * sq = gROOT->GetListOfFiles();
0394     if(!sq->IsEmpty()) {
0395       //  there are already attached files - let's see if we may find something over there
0396       //  now it is really .rootf into runId not what is in the name
0397       TString current = ((TObject *)(sq->First()))->GetName();
0398       //  check the extension first
0399       if(current.Contains(".root")) {
0400     cout<<"<hLabHelper::runToRootFile> :  attached .root file "<<current<<"  found"<<endl;
0401     rootfName   = current;
0402     runId       = current;   // get run number from root file name
0403     runId.ReplaceAll(".root","");
0404     Ssiz_t lastslash = runId.Last('/');
0405     runId.Remove(0,lastslash+1);
0406     runId.ReplaceAll("HClab_","");
0407     TString rn       = runId;
0408     rn.Remove(rn.Last('-'));
0409     runnumber        = rn.Atoi();
0410     return(runnumber);
0411       }
0412     } else    return 0;
0413   } else {
0414     //  TODO:  check if TFile for requested run is already attached
0415     //  let's try to find and attach file with the right runnumber
0416     //  find where I am
0417     cout<<"<runToRootFile> :  looking into .root directory for runnumber "<<run<<endl;
0418     //   ====================================================================================
0419     //  parsing .root file name from runnumber
0420     rootfName  = (host.Contains("rcas"))?   RCFrootDir : HLABrootDir;
0421     TString fname = rootfName;
0422     FileStat_t buf;
0423     runId = "00"; runId += run;  runId += "-0";
0424     fname += "HClab_";  fname += runId;  fname += ".root";
0425     if(gSystem->GetPathInfo(fname, buf)) {
0426       fname = rootfName;
0427       runId = "0"; runId += run;  runId += "-0";
0428       fname += "HClab_";  fname += runId;  fname += ".root";
0429       rootfName  = fname;
0430       if(gSystem->GetPathInfo(fname, buf)) {
0431     cout<<"<runToRootFile> : File  "<<fname<<"  does not exist"<<endl;
0432       } else  cout<<"<runToRootFile> : File  "<<fname<<"  already  exist and will be recreated"<<endl;
0433     }
0434     runnumber  = run;
0435   }
0436   if(!rootTempFileName.IsNull()) {
0437     cout<<"<hLabHelper:: runToRootFile> found temporaly root file "<<rootTempFileName<<endl; } else {
0438     fhcl = new TFile(rootfName,"recreate");
0439     cout<<"hLabHelper:: runToRootFile> "<<rootfName<<"  opened"<<endl;
0440   }
0441   return runnumber;
0442 }
0443
0444 // **************************************************************************
0445
0446 Int_t  hLabHelper::runToPRDFFile(Int_t run){
0447   //   ====================================================================================
0448   //  find where I am
0449   //   ====================================================================================
0450   //  parsing .prdf file name from runnumber
0451   prdfName  = (host.Contains("rcas"))?   RCFdataDir : HLABdataDir;
0452   if(!runKind.IsNull()) {
0453     prdfName += runKind; prdfName += "/";
0454     cout<<"hLabHelper::runToPRDFFile> PRDF Directory name updated "<<prdfName<<endl;
0455   }
0456   TString fname = prdfName;
0457   FileStat_t buf;
0458   // runId = "00"; runId += run;  runId += "-0";
0459   // fname += "rc-";  fname += runId;  fname += ".prdf";
0460   runId = "0000"; runId += run;  runId += "-0000";
0461   fname += runKind; fname += "_";
0462   fname += runId;  fname += ".prdf";
0463   if(gSystem->GetPathInfo(fname, buf)) {
0464     fname = prdfName;
0465     // runId = "00"; runId += run;  runId += "-0";
0466     // fname += "rc-";  fname += runId;  fname += ".prdf";
0467     runId = "0000"; runId += runnumber;  runId += "-0000";
0468     fname += runKind;   fname += "_";
0469     fname += runId;  fname += ".prdf";
0470     if(gSystem->GetPathInfo(fname, buf)) {
0471       cout<<"<runToPRDFFile> File  "<<fname<<"  with runId  "<<runId<<"  does not exist"<<endl;
0472       return 0;
0473     }
0474   }
0475   prdfName  = fname;
0476   //   ====================================================================================
0477   cout<<"<runToPRDFFile> Parsed prdf file name  "<<prdfName<<endl;
0478   //   ====================================================================================
0479
0480   cout << "<runToPRDFFile>  prdfName "<<prdfName << "  RunId  "<<runId<<"  Run Number  "<<runnumber<<endl;
0481
0482   //   ====================================================================================
0483   return  runnumber;
0484 }
0485
0486
0487 // **************************************************************************
0488 void hLabHelper::initPRDFRun(Bool_t bookH){
0489   runToRootFile(runnumber);
0490   //   ====================================================================================
0491   // initialize root file and TTree (the idea is to save all data and histos at the end of run)
0492   if(runnumber<=0) {
0493     cout<<"<hLabHelper::initPRDFRun> Cant instantiate run with runnumber = 0.  Exiting"<<endl;
0494     return;
0495   }
0496   //   ====================================================================================
0497   //  TApplication theApp("theApp",&argc,argv);
0498
0499   gROOT->SetStyle("Plain");
0500   /*
0501     gStyle->SetLabelSize(0.1,"xyz");
0502     gStyle->SetTitleSize(0.08,"xyz");
0503     gStyle->SetStatH(0.4);
0504   */
0505   if(!bookH) return;
0506   //  fit parameters (for templating)
0507   for (int ich = 0; ich<ACTIVECHANNELS; ich ++){
0508     TString fhn  = "ft_"; fhn += ich; TString fht = "Fitted pulse peak time CH="; fht += ich;
0509     ft[ich]    = new TH1F(fhn, fht,   100, 0., NSAMPLES);
0510     fhn        = "fm_"; fhn += ich; fht = "Fitted pulse amplitude CH="; fht += ich;
0511     fm[ich]    = new TH1F(fhn, fht,   4100, 0.,    4100.);
0512     fhn        = "fw_"; fhn += ich; fht = "Fitted pulse width CH="; fht += ich;
0513     fw[ich]    = new TH1F(fhn, fht,   100, 0.,    5.);
0514     fhn        = "fpd_"; fhn += ich; fht = "Fitted pedestal at t=PEAK CH="; fht += ich;
0515     fpd[ich]   = new TH1F(fhn, fht,   200, 2000., 2100.);
0516     fhn        = "fint_"; fhn += ich; fht = "Integral of fitted curve (pedestal suppressed) CH="; fht += ich;
0517     fint[ich]  = new TH1F(fhn, fht,   4000, 0., 4000.);
0518     fhn        = "fchi2_";fhn += ich; fht = "Chi2 of fit  CH="; fht += ich;
0519     fchi2[ich] = new TH1F(fhn, fht,   500, 0., 100.);
0520
0521     for(int ip=0; ip<NPARAMETERS; ip++){      TString phn = "p_"; phn += ip; phn += "_"; phn += ich;
0522       TString pht = "Fit Parameter "; pht += ip; pht += " Channel "; pht += ich;
0523       fitPar[ich][ip] = new TH1F(phn, pht, 100, par0Min[ip], par0Max[ip]);
0524     }
0525
0526     TString rvn = "rv_"; rvn += ich; TString rvt = "Raw peak value CH = "; rvt += ich;
0527     rpeak[ich] = new TH1F(rvn, rvt, 4100, 0., 4100.);
0528     TString rtn = "rt_"; rtn += ich; TString rtt = "Raw peak time CH = "; rtt += ich;
0529     rtm[ich]   = new TH1F(rtn, rtt, NSAMPLES, 0., NSAMPLES);
0530   }
0531   for (int ich=0; ich<ACTIVECHANNELS; ich++) {
0532     TString fn = "S_"; fn += ich;
0533     sigFit[ich] =  new TF1(fn, &signalShape, 0., (Double_t)NSAMPLES, 6);
0534     //    sigFit[ich]->AbsValue(true);
0535   }  //  This version of plotting is for High/Low gain preamplifier sending its outputs into two sequential readout channels
0536   //  int nx_c = 2;
0537   //  int ny_c = ACTIVECHANNELS/2+ACTIVECHANNELS%2;
0538   //  ev_display = new TCanvas("sphenix","sphenix display",400*nx_c,200*ny_c);
0539
0540   //  fitFunc   = new TCanvas("fitFunc","Fit Functions",  400*nx_c,200*ny_c);
0541   //  ev_display->SetEditable(kTRUE);
0542   //  ev_display->Divide(nx_c,ny_c);
0543   //  ev_display->SetBit(TPad::kClearAfterCR);
0544   //  fitFunc->SetEditable(kTRUE);
0545   //  fitFunc->Divide(nx_c,ny_c);
0546
0547   TString name, title;
0548   for (int ich=0; ich<ACTIVECHANNELS; ich++)
0549     {
0550       gpulse[ich] = new TGraph(NSAMPLES);
0551       name = "gch"; name += ich;
0552       gpulse[ich]->SetName(name);
0553       gpulse[ich]->SetMarkerStyle(20);
0554       gpulse[ich]->SetMarkerSize(0.4);
0555       //gpulse[ich]->SetMarkerColor(ich+1);
0556       //gpulse[ich]->SetLineColor(ich+1);
0557     }
0558
0559 }
0560
0561 //  -----------------------------------------------------------------------------------------------
0562 void hLabHelper::renameAndCloseRF(){
0563   if(rootTempFileName.IsNull()) return;
0564   fhcl -> Close();
0565   gSystem->Rename(rootTempFileName, rootfName);
0566   cout<<"<hLabHelper::copyAndCloseRootFile> Renamed  "<<endl;
0567   cout<<"           "<<rootTempFileName<<endl;
0568   cout<<"into  "<<endl;
0569   cout<<"           "<<rootfName<<endl;
0570 }
0571
0572
0573 //  -----------------------------------------------------------------------------------------------
0574 void hLabHelper::makeCanvasDirectory(){
0575   rootdirname  = (host.Contains("rcas"))?   RCFrootDir : HLABrootDir;
0576   cDirectory = rootdirname; cDirectory += "canvases/r";
0577   cDirectory += runnumber;  cDirectory += "/";
0578   FileStat_t buf;
0579   if(!gSystem->GetPathInfo(cDirectory, buf)) {
0580     cout<<"<hLabHelper::makeCanvasDirectory>   Directory "<<cDirectory<<" already exists "<<endl;
0581   } else {
0582     //    gSystem->MakeDirectory(cDirectory);
0583     TString command = "mkdir "; command += cDirectory;
0584     gSystem->Exec(command);
0585     cout<<"<hLabHelper::makeCanvasDirectory>  Directory "<<cDirectory<<" created "<<endl;
0586
0587   }
0588 }
0589 //  -----------------------------------------------------------------------------------------------
0590
0591 void hLabHelper::fitShape(int ich, Double_t & rVal, Int_t & rTime, Int_t mode){
0592   if(!shapes[ich]){
0593     TString fitName  = active[ich]%2? "HG_" : "LG_";
0594     fitName += ich/2;
0595     fitName += "(";
0596     fitName += ich;
0597     fitName += ")";
0598     shapes[ich] = new TF1(fitName, &signalShape, 0., (Double_t)NSAMPLES, 6);
0599   }
0600   //  GAIN range
0601   Int_t ig = feechinsp[ich]%2==HIGH? HIGH  :  LOW;
0602   //  zero approximation to fit parameters
0603   //  use pulse data to find location of maximum
0604   //  cout<<"<fitShape> Event "<<hclt->eventseq<<" channel "<<ich<<endl;
0605   rTime = 0;
0606   rVal  = 0.;
0607   Int_t iss(RISETIME-2), ise(NSAMPLES-FALLTIME);
0608   Int_t rMaxVal  = adc[ich][iss]; Int_t maxsmpl = iss;
0609   Int_t rMinVal  = adc[ich][iss]; Int_t minsmpl = iss;
0610
0611   for (int is = iss; is<=ise; is++) {
0612     if(adc[ich][is]>rMaxVal) {
0613       rMaxVal  = adc[ich][is];
0614       maxsmpl  = is;
0615     }
0616     if(adc[ich][is]<rMinVal) {
0617       rMinVal  = adc[ich][is];
0618       minsmpl  = is;
0619     }
0620     //    cout<<ich<<"  "<<is<<"  "<<adc[ich][is]<<"  "<<rVal<<"  "<<rTime<<endl;
0621   }
0622   rMinVal  -= PEDESTAL;   rMaxVal -= PEDESTAL;
0623   if(abs(rMinVal) > abs(rMaxVal)) {
0624     rVal = rMinVal;  rTime = minsmpl;
0625   } else {
0626     rVal = rMaxVal;  rTime = maxsmpl;
0627   }
0628
0629   //  cout<<"  rVal-ped  "<<rVal<<endl;
0630
0631   par0[0] = rVal/3.;
0632   if(mode == 0 ) {
0633     par0[1] = max(0.,(Double_t)(rTime-RISETIME));
0634     par0[2] = 4.;
0635     par0[3] = 1.6;
0636     par0[4] = PEDESTAL;
0637   //   we do not do pedestal subtrastion at this time
0638     par0[5] = 0.;      //   slope of the pedestals is initially set to "0"
0639     shapes[ich]->SetParameters(par0);
0640     for(int parn=0; parn<6; parn++) shapes[ich]->SetParLimits(parn, par0Min[parn], par0Max[parn]);
0641   } else {
0642     //  channels which belong to detector
0643     par0[1] = evtfitpar[ig][1];
0644     par0[2] = evtfitpar[ig][2];
0645     par0[3] = evtfitpar[ig][3];
0646     par0[4] = PEDESTAL;
0647     par0[5] = 0.;
0648     shapes[ich]->SetParameters(par0);
0649     shapes[ich]->SetParLimits(2, evtfitpar[ig][2], evtfitpar[ig][2]);
0650     shapes[ich]->SetParLimits(3, evtfitpar[ig][3], evtfitpar[ig][3]);
0651   }
0652   //  fitFunc->cd(ich+1);  //  shapes[ich]->Draw();
0653   //  gPad->Modified();
0654   gpulse[ich]->Fit(shapes[ich], "RQ", "NQ", 0., (Double_t)NSAMPLES);
0655   shapes[ich]->GetParameters(fitResults[ich]);
0656   //  ev_display->cd(ich+1);
0657   //  shapes[hHelper->ich]->Draw("same");
0658
0659 }
0660
0661
0662 // **************************************************************************
0663 //   The next few functions are to process data from HCal Lab recorded via HBD electronics
0664
0665 void hLabHelper::collect(Bool_t fitshape){
0666   //  Fit parameters
0667
0668   //  Double_t uSum(0.), cSum(0.), luSum(0.), lcSum(0.), ruSum(0.), rcSum(0.);
0669   //  Double_t YU(0.), YC(0.),  XU(0.), XC(0.);
0670   //  hits above threshold in this event depending on the threshold (binned in Chi units - 8 counts)
0671   Double_t  peakVal;
0672   Int_t     peakPos;
0673   for (int ich=0; ich<ACTIVECHANNELS; ich++) {
0674     for (int is=0; is<NSAMPLES; is++){
0675       gpulse[ich]->SetPoint(is,(Double_t)is, adc[ich][is]);
0676     }
0677     if(fitshape) fitShape(ich, peakVal, peakPos, (ich<detchannels? 1 : 0));
0678     // fitShape(ich, peakVal, peakPos, 1);
0679   }
0680
0681   for (int ich=0; ich<ACTIVECHANNELS; ich++){
0682     int      rTime = 0;
0683     Double_t rVal  = 0.;
0684     //    Int_t iss(RISETIME-2), ise(NSAMPLES-FALLTIME);
0685     Int_t iss(0), ise(NSAMPLES-1);
0686     Int_t rMaxVal  = adc[ich][iss]; Int_t maxsmpl = iss;
0687     Int_t rMinVal  = adc[ich][iss]; Int_t minsmpl = iss;
0688
0689     for (int is = iss; is<=ise; is++) {
0690       if(adc[ich][is]>rMaxVal) {
0691     rMaxVal  = adc[ich][is];
0692     maxsmpl  = is;
0693       }
0694       if(adc[ich][is]<rMinVal) {
0695     rMinVal  = adc[ich][is];
0696     minsmpl  = is;
0697       }
0698     }
0699     rMinVal  -= PEDESTAL;   rMaxVal -= PEDESTAL;
0700     if(abs(rMinVal) > abs(rMaxVal)) {
0701       rVal = rMinVal;  rTime = minsmpl;
0702     } else {
0703       rVal = rMaxVal;  rTime = maxsmpl;
0704     }
0705     //   cout<<"<collect>  "<<ich<<"  "<<rTime<<"  "<<adc[ich][(int)rTime]<<endl;
0706     rawTime[ich] = rTime;
0707     rawPeak[ich] = rVal;
0708     rawPeak[ich] -= (fitshape? (fitResults[ich][4]+fitResults[ich][5]*rTime) : PEDESTAL);   //
0709     rpeak[ich]  -> Fill(rawPeak[ich]);    rtm[ich]  -> Fill(rawTime[ich]);
0710     rdata       -> Fill(ich, rawPeak[ich]);
0711     //  Below only in case if fitshape is required
0712     if(fitshape){
0713       sigFit[ich]->SetParameters(fitResults[ich]);
0714       sigFit[ich]->SetParameter(4, 0.);
0715       sigFit[ich]->SetParameter(5, 0.);
0716
0717       //      Double_t pVal  = sigFit[ich]->GetMaximum(par0Min[1], par0Max[1]);
0718       //      Double_t pTime = sigFit[ich]->GetMaximumX(par0Min[1], par0Max[1]);
0719       Double_t sgp  =  sigFit[ich]->GetMaximum(par0Min[1], par0Max[1]);
0720       Double_t tmp  =  sigFit[ich]->GetMaximumX(par0Min[1], par0Max[1]);
0721       Double_t sgn  =  sigFit[ich]->GetMinimum(par0Min[1], par0Max[1]);
0722       Double_t tmn  =  sigFit[ich]->GetMinimumX(par0Min[1], par0Max[1]);
0723       Double_t pVal    = rVal;
0724       Double_t pTime   = rTime;
0725       if(tmp!=0&&abs(sgp)>abs(sgn)) {
0726     pVal   = sgp;
0727     pTime  = tmp;
0728     //  cout<<"<getDetectorTimingT>: Positive peak selected"<<endl;
0729       } else {
0730     pVal   = sgn;
0731     pTime  = tmn;
0732     //  cout<<"<getDetectorTimingT>: Negative peak selected"<<endl;
0733       }
0734       Double_t xmin    = max(0.,pTime-RISETIME);
0735       Double_t xmax    = min(Double_t(NSAMPLES), pTime+FALLTIME);
0736       Double_t pInt    = (pVal!=0?  sigFit[ich]->Integral(xmin, xmax)  : 0);
0737       // if(pInt==0) {
0738       //   cout<<"<collect> Integral zero Event "<<hclt->eventseq<<"  Ch "<<ich<<" pVal/pTime/xmin/xmax "<<pVal<<"  "<<pTime<<"  "<<xmin<<"  "<<xmax<<endl;
0739       // }
0740       //    pVal          -= fped;
0741       fitPeak[ich]   = pVal;
0742       fitTime[ich]   = pTime;
0743       fitInt[ich]    = pInt;
0744       fitChi2[ich]   = shapes[ich]->GetChisquare()/shapes[ich]->GetNDF();
0745     //  trigger contribution
0746       ft[ich]      -> Fill(pTime);
0747       fm[ich]      -> Fill(pVal);
0748     //    Double_t mean = sigFit[ich]->Mean(xmin, xmax);
0749       Double_t mom2 = (pVal!=0?  sigFit[ich]->CentralMoment(2, xmin, xmax)  : 0.);
0750       // if(mom2==0) {
0751       //   cout<<"<collect> CentralMoment zero Event "<<hclt->eventseq<<"  Ch "<<ich<<" pVal/pTime/pInt/xmin/xmax "<<pVal<<"  "<<pTime<<"  "<<pInt<<"  "<<xmin<<"  "<<xmax<<endl;
0752       // }
0753       Double_t rms  = sqrt(mom2);
0754       fw[ich]      -> Fill(rms);
0755       fpd[ich]     -> Fill(fitResults[ich][4]+pTime*fitResults[ich][5]);
0756       fint[ich]    -> Fill(pInt);
0757       fchi2[ich]   -> Fill(fitChi2[ich]);
0758
0759       for (int ip=0; ip<NPARAMETERS; ip++)  fitPar[ich][ip]->Fill(fitResults[ich][ip]);
0760       //    Double_t fVal = GetMaximum( par0Min[2], par0Max[2]);
0761     }
0762   }
0763
0764   if(!reject()){
0765     //  few histogramms for clean sample
0766   };
0767 }
0768
0769 // **************************************************************************
0770
0771 void hLabHelper::getDetectorTiming(){
0772   //  build eventsum (adc sum without pedestal subtruction)
0773   //  do fitting etc (HG/LG separately)
0774
0775   for(int is = 0; is < NSAMPLES; is++){
0776     evtadcsum[0][is] = 0.;
0777     evtadcsum[1][is] = 0.;
0778     for (int ich = 0; ich < detchannels; ich++) {
0779       Int_t ig = feechinsp[ich]%2==HIGH? HIGH  :  LOW;
0780       evtadcsum[ig][is] += adc[ich][is];
0781     }
0782   }
0783   for (int is = 0; is < NSAMPLES; is++) {
0784     evtGraph[HIGH]->SetPoint(is,(Double_t)is,evtadcsum[HIGH][is]);
0785     evtGraph[LOW] ->SetPoint(is,(Double_t)is,evtadcsum[LOW][is]);
0786   }
0787   Double_t rVal, rTime;
0788   for(int ig = 0; ig<2; ig++) {
0789     Int_t iss(RISETIME-2), ise(NSAMPLES-FALLTIME);
0790     Int_t rMaxVal  = evtadcsum[ig][iss]; Int_t maxsmpl = iss;
0791     Int_t rMinVal  = evtadcsum[ig][iss]; Int_t minsmpl = iss;
0792
0793     for (int is = iss; is<=ise; is++) {
0794       if(evtadcsum[ig][is]>rMaxVal) {
0795     rMaxVal  = evtadcsum[ig][is];
0796     maxsmpl  = is;
0797       }
0798       if(evtadcsum[ig][is]<rMinVal) {
0799     rMinVal  = evtadcsum[ig][is];
0800     minsmpl  = is;
0801       }
0802       //    cout<<ich<<"  "<<is<<"  "<<adc[ich][is]<<"  "<<rVal<<"  "<<rTime<<endl;
0803     }
0804     cout<<"getDetectorTiming:  Event "<<eventseq<<"  Gain "<<ig<<"  Min "<<minsmpl<<" / "<<rMinVal<<"  Max "<<maxsmpl<<" / "<<rMaxVal<<endl;
0805     rMinVal  -= (ig==HIGH?  PEDESTAL*hgDetChannels   :  PEDESTAL*lgDetChannels);;
0806     rMaxVal  -= (ig==HIGH?  PEDESTAL*hgDetChannels   :  PEDESTAL*lgDetChannels);;
0807     if(abs(rMinVal) > abs(rMaxVal)) {
0808       rVal = rMinVal;  rTime = minsmpl;
0809     } else {
0810       rVal = rMaxVal;  rTime = maxsmpl;
0811     }
0812     cout<<"getDetectorTimingT:  Signal at "<<rTime<<"  Amplitude "<<rVal<<endl;
0813
0814     //    if(rVal<0.) rVal = 0.;
0815     Double_t  par0[6];
0816     par0[0] = rVal;      //3.;
0817     par0[1] = max(0.,(Double_t)(rTime-RISETIME));
0818     par0[2] = 4.;
0819     par0[3] = 1.6;
0820     par0[4] = (ig==HIGH?  PEDESTAL*hgDetChannels   :  PEDESTAL*lgDetChannels);  //   we do not do pedestal subtrastion at this time
0821     par0[5] = 0;      //   slope of the pedestals is initially set to "0"
0822     evtShape[ig]->SetParameters(par0);
0823     for(int parn=0; parn<6; parn++) evtShape[ig]->SetParLimits(parn, par0Min[parn], par0Max[parn]);
0824     evtShape[ig]->SetParLimits(4, par0[4]-PEDESTAL, par0[4]+PEDESTAL);
0825     //  fitFunc->cd(chan+1);  //  shapes[chan]->Draw();
0826     //  gPad->Modified();
0827     evtGraph[ig]->Fit(evtShape[ig], "RQWM", "NQ", 0., (Double_t)NSAMPLES);
0828     evtShape[ig]->GetParameters(evtfitpar[ig]);
0829     cout<<"<getDetectorTiming>: Gain "<<ig<<"  Parameters (0) "<<evtfitpar[ig][0]<<" (1) "<<evtfitpar[ig][1]<<" (2) "<<evtfitpar[ig][2]<<" (3) "<<evtfitpar[ig][3]<<" (4) "<<evtfitpar[ig][4]<<" (5) "<<evtfitpar[ig][5]<<endl;
0830     evtpedestal[ig]  = (evtfitpar[ig][4]+evtfitpar[ig][5]*rTime);   //
0831     evtSignal[ig]->SetParameters(evtfitpar[ig]);
0832     evtSignal[ig]->SetParameter(4, 0.);
0833     evtSignal[ig]->SetParameter(5, 0.);
0834     Double_t sgp  = evtSignal[ig]->GetMaximum(par0Min[1], par0Max[1]);
0835     Double_t tmp  = evtSignal[ig]->GetMaximumX(par0Min[1], par0Max[1]);
0836     Double_t sgn  = evtSignal[ig]->GetMinimum(par0Min[1], par0Max[1]);
0837     Double_t tmn  = evtSignal[ig]->GetMinimumX(par0Min[1], par0Max[1]);
0838     evtpeak[ig]  = rVal;
0839     evttime[ig]  = rTime;
0840     cout<<"<getDetectorTimingT>: Gain "<<ig<<"  Peak0 "<<  evtpeak[ig]<<"  Time0 "<< evttime[ig]<<endl;
0841     if(tmp!=0&&abs(sgp)>abs(sgn)) {
0842       evtpeak[ig]  = sgp;
0843       evttime[ig]  = tmp;
0844       cout<<"<getDetectorTimingT>: Positive peak selected"<<endl;
0845     } else {
0846       evtpeak[ig]  = sgn;
0847       evttime[ig]  = tmn;
0848       cout<<"<getDetectorTimingT>: Negative peak selected"<<endl;
0849     }
0850     tChi2[ig]    = evtShape[ig]->GetChisquare()/evtShape[ig]->GetNDF();
0851     cout<<"<getDetectorTimingT>: Gain "<<ig<<"  Peak "<<  evtpeak[ig]<<"  Time "<< evttime[ig]<<"  Chi2 "<<tChi2[ig]<<endl;
0852   }
0853
0854 }
0855 // **************************************************************************
0856
0857 void hLabHelper::displaySumEvent(){
0858
0859   Int_t ng           = (lgDetChannels!=0&&hgDetChannels!=0)? 2  :  1;
0860   TCanvas * eventSum = (TCanvas *)(gROOT->FindObject("eventSum"));
0861   if(!eventSum) eventSum  = new TCanvas("eventSum","sPhenix eventSum display", 300, 300, 800, 500);
0862   else eventSum->Clear();
0863   eventSum -> Divide(1,ng);
0864   for(int ig = 0; ig < ng; ig++){
0865     eventSum->cd(ig+1);
0866     //    for (int is = 0; is < NSAMPLES; is++)  evtGraph[ig]->SetPoint(is,(Double_t)is,evtadcsum[ig][is]);
0867     //  cout<<"<tileDisplay>   call to fitTileSignal"<<endl;
0868     //  eventSum->Divide(1,2);
0869     //  eventSum->cd(1);
0870     //  eventSum->SetEditable(kTRUE);
0871     evtGraph[ig]->SetMarkerStyle(20);
0872     evtGraph[ig]->SetMarkerSize(0.4);
0873     eventSum->cd(ig+1);
0874     evtGraph[ig]->Draw("acp");
0875     //  eventSum->cd(2);
0876     evtShape[ig]->Draw("same");
0877     // evtShape[ig]->Draw();
0878     // evtGraph[ig]->Draw("same");
0879     eventSum->Update();
0880   }
0881   eventSum->Update();
0882   //  eventSum->Clear();
0883 }
0884 // **************************************************************************
0885 //  FIXES to data
0886 void hLabHelper::dofixes(){
0887   // if(runnumber>=1125 && runnumber<1152){
0888   //   //  invert pulse in channel 134 (number 9 in todays parlance
0889   //   for (int is = 0; is < NSAMPLES; is++) {adc[8][is] -= 2*PEDESTAL; adc[8][is] = abs(adc[8][is]);}
0890   // }
0891   ;
0892 }
0893
0894 // **************************************************************************
0895
0896 Bool_t hLabHelper::reject(){
0897   return false;
0898 }
0899 // **************************************************************************