Ignore:
Timestamp:
11/04/11 15:12:50 (13 years ago)
Author:
neise
Message:
produces TH2F with an amplitude spektrum for all pixels
File:
1 edited

Legend:

Unmodified
Added
Removed
  • fact/tools/rootmacros/fpeak_cdf.C

    r12369 r12371  
    8282                                                // x = pixel id, y = DRS cell id
    8383
    84 TH1F *hBaseline[ NPIX ]; // histograms for baseline extraction
    85 TH1F *hMeanBsl, *hpltMeanBsl;
    86 TH1F *hRmsBsl, *hpltRmsBsl;
    8784TH2F * hAmplSpek_cfd;
    88 TH2F * hAmplSpek_discri;
    89 TH2F * hTemp_Array[1441];
    9085TObjArray hList;
    91 TObjArray hListBaseline, hListTemplates;
    9286
    9387void BookHistos( );
     
    10397//      * compute pulse height and pulse integral spektrum of the peaks
    10498int fpeak(
    105         char *datafilename              = "../../20111011_055.fits.gz",
    106         const char *drsfilename = "../../20111011_054.drs.fits.gz",
    107         const char *OutRootFileName = "./fpeak_cdf.Coutput.root",
     99        char *datafilename              = "data/20111016_013.fits.gz",
     100        const char *drsfilename = "../../20111016_011.drs.fits.gz",
     101        const char *OutRootFileName = "../analysis/fpeak_cdf.Coutput.root",
    108102        int firstevent                  = 0,
    109103        int nevents                     = -1,
     
    111105        int npixel                              = -1,
    112106        bool spikeDebug = false,
    113         int verbosityLevel = 1 // different verbosity levels can be implemented here
     107        int avg1                = 14,
     108        int avg2                = 8,
     109        int verbosityLevel = 1, // different verbosity levels can be implemented here
     110        bool ProduceGraphic = true
    114111        )
    115112{
     
    117114gROOT->SetStyle("Plain");
    118115
    119 // Create (pointer to) Canvases, which are used in every run,
    120 // also in 'non-debug' runs
    121         TCanvas * cSpektrum;
    122         TCanvas * cStartCell;
    123         TCanvas *cTemplate;
    124         cSpektrum = new TCanvas("cSpektrum","Amplitude Spektra of different discriminators",10,10,400,400);
    125         cSpektrum->Divide(1,2);
    126         cStartCell = new TCanvas("cStartCell ", "The Startcells of this run", 10,410,400,400);
    127         cTemplate = new TCanvas("cTemplate","Template of current Pixel",1,1,1600,1000);
    128 
    129         // Canvases only need if spike Debug, but I want to deklare
    130         // the pointers anyway ...
    131         TCanvas *cRawAndSpikeRemoval = NULL;
    132116        TCanvas *cFiltered = NULL;
    133 
    134 
    135117        if (spikeDebug){
    136                 cRawAndSpikeRemoval = new TCanvas("cRawAndSpikeRemoval","DRS Waveform",410,10,400,400);
    137                 cRawAndSpikeRemoval->Divide(1, 2);
    138118                cFiltered = new TCanvas("cFiltered","filtered DRS Waveforms",410,410,400,400);
    139                 cFiltered->Divide(1, 2);
     119                cFiltered->Divide(1, 3);
    140120        }
    141121
     
    188168                for ( int pix = firstpixel; pix < firstpixel+npixel; pix++ ){
    189169                        if (verbosityLevel > 0){
    190                                 if (pix % 20 ==0){
    191                                         cout << "Processing Event number: " << CurrentEventID << "\t"
    192                                                 << "Pixel number: "<< pix << endl;
     170                                if (pix == firstpixel){
     171                                        cout << "Processing Event: " << CurrentEventID << "/" << nevents << endl;
    193172                                }
    194173                        }
     
    205184
    206185                        // filter Vcorr with sliding average using FIR filter function
    207                         //factfir(b_slide , a_slide, k_slide, Vcorr, Vslide);
    208                         sliding_avg(Vcorr, Vslide, 8);
     186                        sliding_avg(Vcorr, Vslide, avg1);
    209187                        // filter Vslide with CFD using FIR filter function
    210188                        factfir(b_cfd , a_cfd, k_cfd, Vslide, Vcfd);
    211189                        // filter Vcfd with sliding average using FIR filter function
    212                         //factfir(bs2 , as2, ks2, Vcfd, Vcfd2);
    213                         sliding_avg(Vcfd, Vcfd2, 8);
     190                        sliding_avg(Vcfd, Vcfd2, avg2);
    214191
    215192
     
    220197                        vector<Region> * zXings = zerosearch( Vcfd2 , 1 , 8);
    221198                        // zXings means "zero cross ings"
    222                         //ShiftRegionBy(*zXings, -ks2/2);
    223                         EnlargeRegion(*zXings, 10, 10);
    224                         findAbsMaxInRegions(*zXings, Vslide);
    225                         //removeMaximaBelow( *zXings, 11.0, 0);
    226                         //removeMaximaAbove( *zXings, 14.0, 0);
     199                        EnlargeRegion( *zXings, 10, 10);
     200                        findAbsMaxInRegions( *zXings, Vslide);
     201                        removeMaximaBelow( *zXings, 3.0, 0);
     202                        removeRegionOnFallingEdge( *zXings, 100);
    227203
    228204                        // fill maxima in Histogram
     
    234210                                        }
    235211                                        hAmplSpek_cfd->Fill(pix, zXings->at(i).maxVal);
    236 
    237                                         for (int j=-100; j<150; j++){
    238                                                 if (zXings->at(i).maxPos + j >= 0 && zXings->at(i).maxPos + j <= 1023)
    239                                                         hTemp_Array[pix]->Fill(j+101, Vslide[zXings->at(i).maxPos + j]);
    240                                         }
    241212                                }
    242213                        }
     
    250221                                }
    251222
    252                                 cRawAndSpikeRemoval->cd( 1);
    253                                 gPad->SetGrid();
    254                                 debugHistos[Ameas_].Draw();
    255 
    256                                 cRawAndSpikeRemoval->cd( 2);
     223                                cFiltered->cd(1);
    257224                                gPad->SetGrid();
    258225                                debugHistos[Vcorr_].Draw();
    259226
    260                                 cFiltered->cd(1);
     227                                cFiltered->cd(2);
    261228                                gPad->SetGrid();
    262229                                debugHistos[Vslide_].Draw();
     
    278245                                }
    279246
    280                                 cFiltered->cd(2);
     247                                cFiltered->cd(3);
    281248                                gPad->SetGrid();
    282249                                debugHistos[Vcfd2_].Draw();
     
    285252                                zeroline->Draw();
    286253
    287                                 cRawAndSpikeRemoval->Update();
    288254                                cFiltered->Update();
    289 
    290                                 cTemplate->cd();
    291                                 hTemp_Array[pix]->Draw("COLZ");
    292                                 cTemplate->Modified();
    293                                 cTemplate->Update();
    294 
    295255                                //Process gui events asynchronously during input
    296256                                TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
     
    304264                                //TODO!!!!!!!!!
    305265                                // do some Garbage collection ...
    306                                 // all the Objects on the heap should be deleted here.
    307 
    308266                        }// end of if(spikeDebug)
    309267
     
    314272                } // end of loop over pixels
    315273
    316 if (ev % 10 == 0){
    317                 cSpektrum->cd(1);
    318                 hAmplSpek_cfd->Draw("COLZ");
    319                 cSpektrum->cd(2);
    320                 hAmplSpek_discri->Draw("COLZ");
    321                 cSpektrum->Modified();
    322                 cSpektrum->Update();
    323 
    324                 // updating seems not to work ..
    325                 // debug cout
    326                 cStartCell->cd();
    327                 hStartCell->Draw();
    328                 cStartCell->Modified();
    329                 cStartCell->Update();
    330 
    331                 cTemplate->cd();
    332                 hTemp_Array[firstpixel]->GetYaxis()->SetRangeUser(-10,40);
    333                 hTemp_Array[firstpixel]->Draw("COLZ");
    334                 cTemplate->Modified();
    335                 cTemplate->Update();
    336 }
    337 
    338 
    339274                if (breakout)
    340275                        break;
    341276        }       // end of loop over pixels
    342 
    343 
     277if (ProduceGraphic){
     278        TCanvas * cSpektrum;
     279        cSpektrum = new TCanvas("cSpektrum","Amplitude Spektra of different discriminators",10,10,400,400);
     280        cSpektrum->Divide(1,1);
     281  cSpektrum->cd(1);
     282  hAmplSpek_cfd->Draw("COLZ");
     283  cSpektrum->Modified();
     284  cSpektrum->Update();
     285}
    344286        SaveHistograms( OutRootFileName );
    345 
    346287        return( 0 );
    347288}
     
    412353void BookHistos( ){
    413354        // histograms for baseline extraction
    414         char hName[500];
    415         char hTitle[500];
    416355
    417356        TString name,title;
    418357        title = "all events all slices of pixel ";
    419358        name = "base";
    420         TH1F *h;
    421 
    422         for( int i = 0; i < NPIX; i++ ) {
    423                 title += i;
    424                 name += i;
    425 //              sprintf(&hTitle[0],"all events all slices of pixel %d", i);
    426 //              sprintf(&hName[0],"base%d", i);
    427 
    428                 h = new TH1F( name, title, 400, -99.5 ,100.5 );
    429 
    430                 h->GetXaxis()->SetTitle( "Sample value (mV)" );
    431                 h->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
    432                 hListBaseline.Add( h );
    433                 hBaseline[i] = h;
    434         }
    435 
    436     hMeanBsl = new TH1F("histo_mean","Value of maximal probability",400,-99.5,100.5);
    437     hMeanBsl->GetXaxis()->SetTitle( "max value (mV)" );
    438     hMeanBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
    439     hList.Add( hMeanBsl );
    440 
    441     hpltMeanBsl = new TH1F("hplt_mean","Value of maximal probability",1440,-0.5,1439.5);
    442     hpltMeanBsl->GetXaxis()->SetTitle( "pixel" );
    443     hpltMeanBsl->GetYaxis()->SetTitle( "max value in mV" );
    444     hList.Add( hpltMeanBsl );
    445 
    446     hRmsBsl = new TH1F("histo_rms","RMS in mV",2000,-99.5,100.5);
    447     hRmsBsl->GetXaxis()->SetTitle( "RMS (mV)" );
    448     hRmsBsl->GetYaxis()->SetTitle( "Entries / 0.5 mV" );
    449     hList.Add( hRmsBsl );
    450 
    451     hpltRmsBsl = new TH1F("hplt_rms","Value of maximal probability",1440,-0.5,1439.5);
    452     hpltRmsBsl->GetXaxis()->SetTitle( "pixel" );
    453     hpltRmsBsl->GetYaxis()->SetTitle( "RMS in mV" );
    454     hList.Add( hpltRmsBsl );
    455359
    456360        hAmplSpek_cfd = new TH2F("hAmplSpek_cfd","amplitude spektrum - CFD",1440,-0.5,1439.5, 256, -27.5, 100.5);
     
    459363        hList.Add( hAmplSpek_cfd );
    460364
    461         hAmplSpek_discri = new TH2F("hAmplSpek_discri","amplitude spektrum - std discriminator",1440,-0.5,1439.5, 256, -27.5, 100.5);
    462         hAmplSpek_discri->GetXaxis()->SetTitle( "pixel" );
    463         hAmplSpek_discri->GetXaxis()->SetTitle( "amplitude in mV" );
    464         hList.Add( hAmplSpek_discri );
    465365
    466366        debugHistos = new TH1F[ NumberOfDebugHistoTypes ];
     
    482382                debugHistos[ type ].GetYaxis()->SetTitle("Amplitude (a.u.)");
    483383        }
    484 
    485         TH2F *temp;
    486         for ( int type = 0; type < 1441; type++){
    487                 sprintf(&hTitle[0],"pulse template of pixel %d", type );
    488                 sprintf(&hName[0],"template_%d", type );
    489 
    490                 temp = new TH2F( hName, hTitle, 256, 0, 255, 256, -10.5, 117.5);
    491 
    492                 temp->GetXaxis()->SetTitle( "Time slice (%.1f ns/slice)" );
    493                 temp->GetYaxis()->SetTitle( "Amplitude (about mV)" );
    494                 hListTemplates.Add( temp );
    495                 hTemp_Array[ type ] = temp;
    496         }
    497 
    498 
    499384}
    500385
     
    504389
    505390        hList.Write(); // write the major histograms into the top level directory
    506         tf.mkdir("BaselineHisto");
    507         tf.cd("BaselineHisto"); // go to new subdirectory
    508         hListBaseline.Write(); // write histos into subdirectory
    509         tf.cd("..");
    510         tf.mkdir("PulseTempplates");
    511         tf.cd("PulseTemplates");
    512         hListTemplates.Write(); // write histos into subdirectory
    513391
    514392        tf.Close(); // close the file
    515 } // end of SaveHistograms(const char * loc_fname )
     393}
Note: See TracChangeset for help on using the changeset viewer.