TPC/groot/MyFavoriteMartin.C

0001 #include <iostream>
0002 #include <cmath>
0003 #include <pmonitor.h>
0004 #include "MyFavoriteMartin.h"
0005
0006 #include "groot.h"
0007 #include "ATrace.h"
0008 #include "AZigzag.h"
0009 #include "FindBlobs.h"
0010 #include "FillBlobHist.h"
0011 #include "FillHoughHist.h"
0012
0013 #include "TH1.h"
0014 #include "TH2.h"
0015 #include "TH3.h"
0016 #include "TNtuple.h"
0017 #include "TFile.h"
0018
0019 #include "prdfoStream.h"
0020 #include "Quiver.h"
0021
0022 #include "params.h"
0023
0024 //#include <stdlib.h>
0025 #include <cstdlib>
0026 #include <ctime>
0027 #include <map>
0028
0029 prdfoStream *os = 0;
0030 int init_done = 0;
0031
0032 //KD
0033 int first_event_nr = 0;
0034 //KD
0035
0036 using namespace std;
0037
0038 //  These are the waveforms...
0039 //  They are the most basic monitor of the SRS system.
0040 TH1F *h1a;
0041 TH1F *h1b;
0042 TH1F *h1c;
0043 TH1F *h1d;
0044
0045 TH1F *h2a;
0046 TH1F *h2b;
0047 TH1F *h2c;
0048 TH1F *h2d;
0049
0050 TH1F *h3a;
0051 TH1F *h3b;
0052 TH1F *h3c;
0053 TH1F *h3d;
0054
0055 TH1F *h4a;
0056 TH1F *h4b;
0057 TH1F *h4c;
0058 TH1F *h4d;
0059
0060 //  These "scope" histograms are the most basic monitor of the DRS4 system
0061 TH1F *scope0, *scope1, *scope2, *scope3;
0062
0063 //  These hnl_D histograms are used to phase the PLL...
0064 TH1F *h1a_D;
0065 TH1F *h1b_D;
0066 TH1F *h1c_D;
0067 TH1F *h1d_D;
0068 TH1F *h2a_D;
0069 TH1F *h2b_D;
0070 TH1F *h2c_D;
0071 TH1F *h2d_D;
0072 TH1F *h3a_D;
0073 TH1F *h3b_D;
0074 TH1F *h3c_D;
0075 TH1F *h3d_D;
0076 TH1F *h4a_D;
0077 TH1F *h4b_D;
0078 TH1F *h4c_D;
0079 TH1F *h4d_D;
0080
0081 int pinit()
0082 {
0083   // Here we will book all the histograms, but NONE of the
0084   // Ntuples.  The reason is that we only want to make the NTuples
0085   // when we've made an output file.
0086
0087   // Generate the First groot and initialize the Zigzags.
0088   groot *Tree = groot::instance();
0089
0090   scope0 = new TH1F ( "scope0","waveform Ch 0", 1024, -0.5, 1023.5);
0091   scope1 = new TH1F ( "scope1","waveform Ch 1", 1024, -0.5, 1023.5);
0092   scope2 = new TH1F ( "scope2","waveform Ch 2", 1024, -0.5, 1023.5);
0093   scope3 = new TH1F ( "scope3","waveform Ch 3", 1024, -0.5, 1023.5);
0094
0095   scope0->SetLineColor(kYellow+2);
0096   scope1->SetLineColor(kBlue);
0097   scope2->SetLineColor(kViolet);
0098   scope3->SetLineColor(kGreen);
0099
0100   h1a = new TH1F( "h1a", "Waveform", 4500, -0.5, 4499.5);
0101   h1b = new TH1F( "h1b", "Waveform", 4500, -0.5, 4499.5);
0102   h1c = new TH1F( "h1c", "Waveform", 4500, -0.5, 4499.5);
0103   h1d = new TH1F( "h1d", "Waveform", 4500, -0.5, 4499.5);
0104
0105   h2a = new TH1F( "h2a", "Waveform", 4500, -0.5, 4499.5);
0106   h2b = new TH1F( "h2b", "Waveform", 4500, -0.5, 4499.5);
0107   h2c = new TH1F( "h2c", "Waveform", 4500, -0.5, 4499.5);
0108   h2d = new TH1F( "h2d", "Waveform", 4500, -0.5, 4499.5);
0109
0110   h3a = new TH1F( "h3a", "Waveform", 4500, -0.5, 4499.5);
0111   h3b = new TH1F( "h3b", "Waveform", 4500, -0.5, 4499.5);
0112   h3c = new TH1F( "h3c", "Waveform", 4500, -0.5, 4499.5);
0113   h3d = new TH1F( "h3d", "Waveform", 4500, -0.5, 4499.5);
0114
0115   h4a = new TH1F( "h4a", "Waveform", 4500, -0.5, 4499.5);
0116   h4b = new TH1F( "h4b", "Waveform", 4500, -0.5, 4499.5);
0117   h4c = new TH1F( "h4c", "Waveform", 4500, -0.5, 4499.5);
0118   h4d = new TH1F( "h4d", "Waveform", 4500, -0.5, 4499.5);
0119
0120   h1a_D = new TH1F( "h1a_D", "Waveform", 4001, -0.5, 4000.5);
0121   h1b_D = new TH1F( "h1b_D", "Waveform", 4001, -0.5, 4000.5);
0122   h1c_D = new TH1F( "h1c_D", "Waveform", 4001, -0.5, 4000.5);
0123   h1d_D = new TH1F( "h1d_D", "Waveform", 4001, -0.5, 4000.5);
0124   h2a_D = new TH1F( "h2a_D", "Waveform", 4001, -0.5, 4000.5);
0125   h2b_D = new TH1F( "h2b_D", "Waveform", 4001, -0.5, 4000.5);
0126   h2c_D = new TH1F( "h2c_D", "Waveform", 4001, -0.5, 4000.5);
0127   h2d_D = new TH1F( "h2d_D", "Waveform", 4001, -0.5, 4000.5);
0128   h3a_D = new TH1F( "h3a_D", "Waveform", 4001, -0.5, 4000.5);
0129   h3b_D = new TH1F( "h3b_D", "Waveform", 4001, -0.5, 4000.5);
0130   h3c_D = new TH1F( "h3c_D", "Waveform", 4001, -0.5, 4000.5);
0131   h3d_D = new TH1F( "h3d_D", "Waveform", 4001, -0.5, 4000.5);
0132   h4a_D = new TH1F( "h4a_D", "Waveform", 4001, -0.5, 4000.5);
0133   h4b_D = new TH1F( "h4b_D", "Waveform", 4001, -0.5, 4000.5);
0134   h4c_D = new TH1F( "h4c_D", "Waveform", 4001, -0.5, 4000.5);
0135   h4d_D = new TH1F( "h4d_D", "Waveform", 4001, -0.5, 4000.5);
0136
0137   return 0;
0138 }
0139
0140
0141 int process_event (Event * e)
0142 {
0143   // Set the run number if event is from new run...
0144   int thisRun = e->getRunNumber();
0145   if( Quiver::RunNumber != thisRun ) Quiver::RunNumber = thisRun;
0146
0147   // Line below is used to put NON-DATA events into output stream
0148   if ( os && e->getEvtType() >7 ) os->addEvent(e);
0149
0150   //  Don't run physics analysis code on NON-DATA events!
0151   if (e->getEvtType() > 7) return 0;
0152
0153   //Use nr-th event as starting point of process_event
0154   // and skip every event before that
0155   static int skipentries=0;
0156   skipentries++;
0157   static bool SetSkipEntry = true;
0158   if(SetSkipEntry)
0159     {
0160       cout << " Skipping " << first_event_nr << " events."
0161        << endl;
0162       SetSkipEntry = false;
0163     }
0164   if(skipentries%100==0 && skipentries<first_event_nr)
0165     cout << " Skipped " << skipentries << " events." << endl;
0166   if(e->getEvtSequence() < first_event_nr && e->getEvtType()==1)
0167     return 0;
0168
0169   // Here is a data integrity check for the SRS system
0170   // Use this for offline analysis at a mature level only.
0171   static bool PriorPassed = true;
0172   Packet *p = e->getPacket(1010);
0173   if (p)
0174     {
0175       int nhb = p->iValue(0,"NHYBRIDS");
0176
0177       //  DIagnostic...
0178       // cout << "Number of Hybrids:  " << nhb << endl;
0179       // cout << "Number of Samples:  " << endl;
0180       // for (int i=0; i<nhb; i++)
0181       //    {
0182       //      int nsmp = p->iValue(i,"NSAMPLES");
0183       //      cout << "  " << nsmp;
0184       //    }
0185       // cout << endl;
0186
0187       //  Check that the data has the right number of hybrids.
0188       if (nhb != Nhybrid)
0189     {
0190       if (PriorPassed)
0191         {
0192           cout << "Wrong number of hybrids: "
0193            << nhb << " found but expected "
0194            << Nhybrid << ", abort event."
0195            << endl;
0196         }
0197       else
0198         {
0199           cout << ".";
0200         }
0201       delete p;
0202       PriorPassed = false;
0203       return 0;
0204     }
0205
0206       //  Check that all hybrids have the same data length as the first one.
0207       for (int i=1; i<nhb-1; i++)
0208     {
0209       if (p->iValue(0,"NSAMPLES") != p->iValue(i,"NSAMPLES"))
0210         {
0211           if (PriorPassed)
0212         {
0213           cout << "Data length mismatch: "
0214                << p->iValue(i,"NSAMPLES")
0215                << " out of " << p->iValue(0,"NSAMPLES")
0216                << " on hybrid " << i
0217                << " abort event."
0218                << endl;
0219         }
0220           else
0221         {
0222           cout << ".";
0223         }
0224           delete p;
0225           PriorPassed = false;
0226           return 0;
0227         }
0228     }
0229       if (!PriorPassed) cout << endl;
0230       PriorPassed = true;
0231     }
0232
0233   // Fancy event counter written to the screen...
0234   static int entries=0;
0235   entries++;
0236   if (entries%10==0)
0237     {
0238       cout << "Processed " << entries;
0239       time_t now = time(0);
0240       char* dt = ctime(&now);
0241       cout << " at " << dt << endl;
0242     }
0243
0244   groot *Tree = groot::instance();
0245   Tree->ClearTheDetector();
0246   Tree->event = e;
0247
0248   scope0->Reset();
0249   scope1->Reset();
0250   scope2->Reset();
0251   scope3->Reset();
0252
0253   h1a->Reset();
0254   h1b->Reset();
0255   h1c->Reset();
0256   h1d->Reset();
0257
0258   h2a->Reset();
0259   h2b->Reset();
0260   h2c->Reset();
0261   h2d->Reset();
0262
0263   h3a->Reset();
0264   h3b->Reset();
0265   h3c->Reset();
0266   h3d->Reset();
0267
0268   h4a->Reset();
0269   h4b->Reset();
0270   h4c->Reset();
0271   h4d->Reset();
0272
0273   // Reading the PSI 4-channel "oscilloscope"...
0274   Packet *p1 = e->getPacket(1020);
0275   if (p1 && Quiver::AnalyzeScope)
0276     {
0277       int samples = p1->iValue(0, "SAMPLES");
0278       for ( int i = 0; i < samples; i++)
0279     {
0280
0281       scope0->Fill ( i, -p1->rValue(i,0) );
0282       scope1->Fill ( i, -p1->rValue(i,1) );
0283       scope2->Fill ( i, -p1->rValue(i,2) );
0284       scope3->Fill ( i, -p1->rValue(i,3) );
0285
0286       Tree->theTraces[0]->voltage.push_back(-p1->rValue(i,0));
0287       Tree->theTraces[1]->voltage.push_back(-p1->rValue(i,1));
0288       Tree->theTraces[2]->voltage.push_back(-p1->rValue(i,2));
0289       Tree->theTraces[3]->voltage.push_back(-p1->rValue(i,3));
0290
0291     }
0292     }
0293   if (p1) delete p1;
0294
0295   //  Now monitor and fill SRS to groot.
0296   if (p && Quiver::AnalyzeSRS)
0297     {
0298
0299       int nhb = p->iValue(0,"NHYBRIDS");
0300       for (int i=0; i<nhb; i++)
0301     {
0302       if (p->iValue(i,"NSAMPLES")<24)
0303         {
0304           cout << "ERROR in Sample Length, "<< p->iValue(i,"NSAMPLES") <<" abort event" << endl;
0305           return 0;
0306         }
0307     }
0308
0309       static bool FirstEvent = true;
0310       if (FirstEvent)
0311     {
0312       cout << "Reports for " << nhb << " Hybrids:" << endl;
0313       for (int i=0; i<nhb; i++)
0314         {
0315           cout << i << ": " << p->iValue(i,"NSAMPLES") << endl;
0316         }
0317       FirstEvent = false;
0318     }
0319
0320
0321       //  Check the hybrid count for sensibility...
0322       if (nhb>=1)
0323     {
0324       int idigi;
0325       for (idigi = 0; idigi< 4500; idigi++)
0326         {
0327           h1a->Fill( idigi, p->iValue(idigi, 0, "RAWSAMPLES"));
0328           h1b->Fill( idigi, p->iValue(idigi, 1, "RAWSAMPLES"));
0329           h1c->Fill( idigi, p->iValue(idigi, 2, "RAWSAMPLES"));
0330           h1d->Fill( idigi, p->iValue(idigi, 3, "RAWSAMPLES"));
0331
0332           h2a->Fill( idigi, p->iValue(idigi, 4, "RAWSAMPLES"));
0333           h2b->Fill( idigi, p->iValue(idigi, 5, "RAWSAMPLES"));
0334           h2c->Fill( idigi, p->iValue(idigi, 6, "RAWSAMPLES"));
0335           h2d->Fill( idigi, p->iValue(idigi, 7, "RAWSAMPLES"));
0336
0337           h3a->Fill( idigi, p->iValue(idigi, 8, "RAWSAMPLES"));
0338           h3b->Fill( idigi, p->iValue(idigi, 9, "RAWSAMPLES"));
0339           h3c->Fill( idigi, p->iValue(idigi,10, "RAWSAMPLES"));
0340           h3d->Fill( idigi, p->iValue(idigi,11, "RAWSAMPLES"));
0341
0342           h4a->Fill( idigi, p->iValue(idigi, 12, "RAWSAMPLES"));
0343           h4b->Fill( idigi, p->iValue(idigi, 13, "RAWSAMPLES"));
0344           h4c->Fill( idigi, p->iValue(idigi, 14, "RAWSAMPLES"));
0345           h4d->Fill( idigi, p->iValue(idigi, 15, "RAWSAMPLES"));
0346
0347           if (idigi>500 && idigi<3500)
0348         {
0349           h1a_D->Fill( p->iValue(idigi, 0, "RAWSAMPLES"));
0350           h1b_D->Fill( p->iValue(idigi, 1, "RAWSAMPLES"));
0351           h1c_D->Fill( p->iValue(idigi, 2, "RAWSAMPLES"));
0352           h1d_D->Fill( p->iValue(idigi, 3, "RAWSAMPLES"));
0353           h2a_D->Fill( p->iValue(idigi, 4, "RAWSAMPLES"));
0354           h2b_D->Fill( p->iValue(idigi, 5, "RAWSAMPLES"));
0355           h2c_D->Fill( p->iValue(idigi, 6, "RAWSAMPLES"));
0356           h2d_D->Fill( p->iValue(idigi, 7, "RAWSAMPLES"));
0357           h3a_D->Fill( p->iValue(idigi, 8, "RAWSAMPLES"));
0358           h3b_D->Fill( p->iValue(idigi, 9, "RAWSAMPLES"));
0359           h3c_D->Fill( p->iValue(idigi,10, "RAWSAMPLES"));
0360           h3d_D->Fill( p->iValue(idigi,11, "RAWSAMPLES"));
0361           h4a_D->Fill( p->iValue(idigi,12, "RAWSAMPLES"));
0362           h4b_D->Fill( p->iValue(idigi,13, "RAWSAMPLES"));
0363           h4c_D->Fill( p->iValue(idigi,14, "RAWSAMPLES"));
0364           h4d_D->Fill( p->iValue(idigi,15, "RAWSAMPLES"));
0365         }
0366         }
0367
0368       // Read the tracker raw data into the raw array...
0369       int rindex=0;
0370       for (int JINX =0; JINX<nhb; JINX++)
0371         {
0372           for ( int j = 0; j< CHhybrid; j++)
0373         {
0374           for (int i=0; i<min(p->iValue(JINX,"NSAMPLES"),24); i++) // limit to 24 samples...
0375             {
0376               //  For now we'll use a trivial routine to analyze waveform.
0377               //  We'll just sum across the samples...
0378               AZigzag::Raw[rindex].push_back(p->iValue(j,i,JINX));
0379             }
0380           rindex++;
0381         }
0382         }
0383       AZigzag::DetermineCommonMode();
0384       AZigzag::ApplyCalibration();
0385       for (int i=0; i<Tree->theZigzags.size(); i++)
0386         {
0387           Tree->theZigzags[i]->DetermineQ();
0388         }
0389       FindBlobs();
0390       FillBlobHist();
0391       //FillHoughHist();
0392     }
0393       else
0394     {
0395       cout << "ERROR " << nhb << " Hybrids Found:" << endl;
0396       for (int i=0; i<nhb; i++)
0397         {
0398           cout << i << ": " << p->iValue(i,"NSAMPLES") << endl;
0399         }
0400     }
0401     }
0402
0403
0404   // Here select whether to write the event to the output stream.
0405   static int nSaved=0;
0406   if ( os && Tree->theTracks.size() >= 1 )
0407     {
0408       nSaved++;
0409       os->addEvent(e);
0410       if (nSaved%100==0)
0411     {
0412       cout << "Saved " << nSaved;
0413       cout << " percentage "
0414            << 100.0*(double)nSaved/(double)entries;
0415       cout << endl;
0416     }
0417     }
0418   if (p) delete p;
0419
0420   return 0;
0421 }
0422
0423
0424
0425 int prdfopen ( const char *filename)
0426 {
0427   if (os)
0428     {
0429       cout << "file is already open " << endl;
0430       return -1;
0431     }
0432
0433   os = new prdfoStream (filename);
0434   if ( os->is_defunct() )
0435     {
0436       delete os;
0437       os = 0;
0438       return -1;
0439     }
0440
0441   return 0;
0442 }
0443
0444 int prdfclose ()
0445 {
0446   if ( !os)
0447     {
0448       cout << "no file open " << endl;
0449       return -1;
0450     }
0451   delete os;
0452   os = 0;
0453   return 0;
0454 }
0455
0456
0457
0458 //Function that allows to start prun() from the nr-th event
0459 int getFirstEventNr() { return first_event_nr; }
0460
0461 void setFirstEventNr(const int nr) { first_event_nr = nr; }
0462
0463 // Binary number to GrayCode (courtesy of Martin)
0464 int b2G(int num) { return (num >> 1) ^ num; }