ROCsoft: beamtime/cern-oct12/go4/RICH/TRICHProc.cxx (r4864/r4735) Source File

Go to the documentation of this file.
00001 /* Generated by Together */
00002 
00003 #include "TRICHProc.h"
00004 #include "TCBMBeamtimeEvent.h"
00005 #include "TGo4Log.h"
00006 
00007 #include "TSystem.h"
00008 #include "TGo4UserException.h"
00009 #include "TGo4MbsEvent.h"
00010 #include "roc/Message.h"
00011 #include "roc/Board.h"
00012 
00013 #include "TGo4Analysis.h"
00014 
00015 #include "CbmRichRingFitterCOP.h"
00016 #include "CbmRichRingFitterEllipseTau.h"
00017 #include "CbmRichRingFinderHoughImpl.h"
00018 
00019 #include "TEllipse.h"
00020 #include "TLatex.h"
00021 #include <sstream>
00022 #include <fstream>
00023 
00024 TRICHProc::TRICHProc() :
00025    TCBMBeamtimeProc(),
00026    fRocInputEvent(0),
00027    fTrbInputEvent(0),
00028    fOutputEvent(0),
00029    fHodo1(0),
00030    fHodo2(0),
00031    fTriglogEvent(0),
00032    fBeamEvent(0),
00033    fEventNr(0),
00034    fNrSED(0),
00035    fEventNr_1stloop(0),
00036    fEventNr_2ndloop(0),
00037    nof2hitevents(0)
00038 {
00039 }
00040 
00041 
00042 TRICHProc::TRICHProc(const char* name) :
00043    TCBMBeamtimeProc(name),
00044    fRocInputEvent(0),
00045    fTrbInputEvent(0),
00046    fOutputEvent(0),
00047    fHodo1(0),
00048    fHodo2(0),
00049    fTriglogEvent(0),
00050    fBeamEvent(0),
00051    fEventNr(0),
00052    fNrSED(0),
00053    fEventNr_1stloop(0),
00054    fEventNr_2ndloop(0),
00055    nof2hitevents(0)
00056 {
00057    TGo4Log::Info("TRICHProc: Create instance %s", name);
00058       
00059    fPar = (TRICHParam*) MakeParameter("RICHPar", "TRICHParam", "set_RICHPar.C");
00060       
00061    if (!ReadMapmtGeometry(fPar->mapFileName))
00062    {
00063       TGo4Log::Info("First step of reading MAPMT geometry failed");
00064       char newMapFileName[256];
00065       sprintf (newMapFileName, "RICH/%s", fPar->mapFileName);
00066       if (!ReadMapmtGeometry(newMapFileName)) {
00067          TGo4Log::Info("Second step of reading MAPMT geometry failed\n");
00068       }
00069    }
00070    
00071    hEntries2dIntegral = MakeTH2('D', "MAPMT/Entries2dIntegral", "Integral of ADC values per MAPMT channels", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00072    hEntries2dMean = MakeTH2('D', "MAPMT/Entries2dMean", "Mean ADC values per MAPMT channels", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00073    
00074    hEntriesperPMT = MakeTH1('I', "MAPMT/EntriesperPMT", "Number of hits per MAPMT", 24, 0., 24., "mapmt", "hits");
00075    hEntriesperPMT2D = MakeTH2('I', "MAPMT/EntriesperPMT2D", "Number of hits per MAPMT", 4, 0., 32., 4, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00076    
00077    hADCall = MakeTH1('I', "MAPMT/ADCall", "Distribution of corrected ADC values", 100, 0., 4096.);
00078    hHittime = MakeTH1('I', "MAPMT/Hittime", "Distribution of difference between AUX2 and MAPMT hit", 1024, -1500, 4096.);
00079 
00080    for (Int_t iPMT=1; iPMT<=16; iPMT++) {
00081       hHitMultperPMT[iPMT]=MakeTH1('I', Form("MAPMT/HitMultperPMT/HitMultperPMT%i",iPMT), Form("Hit Multiplicity on PMT %i",iPMT),50,-0.5,49.5, "x (mm)", "y (mm)");
00082    }
00083       
00084    hNofHitsInHodo = MakeTH1('F', "MAPMT/hNofHitsInHodo", "Number of hits in hodoscope", 4, -0.5, 3.5, "Number of hits in hodoscope");
00085    hHodoADCvsNofHitInEv = MakeTH2('F', "MAPMT/hHodoADCvsNofHitInEv", "ADC value in hodoscope vs. number of hits in event", 6000, 0., 6000., 60, 0., 60., "ADC in hodoscope", "Number of hits in event");
00086    hNotAssignedHitsXY = MakeTH2('F', "MAPMT/hNotAssignedHitsXY", "Not assigned hits", 60, 0., 300., 48, 0., 240., "x-pos (mm)", "y-pos (mm)");
00087 
00088    // book additional 2d histograms for position plots with condition on fiber hod
00089    hMAPMT_hcor2d_s1 = MakeTH2('I', "MAPMT/mapmt_detector_h1", "Number of hits per MAPMT channels, hcut1", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00090    hMAPMT_hcor2d_s2 = MakeTH2('I', "MAPMT/mapmt_detector_h2", "Number of hits per MAPMT channels, hcut2", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00091    hMAPMT_hcor2d_s3 = MakeTH2('I', "MAPMT/mapmt_detector_h3", "Number of hits per MAPMT channels, hcut3", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00092    hMAPMT_hcor2d_s4 = MakeTH2('I', "MAPMT/mapmt_detector_h4", "Number of hits per MAPMT channels, hcut4", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00093    hMAPMT_hcor2d_s5 = MakeTH2('I', "MAPMT/mapmt_detector_h5", "Number of hits per MAPMT channels, hcut5", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00094    
00095    hNofFoundRings= MakeTH1('F', "MAPMT/hNofFoundRings", "Number of found rings", 3, -0.5, 2.5, "Number of found rings");
00096 
00097    hTimespread=MakeTH1('I',"MAPMT/Timespread","Timespread of hits",50,-25,25,"time (ns)","hits");
00098    
00099    //[0] is electron, [1] is pion, [2] is muon, [3] is all
00100    string p[] = {"e", "pi", "mu", "all"};
00101    for (Int_t i = 0; i < 4; i++){
00102       //book circle fitting histograms
00103       hCircleFitRadius[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/CircleFitRadius_")+p[i]).c_str(), "Radius of fitted rings", 501, -0.5, 100.5, "ring radius [mm]");
00104       hCircleFitRadius_cor[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/CircleFitRadius_cor_")+p[i]).c_str(), "Radius of fitted rings (corrected)", 501, -0.5, 100.5, "ring radius [mm]"); // temp and press corrected
00105       hCircleFitCenter[i] = MakeTH2('F', (string("MAPMT/fitting/")+p[i]+string("/CircleFitCenter_")+p[i]).c_str(), "Center of fitted rings", 208, 0., 208., 208, 0., 208., "x-pos [mm]", "y-pos [mm]");
00106       hCircleFitChi2[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/CircleFitChi2_")+p[i]).c_str(), "Chi2 of fitted rings", 200, 0., 40., "chi2");
00107       hCircleFitDR[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/CircleFitDR_")+p[i]).c_str(), "dR of fitted rings", 200, -20., 20., "dR [mm]");
00108       
00109       //book ellipse fitting histograms
00110       hEllipseFitAaxis[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitAaxis_")+p[i]).c_str(), "A axis of fitted ellipse", 501, -0.5, 100.5, "A axis [mm]");
00111       hEllipseFitBaxis[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitBaxis_")+p[i]).c_str(), "B axis of fitted ellipse", 501, -0.5, 100.5, "B axis [mm]");
00112       hEllipseFitPhi[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitPhi_")+p[i]).c_str(), "Phi of fitted ellipses", 400, -2., 2., "Phi [rad]");
00113       hEllipseFitBoverA[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitBoverA_")+p[i]).c_str(), "B/A of fitted ellipses", 100, 0., 1., "B/A");
00114       hEllipseFitCenter[i] = MakeTH2('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitCenter_")+p[i]).c_str(), "Center of fitted ellipses", 208, 0., 208., 208, 0., 208., "x-pos [mm]", "y-pos [mm]");
00115       hEllipseFitChi2[i] = MakeTH1('F', (string("MAPMT/fitting/")+p[i]+string("/EllipseFitChi2_")+p[i]).c_str(), "Chi2 of fitted ellipses", 200, 0., 40., "chi2");
00116 
00117       //number of hits per event for all events, for events with good circle fit, for event with good ellipse fit
00118       hNofHitsEvAll[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvAll_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
00119       hNofHitsEvRingFinder[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvRingFinder_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
00120       hNofHitsEvCircleFit[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvCircleFit_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
00121       hNofHitsEvEllipseFit[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvEllipseFit_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
00122       hNofHitsEvRingFinderEff[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvRingFinderEff_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
00123       hNofHitsEvCircleFitEff[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvCircleFitEff_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
00124       hNofHitsEvEllipseFitEff[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsEvEllipseFitEff_")+p[i]).c_str(), "Nof hits in event", 60, 0, 60, "Nof hits in event");
00125 
00126       hNofHitsRingFinder[i] = MakeTH1('F', (string("MAPMT/efficiency/")+p[i]+string("/NofHitsRingFinder_")+p[i]).c_str(), "Nof hits in ring", 60, 0, 60, "Nof hits in ring");
00127 
00128       hEntries2d[i]=MakeTH2('I', (string("MAPMT/Entries2d_")+p[i]).c_str(), "Hits per Pixel", 40, 0., 40., 32, 0., 32., "x-pos [pixel]", "y-pos [pixel]");
00129    }
00130    
00131    hNofHitsRingFinder_cor = MakeTH1('F', "MAPMT/efficiency/all/NofHitsRingFindercor_all", "Nof hits in ring (corrected)", 60, 0, 60, "Nof hits in ring (corrected)");
00132    //hNofHitsRingFinder_cor2 = MakeTH1('F', "MAPMT/efficiency/all/NofHitsRingFindercor2_all", "Nof hits in ring (corrected)", 60, 0, 60, "Nof hits in ring (corrected)");
00133    //hNofHitsRingFinder_cor3 = MakeTH1('F', "MAPMT/efficiency/all/NofHitsRingFindercor3_all", "Nof hits in ring (corrected)", 60, 0, 60, "Nof hits in ring (corrected)");
00134    
00135    
00137    hMAP_from_TRB = MakeTH2 ('I', "MAPMT/TRB_Entries2d", "", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00138    hEntries_TRBROC = MakeTH2 ('I', "MAPMT/TRBROC_Entries2d", "", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00140 
00142    hLedPulserEntriesPerChannel = MakeTH1 ('I', "MAPMT/LedPulser/LPEntriesPerChannel", "", 1280, 0., 1280.);
00143    hLedPulserEntries2d = MakeTH2 ('I', "MAPMT/LedPulser/LPEntries2d", "", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00144    hLedPulserHitVsMapmt = MakeTH2 ('I', "MAPMT/LedPulser/LPHitVsMapmt", "", 24, 0, 24, 64, 0., 64., "mapmt", "hits");
00145    for (UInt_t i = 0; i<24; i++) {
00146       hLedPulserADCvsPixel[i] = MakeTH2 ('I', Form("MAPMT/LedPulser/LPADCvsPixel/LP_mapmt%i", i+1), "", 65, -0.5, 64.5, 100, 0., 1500., "pixel", "ADC");
00147    }
00149 
00150    // book histograms for various single photon spectra
00151    for(int i=1; i<=24; i++) {
00152         hSinglePhotonSpectra[i] = MakeTH1('I', Form("MAPMT/LedPulser/SP_allEvents/SP_ADC_mapmt%i",i), "single photon spectra of all pixels (all-hit-events)", 160, -100., 1500.);
00153         hLEDPulserADC_1hit[i] = MakeTH1('I', Form("MAPMT/LedPulser/SP_1-Hit-Events/ADC_1Hit_mapmt%i",i), "ADC spectra of all pixels for 1-hit-events", 160, -100., 1500.); // ADC histogram for 1-hit-events
00154         hLEDPulserADC_nofph[i] = MakeTH1('I', Form("MAPMT/LedPulser/SPnorm_nOfPh/SP_nofph_mapmt%i",i), "Single-photon spectra in dep. of nof photons", 200, 0., 4.); // single-photon spectra in dep. of nof photons
00155         hLEDPulserADC_nofph_mct[i] = MakeTH1('I', Form("MAPMT/LedPulser/SPnorm_nOfPh/MinusCT/SP_nofph_mct_mapmt%i",i), "Single-photon spectra in dep. of nof photons, minus ct", 200, 0., 4.); // minus crosstalk
00156         hLEDPulserADC_nofphPP[i] = MakeTH2('I', Form("MAPMT/LedPulser/SPnorm_nOfPh/SP_nofph_PP_mapmt%i",i), "Single-photon spectra in dep. of nof photons", 65., -0.5, 64.5, 200, 0., 4.); // single-photon spectra in dep. of nof photons
00157         hLEDPulserADC_nofph_1hitpMA[i] = MakeTH1('I', Form("MAPMT/LedPulser/SPnorm_nOfPh/1hitpMA/SP_nofph_1hitpMA_mapmt%i",i), "Single-photon spectra in dep. of nof photons, 1 hit per mapmt", 100, 0., 4.); // 1 hit per mapmt
00158         hLEDPulserADC_1hit_perPixel[i] = MakeTH2('I', Form("MAPMT/LedPulser/SP_1-Hit-perPixel/ADC_1Hit_perPixel_mapmt%i",i), "ADC spectra of all (separated) pixel for 1-hit-events", 64, 0., 64., 160, -100., 1500., "Pixel", "ADC value");
00159         hLEDPulserADC_ctHits[i] = MakeTH1('I', Form("MAPMT/LedPulser/SP_Crosstalk-Hits/ADC_ctHit_mapmt%i",i), "ADC spectra of all pixels only for crosstalk-hits", 160, -100., 1500.); // ADC histogram for crosstalk-hits only 
00160         //hLEDPulserADC_ctHits[i] = MakeTH1('I', Form("MAPMT/LedPulser/SP_Crosstalk-Hits/ADC_ctHit_mapmt%i",i), "ADC spectra of all pixels only for crosstalk-hits", 200, 0., 4.); // ADC histogram for crosstalk-hits only 
00161         hLEDPulserADCcor_1hit[i] = MakeTH1('I', Form("MAPMT/LedPulser/SPcor_1-Hit-Events/ADCcor_1Hit_mapmt%i",i), "ADC spectra of all pixels for 1-hit-events (cor)", 160, -100., 1500.); // ADC histogram for 1-hit-events, gainPerPixel-corrected
00162    }
00163    
00164    // book test histograms
00165    for(int i=1; i<=24; i++) {
00166         SP_ADC_all[i] = MakeTH1('I', Form("MAPMT/LedPulser/test/mapmt%i_all",i), "ADC spectra: all hits", 160, -100., 1500.); 
00167         SP_ADC_ct[i] = MakeTH1('I', Form("MAPMT/LedPulser/test/mapmt%i_ct",i), "ADC spectra: crosstalk hits", 160, -100., 1500.); 
00168         SP_ADC_1hit[i] = MakeTH1('I', Form("MAPMT/LedPulser/test/mapmt%i_1hit",i), "ADC spectra: 1-hit", 160, -100., 1500.); 
00169         SP_ADC_ctmain[i] = MakeTH1('I', Form("MAPMT/LedPulser/test/mapmt%i_ctmain",i), "ADC spectra: hits with crosstalk in neighbours", 160, -100., 1500.); 
00170         SP_ADC_noct[i] = MakeTH1('I', Form("MAPMT/LedPulser/test/mapmt%i_noct",i), "ADC spectra: all hits except crosstalk hits", 160, -100., 1500.); 
00171         
00172         SP_norm_all[i] = MakeTH1('I', Form("MAPMT/LedPulser/test_norm/mapmt%i_norm_all",i), "ADC spectra: all hits", 250, -1., 4.); 
00173         SP_norm_ct[i] = MakeTH1('I', Form("MAPMT/LedPulser/test_norm/mapmt%i_norm_ct",i), "ADC spectra: crosstalk hits", 250, -1., 4.); 
00174         SP_norm_1hit[i] = MakeTH1('I', Form("MAPMT/LedPulser/test_norm/mapmt%i_norm_1hit",i), "ADC spectra: 1-hit", 250, -1., 4.); 
00175         SP_norm_ctmain[i] = MakeTH1('I', Form("MAPMT/LedPulser/test_norm/mapmt%i_norm_ctmain",i), "ADC spectra: hits with crosstalk in neighbours", 250, -1., 4.); 
00176         SP_norm_noct[i] = MakeTH1('I', Form("MAPMT/LedPulser/test_norm/mapmt%i_norm_noct",i), "ADC spectra: all hits except crosstalk hits", 250, -1., 4.); 
00177         crosstalk[i] = {0};
00178    }
00179    
00180    
00181    
00182    // book histograms for crosstalk analysis
00183    if(fPar->AnalyseCrosstalk) {
00184         hNEdge = MakeTH1('I', "MAPMT/Crosstalk/hNEdge", "n_edge", 11, -0.5, 10.5, "hits per event");
00185         hNCorner = MakeTH1('I', "MAPMT/Crosstalk/hNCorner", "n_corner", 11, -0.5, 10.5, "hits per event");
00186         hNCenter = MakeTH1('I', "MAPMT/Crosstalk/hNCenter", "n_center", 11, -0.5, 10.5, "hits per event");
00187         for (Int_t iPMT=1; iPMT<=24; iPMT++) {
00188                 hHitMultperPMT_CrTalk[iPMT]=MakeTH1('I', Form("MAPMT/Crosstalk/Crosschecks/HitMultperPMT_CrTalk%i",iPMT), Form("Hit Multiplicity on PMT%i CrTalk",iPMT),50,-0.5,49.5, "hit multiplicity", "#");
00189                 hEntries2dperPMT[iPMT] = MakeTH2('I', Form("MAPMT/Crosstalk/Crosschecks/Entries2dperPMT%i",iPMT), "Number of hits per MAPMT channel", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00190                 hHitDist[iPMT] = MakeTH1('D', Form("MAPMT/Crosstalk/HitDist_zones/HitDistPMT%i",iPMT), Form("Distance between hits via zones, PMT%i",iPMT), 15, -0.5, 14.5, "distance between 2 hits (pixel)");
00191                 hHitDist_py[iPMT] = MakeTH1('D', Form("MAPMT/Crosstalk/HitDist_pytha/HitDistPMT%i_py",iPMT), Form("Distance between hits via pythagoras, PMT%i",iPMT), 15, -0.5, 14.5, "distance between 2 hits (pixel)");
00192                 //hCrTalkDist[iPMT] = MakeTH2('I', Form("MAPMT/Crosstalk/hCrTalkDistPMT%i",iPMT), "Distribution of crosstalk hits", 32, 0., 32., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00193                 crosstalk_results[iPMT] = MakeTH1('D', Form("MAPMT/Crosstalk/Results/Ct_results_PMT%i",iPMT), Form("Results for crosstalk, PMT%i",iPMT), 24, -0.5, 23.5, "PMT");
00194                 
00195                 if(fPar->SimulateCrosstalk) {
00196                         hHitDist_sim[iPMT] = MakeTH1('D', Form("MAPMT/Crosstalk/Simulation/HitDistPMT%i_sim",iPMT), Form("Simulated distance between hits via zones, PMT%i",iPMT), 15, -0.5, 14.5, "distance between 2 hits (pixel)");
00197                         hHitDist_sim_py[iPMT] = MakeTH1('D', Form("MAPMT/Crosstalk/Simulation/HitDistPMT%i_py_sim",iPMT), Form("Simulated distance between hits via pythagoras, PMT%i",iPMT), 15, -0.5, 14.5, "distance between 2 hits (pixel)");
00198                 }
00199         }
00200         //hCrTalkDist_allPMTs = MakeTH2('D', "MAPMT/Crosstalk/hCrTalkDist_allPMTs", "Distribution of crosstalk hits", 32, 0., 32., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00201         hitPMT = MakeTH1('I', "MAPMT/Crosstalk/HitPMT_2hits", "PMTs with 2 hits per event", 25, -0.5, 24.5, "PMTs with 2 hits per event");
00202         hitPMTmore = MakeTH1('I', "MAPMT/Crosstalk/HitPMT_!=2hits", "PMTs with != 2 hits per event", 25, -0.5, 24.5, "PMTs with != 2 hits per event");
00203         hEntries2dAll = MakeTH2('D', "MAPMT/Crosstalk/hEntries2dAll", "Number of hits per MAPMT channel for all PMTs", 40, 0., 40., 32, 0., 32., "x-pos (pixel)", "y-pos (pixel)");
00204    }
00205    
00206    // read in individual data of gain for each pixel and pmt
00207    ReadGainPerPixel();
00208    
00209    
00210 
00211    SetParticleIdentificationConditions();
00212    
00213    // book conditions
00214    fMAPMT_cond = MakeWinCond("hADCall", -1000., 10000., "hADCall"); //wide, basically no cut
00215    fMAPMT_cond->SetHistogram("hADCall");
00216    fMAPMT_tcond = MakeWinCond("hHittime", 600., 800., "hHittime"); // cut on beam events only
00217    fMAPMT_tcond->SetHistogram("hHittime");
00218 
00219    Int_t nofPics = 30;
00220    int maxNofRings = 2;
00221    fCirclePicture.resize(nofPics);
00222    fCircle.resize(nofPics);
00223    fCircleText.resize(nofPics);
00224    fEllipse.resize(nofPics);
00225    fEllipseText.resize(nofPics);
00226    hSED.resize(nofPics);
00227    fRingHits.resize(nofPics);
00228    fEventHits.resize(nofPics);
00229    for (Int_t i = 0; i < nofPics; i++){
00230 
00231       hSED[i] = MakeTH2('I', Form("MAPMT/SED/SED%i",i),Form("Single Event %i",i), 32, 0., 208., 32, 0., 208., "x-pos (mm)", "y-pos (mm)");
00232 
00233       fCirclePicture[i] = new TGo4Picture(Form("Event%i",i), "Circle picture");
00234 
00235       fCircle[i].resize(maxNofRings);
00236       fEllipse[i].resize(maxNofRings);
00237       for (int iR = 0; iR < maxNofRings; iR++){
00238                   fCircle[i][iR] = new TEllipse(16., 16., 10.);
00239                   fCircle[i][iR]->SetFillStyle(0);
00240                   fCircle[i][iR]->SetLineWidth(2);
00241 
00242                   fEllipse[i][iR] = new TEllipse(16., 16., 10.);
00243                   fEllipse[i][iR]->SetFillStyle(0);
00244                   fEllipse[i][iR]->SetLineWidth(2);
00245                   fEllipse[i][iR]->SetLineColor(kBlue);
00246 
00247               fCirclePicture[i]->AddSpecialObject(fCircle[i][iR]);
00248               fCirclePicture[i]->AddSpecialObject(fEllipse[i][iR]);
00249       }
00250 
00251       fCircleText[i] = new TLatex(5, 30, "(x,y,r)");
00252       fEllipseText[i] = new TLatex(5, 30, "(x,y,r)");
00253 
00254       fCirclePicture[i]->AddObject(hSED[i]);
00255       fCirclePicture[i]->AddSpecialObject(fCircleText[i]);
00256       //fCirclePicture[i]->AddSpecialObject(fCircle[i]);
00257       fCirclePicture[i]->AddSpecialObject(fEllipseText[i]);
00258       //fCirclePicture[i]->AddSpecialObject(fEllipse[i]);
00259 
00260 
00261       fRingHits[i].resize(fMaxNofHitsInEventDraw);
00262       fEventHits[i].resize(fMaxNofHitsInEventDraw);
00263       for (int iH = 0; iH < fMaxNofHitsInEventDraw; iH++){
00264           fEventHits[i][iH] = new TEllipse(16., 16., 10.);
00265           fEventHits[i][iH]->SetFillColor(kRed);
00266           fCirclePicture[i]->AddSpecialObject(fEventHits[i][iH]);
00267       }
00268       for (int iH = 0; iH < fMaxNofHitsInEventDraw; iH++){
00269           fRingHits[i][iH] = new TEllipse(16., 16., 10.);;
00270           fRingHits[i][iH]->SetFillColor(kYellow+2);
00271           fCirclePicture[i]->AddSpecialObject(fRingHits[i][iH]);
00272       }
00273 
00274       AddPicture(fCirclePicture[i], "CirclePicture");
00275    }
00276 
00277    TGo4Log::Info("RICH Histograms created");
00278    
00279    fRingFitter = new CbmRichRingFitterCOP();
00280    fEllipseFitter = new CbmRichRingFitterEllipseTau();
00281    fRingFinder = new CbmRichRingFinderHoughImpl();
00282    
00283    fEventTime=-1;
00284    fO2level=-1;
00285    fH2Olevel=-1;
00286 
00287   if (GetPicture("Light_pulser_SPspectra")==0) {
00288     TGo4Picture* pic = new TGo4Picture("Light_pulser_SPspectra", "LED measurements - all single photon spectra");
00289     pic->SetDivision(5, 5);
00290     for(int a=0; a<=4; a++) {
00291         for(int b=0; b<=4; b++) {
00292                 if(a==4 && b==4) continue;
00293                 pic->Pic(a, b)->SetAutoScale(kTRUE);
00294                 pic->Pic(a, b)->AddObject(hSinglePhotonSpectra[5*a+b+1]);
00295                 }
00296         }    
00297     AddPicture(pic,"RICH");
00298   }
00299   
00300   for(int i=1; i<=24; i++) {
00301         if (GetPicture(Form("Crosstalk_pmt%i",i))==0) {
00302         pic_ct[i] = new TGo4Picture(Form("Crosstalk_pmt%i",i), Form("Crosstalk from data and simulation, pmt%i",i));
00303         pic_ct[i]->SetAutoScale(kTRUE);
00304         pic_ct[i]->AddObject(hHitDist_py[i]);
00305         pic_ct[i]->AddObject(hHitDist_sim_py[i]);
00306                 AddPicture(pic_ct[i],"RICH/Crosstalk_py");
00307         }  
00308   
00309         if (GetPicture(Form("Crosstalk_zones_pmt%i",i))==0) {
00310         pic_ct_zones[i] = new TGo4Picture(Form("Crosstalk_zones_pmt%i",i), Form("Crosstalk from data and simulation, pmt%i",i));
00311         pic_ct_zones[i]->SetAutoScale(kTRUE);
00312         pic_ct_zones[i]->AddObject(hHitDist[i]);
00313         pic_ct_zones[i]->AddObject(hHitDist_sim[i]);
00314                 AddPicture(pic_ct_zones[i],"RICH/Crosstalk_zones");
00315         } 
00316   }
00317 
00318 }
00319 
00320 TRICHProc::~TRICHProc()
00321 {
00322         if (fPar->AnalyseCrosstalk) {
00323            cout << "==> CrTalkResults() start" << endl;
00324       CrTalkResults(); // destructor is called when "shutdown analysis" is pressed
00325       cout << "==> CrTalkResults() end" << endl;
00326         }
00327         for(int i=1; i<=24; i++) {
00328                 Double_t temp = 1.0*crosstalk[i][1]/crosstalk[i][0];
00329                 cout << "PMT: " << i << "\t all hits: " << crosstalk[i][0] << "\t crosstalk hits: " << crosstalk[i][1] << "\t fraction: " << temp << endl;
00330         }
00331 }
00332 
00333 void TRICHProc::SetParticleIdentificationConditions()
00334 {
00335         //runs 31,32,33
00336         //Double_t electronpntsc1c2[4][2] = { {1540, 400}, {1410, 1660}, {750, 1690}, {720, 415} };
00337         
00338         //runs 18,22,27
00339         Double_t electronpntsc1c2[4][2] = { {1870, 550}, {1870, 2400}, {910, 2400}, {970, 550} };
00340         fElectronCondc1c2 = MakePolyCond("RichElectronCondc1c2", 4, electronpntsc1c2);
00341         fElectronCondc1c2->SetHistogram("Ch1_Ch2");
00342 
00343         Double_t muonpntsc1c2[4][2] = { {500, 500}, {450, 1300}, {1200, 1620}, {1300, 500}  };
00344         fMuonCondc1c2 = MakePolyCond("RichMuonCondc1c2", 4, muonpntsc1c2);
00345         fMuonCondc1c2->SetHistogram("Ch1_Ch2");
00346 
00347         Double_t pionpntsc1c2[4][2] = { {500, 500}, {450, 1300}, {1200, 1620}, {1300, 500} };
00348         fPionCondc1c2 = MakePolyCond("RichPionCondc1c2", 4, pionpntsc1c2);
00349         fPionCondc1c2->SetHistogram("Ch1_Ch2");
00350 }
00351 
00352 
00353 bool TRICHProc::ReadMapmtGeometry(const char* fname)
00354 {
00355   for (int ncell=0; ncell<NUM_MAPMP_CELLS; ncell++)
00356   {
00357     fMAPMTCells[ncell].xbin = -1.;
00358     fMAPMTCells[ncell].ybin = -1.;
00359     fMAPMTCells[ncell].xcoord = 0.;
00360     fMAPMTCells[ncell].ycoord = 0.;
00361     fMAPMTCells[ncell].globalchannel = 0;
00362     fMAPMTCells[ncell].mapmt = 0;
00363     fMAPMTCells[ncell].pixel = 0;
00364     fMAPMTCells[ncell].roc = 0;
00365     fMAPMTCells[ncell].feb = 0;
00366     fMAPMTCells[ncell].febchannel = 0;
00367     fMAPMTCells[ncell].ncell = 0;
00368 
00369     fMAPMTCells[ncell].trb = 0;
00370     fMAPMTCells[ncell].trbchannel = 0;
00371   }
00372   for (int trbindex=0; trbindex<TRB_TDC3_NUMBOARDS*64; trbindex++)
00373     fTRBindex[trbindex]=0;
00374 
00375   std::ifstream f(fname);
00376   if (!f)  {
00377     TGo4Log::Info("MAPMT geometry file not found %s", fname);
00378     return false;
00379   }
00380   
00381   char sbuf[1024];
00382   
00383   // skip first two lines
00384   f.getline(sbuf, sizeof(sbuf));
00385   f.getline(sbuf, sizeof(sbuf));
00386   
00387   int nread = 0;
00388   
00389   while (!f.eof()) {
00390     f.getline(sbuf, sizeof(sbuf));
00391     if (strlen(sbuf)==0) break;
00392     
00393     int globalchannel(0), mapmt(0), pixel(0), roc(0), feb(0), febchannel(0), trb(0), trbchannel(0);
00394     float xcoord(0.), ycoord(0.), xbin(0.), ybin(0.);
00395     
00396     if (sscanf(sbuf, "%d %d %d %d %d %d %f %f %f %f %d %d",
00397         &globalchannel, &roc, &feb, &febchannel, &mapmt, &pixel, &xcoord, &ycoord, &xbin, &ybin, &trb, &trbchannel)!=12)
00398     {
00399       TGo4Log::Error("MAPMT read line %d format error %s", nread, sbuf);
00400       return false;
00401     }
00402     
00403     int ncell = roc*256 + febchannel + (feb==2 ? 128 : 0);    
00404     
00405     if ((ncell<0) || (ncell>=NUM_MAPMP_CELLS))
00406     {
00407       //printf ("roc=%d, feb=%d, febchannel=%d, ncell=%d\n", roc, feb, febchannel, ncell);
00408        TGo4Log::Error("MAPMT read wrong cell address %d in line %d", ncell, nread);
00409       return false;
00410     }
00411     
00412     //printf("ncell = %d xbin = %3.1f ybin = %3.1f  trb %d \n", ncell, xbin, ybin, trb);
00413     
00414     fMAPMTCells[ncell].xbin = xbin;
00415     fMAPMTCells[ncell].ybin = ybin;
00416     fMAPMTCells[ncell].xcoord = xcoord;
00417     fMAPMTCells[ncell].ycoord = ycoord;
00418     fMAPMTCells[ncell].globalchannel = globalchannel;
00419     fMAPMTCells[ncell].mapmt = mapmt;
00420     fMAPMTCells[ncell].pixel = pixel;
00421     fMAPMTCells[ncell].roc = roc;
00422     fMAPMTCells[ncell].feb = feb;
00423     fMAPMTCells[ncell].febchannel = febchannel;
00424     fMAPMTCells[ncell].ncell = ncell;
00425 
00426     fMAPMTCells[ncell].trb = trb;
00427     // ignore trbchannel from file, instead calculate from febchannel
00428     trbchannel = 64- (febchannel%64);
00429     if ((trbchannel % 2)==1) trbchannel=trbchannel+1;
00430     else trbchannel=trbchannel-1;
00431     fMAPMTCells[ncell].trbchannel = trbchannel;
00432 
00433     if (trb > 0) {
00434       fTRBindex[trb*64+fMAPMTCells[ncell].trbchannel] = &(fMAPMTCells[ncell]);
00435        /*    printf("ncell %d mapmt %i, trb: %i   trbch: %i, xy: %f %f %f %f\n",ncell, fTRBindex[trb*64+trbchannel]->mapmt,
00436               fTRBindex[trb*64+trbchannel]->trb,
00437               fTRBindex[trb*64+trbchannel]->trbchannel,
00438               fTRBindex[trb*64+trbchannel]->xbin,
00439               fTRBindex[trb*64+trbchannel]->ybin,
00440               fMAPMTCells[ncell].xbin, fMAPMTCells[ncell].ybin); 
00441        */
00442 
00443     }
00444     nread++;
00445   }
00446   
00447   if (nread != NUM_MAPMP_CELLS)
00448   {
00449      TGo4Log::Info("MAPMT geometry read %d lines, expected %d. STILL PROCESSING!", nread, NUM_MAPMP_CELLS);
00450 //    return false;
00451   }
00452   
00453   TGo4Log::Info("Load MAPMT geometry from file %s done.", fname);
00454 
00455 /*
00456   printf("decoding array: \n");
00457   for (Int_t i=0; i<1280; i++) {
00458     printf("cell: %i   mapmt: %i,  febch: %i,  x: %f,   y: %f,  trb: %i,  trbch: %i \n",
00459            i,
00460            fMAPMTCells[i].mapmt,
00461            fMAPMTCells[i].febchannel,
00462            fMAPMTCells[i].xbin,
00463            fMAPMTCells[i].ybin,
00464            fMAPMTCells[i].trb,
00465            fMAPMTCells[i].trbchannel);
00466   }
00467 */
00468   return true;
00469 }
00470 
00471 
00472 void TRICHProc::ReadGainPerPixel() 
00473 {
00474         // read in fitted mean gain value for each pmt and pixel from file RICH/GainPerPixel.txt
00475         std::ifstream f("RICH/GainPerPixel.txt");
00476         if (!f) {
00477                 TGo4Log::Info("TRICHProc: File GainPerPixel.txt not found, setting all values to 1!");
00478                 for(int pmt=1; pmt<=24; pmt++){
00479                         for(int pixel=1; pixel<=64; pixel++) {
00480                                 gain_per_pixel[pmt][pixel] = 1;
00481                         }
00482                 }
00483         }
00484         else {
00485                 TGo4Log::Info("TRICHProc: File GainPerPixel.txt found and read!");
00486                 char fileinfo[256];
00487                 f.getline(fileinfo,256);
00488                 TGo4Log::Info(Form("TRICHProc: GainPerPixel.txt - %s", fileinfo));              
00489                 for(int pmt=1; pmt<=16; pmt++){
00490                         for(int pixel=1; pixel<=64; pixel++) {
00491                                 f >> gain_per_pixel[pmt][pixel];
00492                         }
00493                 }
00494                 for(int pmt=17; pmt<=24; pmt++){
00495                         for(int pixel=1; pixel<=16; pixel++) {
00496                                 f >> gain_per_pixel[pmt][pixel];
00497                         }
00498                 }
00499         }
00500         f.close();
00501         
00502         // calculate mean gean per pmt from gain per pixel
00503         for(int pmt=1; pmt<=16; pmt++){
00504                 Double_t temp = 0;
00505                 for(int pixel=1; pixel<=64; pixel++) {
00506                         temp += gain_per_pixel[pmt][pixel];
00507                 }
00508                 gain_mean[pmt] = temp/64;
00509         }
00510         for(int pmt=17; pmt<=24; pmt++){
00511                 Double_t temp = 0;
00512                 for(int pixel=1; pixel<=16; pixel++) {
00513                         temp += gain_per_pixel[pmt][pixel];
00514                 }
00515                 gain_mean[pmt] = temp/16;
00516         }
00517 
00518 
00519 }
00520 
00521 
00522 
00523 void TRICHProc::FinalizeEvent()
00524 {
00525   if (!fRocInputEvent->IsValid()) return; // skip incomplete roc events from first step
00526 
00527   // SL: do we need extra processing of pulser events here - probably not
00528   if (fMbsEvent && fMbsEvent->IsPulser())
00529   {
00531     ProcessLedPulser ();
00533     return;
00534   }
00535 
00537    // JAM: here is example how to access data from trbs:
00538 
00539    if (fTrbInputEvent && fTrbInputEvent->IsValid())
00540    {
00541      //printf(" found valid TRB input. \n"); 
00542       for (UInt_t RICHTRB = fPar->RICHFirstTRB; RICHTRB <= fPar->RICHLastTRB; RICHTRB++)
00543       {
00544          // Loop over all TRBs which deliver RICH data
00545 
00546          if (RICHTRB > TRB_TDC3_NUMBOARDS)
00547             continue;
00548 
00549          // printf("  looping over hits in TRB board %i \n",RICHTRB);
00550          TTrbData* theTRB = dynamic_cast<TTrbData*>(fTrbInputEvent->getEventElement(RICHTRB));
00551          if (theTRB == 0)
00552             continue; //no data from this TRB available;
00553 
00554           // JAM TODO TODO
00555 
00556          for (unsigned int i = 0; i<theTRB->fHits.size(); i++)
00557          {
00558             // TODO do something here with every hit
00559 
00560             // printf("***** Evaluate trb %u, hit %u, hits size %d\n", RICHTRB, i, theTRB->fHits.size());
00561            UShort_t tdc = theTRB->fHits[i].GetTDC();
00562            UShort_t ch = theTRB->fHits[i].GetChannel();
00563            UChar_t edge = theTRB->fHits[i].GetEdgeType();
00564            Double_t dtm = theTRB->fHits[i].GetDeltaT();
00565            // printf("     tdc %d, ch %d, edge %d  dt=%e\n", tdc, ch, edge, dtm);
00566            if (ch==0) continue; //skip entries from refergo4ence channel 0 !
00567             // TODO here mapping of trb/tdc/channel to mapmt pixels!
00568 
00569             if (fTRBindex[RICHTRB*64+ch] != 0)
00570             {
00571               // printf("trb: %i, trbch %i,   xbin,ybin: %f %f \n",RICHTRB, ch, fTRBindex[RICHTRB*64+ch]->xbin, fTRBindex[RICHTRB*64+ch]->ybin);
00572               hMAP_from_TRB->Fill (fTRBindex[RICHTRB*64+ch]->xbin, fTRBindex[RICHTRB*64+ch]->ybin);
00573               hEntries_TRBROC->Fill (fTRBindex[RICHTRB*64+ch]->xbin, fTRBindex[RICHTRB*64+ch]->ybin);
00574 
00575             }
00576 
00577          }
00578       } // trb
00579    }
00581 
00582    ProcessRICH();
00583    fEventNr++;  
00584 }
00585 
00586 
00587 void TRICHProc::InitEvent(TGo4EventElement* outevnt)
00588 {
00589   // first assign input event:
00590   // since input event object is never discarded within processor lifetime,
00591   // we just search for subevent by name once to speed up processing  
00592    if (fRocInputEvent==0) {
00593       TCBMBeamtimeEvent* btevent=dynamic_cast<TCBMBeamtimeEvent*>(GetInputEvent());
00594       if (btevent)
00595       {
00596          fRocInputEvent = dynamic_cast<TRocEvent*>(btevent->GetSubEvent("ROC"));
00597          fTrbInputEvent = dynamic_cast<TTrbEvent*>(btevent->GetSubEvent("TRB"));
00598          fTriglogEvent = dynamic_cast<TTriglogEvent*>(btevent->GetSubEvent("TRIGLOG"));
00599          fMbsEvent = dynamic_cast<TMbsCrateEvent*>(btevent->GetSubEvent("MBSCRATE"));
00600          fEpicsEvent = dynamic_cast<TEpicsEvent*>(btevent->GetSubEvent("EPICS"));
00601       } 
00602       else
00603       {
00604          fRocInputEvent = dynamic_cast<TRocEvent*>(GetInputEvent());
00605       }
00606       if (fRocInputEvent==0) {
00607          GO4_STOP_ANALYSIS_MESSAGE("**** TRICHProc: Fatal error: input event is not a TRocEvent!!! STOP GO4");
00608       }
00609    }
00610 
00611    // then assign output event
00612    // since output event object is never discarded within processor lifetime,
00613    // we just search for subevent by name once to speed up processing
00614    if (fOutputEvent==0)
00615    {
00616       TCBMBeamtimeEvent* btevent=dynamic_cast<TCBMBeamtimeEvent*>(outevnt);
00617       if (btevent)
00618       {
00619          fOutputEvent = dynamic_cast<TRICHEvent*>(btevent->GetSubEvent("RICH"));
00620          fHodo1 = dynamic_cast<TFiberHodEvent*>(btevent->GetSubEvent("Hodo1"));
00621          fHodo2 = dynamic_cast<TFiberHodEvent*>(btevent->GetSubEvent("Hodo2"));
00622          fBeamEvent = dynamic_cast<TBeamMonitorEvent*>(btevent->GetSubEvent("BEAM"));
00623       }
00624       else
00625       {
00626          fOutputEvent = dynamic_cast<TRICHEvent*>(outevnt);
00627       }
00628       if (fOutputEvent==0) {
00629          GO4_STOP_ANALYSIS_MESSAGE("**** TRICHProc: Fatal error: output event is not a TRICHEvent!!! STOP GO4");
00630       }
00631    }
00632 }
00633 
00634 
00635 void TRICHProc::ProcessRICH()
00636 {
00637    // Get Event time (in secs) from MBS event
00638    if (fTriglogEvent) {
00639       fEventTime = fTriglogEvent->fMbsTimeSecs;
00640    }
00641    
00642    // Get O2 and H2O concentration from Epics event, interpolate missing data
00643    if (fEpicsEvent && fEpicsEvent->IsValid()) {
00644       fO2level=fEpicsEvent->GetDouble("CBM:RICH:Gas:O2");
00645       fH2Olevel=fEpicsEvent->GetDouble("CBM:RICH:Gas:H2O");
00646       fTT1=fEpicsEvent->GetDouble("CBM:RICH:Gas:TT-1");
00647       fPTB=fEpicsEvent->GetDouble("CBM:RICH:Gas:PTB");
00648       fPress0=fEpicsEvent->GetDouble("CBM:RICH:press_0");
00649       fTemp0=fEpicsEvent->GetDouble("CBM:RICH:temp_0");
00650 
00651    }
00652    if (fEventTime>1318648546 && fEventTime<1318720000) {
00653       fO2level=350;
00654       if (fEventTime>1318655338) fO2level=18844646*pow(0.9999834,(fEventTime-1318000000));
00655    }
00656 
00657    //SaveParameters();
00658       
00659    // Get Cherenkov 1/2 signal
00660    Double_t cherenkov1=fMbsEvent->fData1182[1][0]; //RICH1+2: Electron
00661    Double_t cherenkov2=fMbsEvent->fData1182[0][1]; //RICH2only: Muon
00662                                                    //none: Pion
00663    Double_t leadglass=fMbsEvent->fData1182[1][5];
00664    
00665    //perform particle identification
00666    //Bool_t isElectron = fBeamEvent->fIsElectron;
00667    //Bool_t isMuon = fBeamEvent->fIsMuon;
00668    //Bool_t isPion = fBeamEvent->fIsPion;
00669    bool isElectron = fElectronCondc1c2->Test(cherenkov1,cherenkov2);
00670    bool isMuon = fMuonCondc1c2->Test(cherenkov1,cherenkov2);
00671    bool isPion = fPionCondc1c2->Test(cherenkov1,cherenkov2);
00672 
00673    //access information from fiber hodoscope
00674    Double_t hodx=0;
00675    Double_t hody=0;
00676    if (fHodo1 && (fHodo1->NumHits()>0)) {
00677       hodx = fHodo1->Hit(0).X;
00678       hody = fHodo1->Hit(0).Y;
00679    }
00680    //cout << "hodo x,y " << hodx << " " << hody << endl;
00681    
00682     //printf("\nEventtime in secs: %i \n",fEventTime);
00683     //printf("O2 content: %f \n",fO2level);
00684     //printf("H2O content: %f \n",fH2Olevel);
00685     //printf("Cherenkov 1/2: %f %f \n",cherenkov1,cherenkov2);
00686     //printf("Hodoscope position: %f / %f \n",hodx,hody);
00687 
00688    Int_t NrHitsperPMT[17];
00689    for (Int_t iPMT=0; iPMT<=16; iPMT++) NrHitsperPMT[iPMT]=0;
00690   
00691    vector<Double_t> timearray;
00692    //collect all hits from one event in vector
00693    vector<CbmRichHitLight> hits;
00694   
00695    for (UInt_t RICHROC=fPar->RICHFirstROC; RICHROC<=fPar->RICHLastROC; RICHROC++) {
00696      // cout << "  Processing ROC "<< RICHROC << endl;
00697       // Loop over all ROCs which deliver RICH data
00698       
00699       if (RICHROC>MAX_ROC) continue;
00700       
00701       TRocData* RocData=dynamic_cast<TRocData*>(fRocInputEvent->getEventElement(RICHROC));
00702       if (RocData==0) continue;  //no data from this ROC available;
00703       
00704       for (UInt_t MessageNr=0; MessageNr<RocData->fExtMessages.size(); MessageNr++) {
00705          //cout << "  Processing message " << MessageNr << "from ROC " << RICHROC << endl;
00706          //cout << "Msg.size() = " << RocData->fExtMessages.size() << endl;
00707          // Loop over all messages from particular RICH-ROC
00708          
00709          TRocMessageExtended msg=RocData->fExtMessages.at(MessageNr);
00710          if (msg.GetMessageType() != roc::MSG_HIT) continue;  // only interrested in HIT messages
00711          //cout << "roc::MSG_HIT" << endl;
00712 
00713          // extraxt message information
00714          UInt_t dataROC=msg.GetRocNumber();
00715          UInt_t dataFEBNr=msg.GetNxNumber();  // 0 or 2
00716          UInt_t dataFEBCh=msg.GetNxChNum();   // 0..127
00717          Double_t dataHitTime=msg.GetTriggerDeltaT();
00718          Int_t dataHitADCraw=msg.GetNxADC();
00719          Int_t dataHitADCcor=msg.GetCorrectedNxADC();
00720                   
00721          // calculate ncell from ROC and FEB number
00722          int ncell = (dataROC - fPar->RICHFirstROC)*256 + dataFEBCh + (dataFEBNr==2 ? 128 : 0);
00723 
00724          Float_t dist=msg.GetTriggerDeltaT(); // corrected trigger diff is already in input message!
00725          hHittime->Fill(dist);
00726          
00727          // calculate corrected (inverse) ADC value, mapmtPedestal can be changed in TRICHParams.cxx or parameter editor in the Go4 GUI
00728          int rADC = fPar->mapmtPedestal - msg.GetNxADC();
00729          
00730          // fill histogram with all MAPMT ADC values to have reference for condition
00731          hADCall->Fill(dataHitADCcor);
00732          
00733          // get mapping from geometry file
00734          double xbin = fMAPMTCells[ncell].xbin;
00735          double ybin = fMAPMTCells[ncell].ybin;
00736          double xcoord = fMAPMTCells[ncell].xcoord;
00737          double ycoord = fMAPMTCells[ncell].ycoord;
00738          Int_t mapmt = fMAPMTCells[ncell].mapmt;
00739          
00740          // store mapped info to ROOT tree
00741          fOutputEvent->fMAPM_X[ncell]=xbin;
00742          fOutputEvent->fMAPM_Y[ncell]=ybin;
00743          fOutputEvent->fMAPM_Xcoord[ncell]=xcoord;
00744          fOutputEvent->fMAPM_Ycoord[ncell]=ycoord;
00745          fOutputEvent->fMAPM_Integral[ncell]=dataHitADCcor;
00746          fOutputEvent->fMAPM_Integral_raw[ncell]=rADC;
00747          fOutputEvent->fMAPM_DeltaT[ncell]=dist;
00748          
00749          // printf("GET HIT ncell = %d xbin = %3.1f ybin = %3.1f\n", ncell, xbin, ybin);
00750          
00751          // if condition is ok, fill XY histogram.
00752          // Condition range can be modified in the condition editor in the Go4 GUI
00753          if (fMAPMT_cond->Test(dataHitADCcor) && fMAPMT_tcond->Test(dist)){
00754                  timearray.push_back(dist);
00755 
00756                  // create hit and add it to the vector of hits
00757                  CbmRichHitLight hit(xcoord, ycoord);
00758                  hits.push_back(hit);
00759                  
00760                  if (isElectron) hEntries2d[0]->Fill(xbin, ybin);
00761                  if (isMuon) hEntries2d[2]->Fill(xbin, ybin);
00762                  if (isPion) hEntries2d[1]->Fill(xbin, ybin);
00763                  hEntries2d[3]->Fill(xbin, ybin);
00764                  hEntries_TRBROC->Fill(xbin,ybin);
00765 
00766                  Double_t c1=-0.050;
00767                  Double_t c2=-0.040;
00768                  Double_t c3=-0.030;
00769                  Double_t c4=-0.020;
00770                  Double_t c5=-0.010;
00771                  hMAPMT_hcor2d_s1->Fill(xbin+c1*hodx,ybin+0*hody);
00772                  hMAPMT_hcor2d_s2->Fill(xbin+c2*hodx,ybin+0*hody);
00773                  hMAPMT_hcor2d_s3->Fill(xbin+c3*hodx,ybin+0*hody);
00774                  hMAPMT_hcor2d_s4->Fill(xbin+c4*hodx,ybin+0*hody);
00775                  hMAPMT_hcor2d_s5->Fill(xbin+c5*hodx,ybin+0*hody);
00776 
00777                  hEntriesperPMT->Fill(mapmt);
00778                  hEntriesperPMT2D->Fill(xbin, ybin);
00779                  hEntries2dIntegral->Fill(xbin, ybin, dataHitADCcor);
00780                  
00781                  if (mapmt >=1 && mapmt <= 16) NrHitsperPMT[mapmt]++;
00782          }
00783          
00784          int nbin = hEntries2d[3]->GetBin(hEntries2d[3]->GetXaxis()->FindBin(xbin), hEntries2d[3]->GetYaxis()->FindBin(ybin));
00785          double entries = hEntries2d[3]->GetBinContent(nbin);
00786          double integral = hEntries2dIntegral->GetBinContent(nbin);
00787          hEntries2dMean->SetBinContent(nbin, (entries <= 0.) ? 0. : integral/entries);
00788       } //MessageNr
00789    } //RICHROC
00790    
00791    
00792    
00793    
00794       
00795    
00796       
00797    
00798    
00799    
00800    
00801    
00802    
00803    
00804    
00805 
00806    for (Int_t iPMT=1; iPMT<=16; iPMT++) {
00807       hHitMultperPMT[iPMT]->Fill(NrHitsperPMT[iPMT]);
00808    }
00809    
00810    // time deviation of individual hits from average, not very fast implementation...
00811    Double_t mean=0;
00812    for (Int_t ihit=0; ihit<timearray.size(); ihit++){
00813       // printf("time %i: %f \n",ihit,timearray[ihit]);
00814       mean=mean+timearray[ihit];
00815    }
00816    mean=mean/timearray.size();
00817    // printf("mean: %f \n",mean);
00818    for (Int_t ihit=0; ihit<timearray.size(); ihit++) {
00819       hTimespread->Fill(timearray[ihit]-mean);
00820    }
00821 
00822    if (fHodo1 != NULL && fHodo1->NumHits() >=1){
00823            hNofHitsInHodo->Fill(fHodo1->NumHits());
00824            hHodoADCvsNofHitInEv->Fill(fHodo1->Hit(0).ADC, hits.size());
00825    }
00826 
00827    DoAnalysis(hits, isElectron, isPion, isMuon);
00828 }
00829 
00830 
00831 
00832 // ########################## crosstalk results start ############################
00833 
00834 void TRICHProc::CrTalkResults()
00835 {
00836   Double_t ratioNN[25]; // ratio between hits in next neighbours and all hits
00837   for (Int_t iPMT=1; iPMT<=12; iPMT++) {
00838     ratioNN[iPMT]=hHitDist[iPMT]->GetBinContent(2)/hHitDist[iPMT]->GetEntries()-0.1025; // normalisation to all bin entries, H8500: 0.1025, R11265: 0.3281
00839     //ratioNN[iPMT]=hHitDist[iPMT]->GetBinContent(2)/(hHitDist[iPMT]->GetEntries()-hHitDist[iPMT]->GetBinContent(2))-0.11421;  // normalisation to all bin entries w/o bin2
00840     //ratioNN[iPMT]=hHitDist[iPMT]->GetBinContent(2)/hHitDist[iPMT]->GetBinContent(5)-0.476455;
00841     cout << "ratioNN[" << iPMT << "] = " << ratioNN[iPMT] << endl;
00842   }
00843   for (Int_t iPMT=13; iPMT<=16; iPMT++) {
00844     ratioNN[iPMT]=hHitDist[iPMT]->GetBinContent(2)/hHitDist[iPMT]->GetEntries()-0.1025; // normalisation to all bin entries, H8500: 0.1025, R11265: 0.3281
00845     //ratioNN[iPMT]=hHitDist[iPMT]->GetBinContent(2)/(hHitDist[iPMT]->GetEntries()-hHitDist[iPMT]->GetBinContent(2))-0.11421;  // normalisation to all bin entries w/o bin2
00846     //ratioNN[iPMT]=hHitDist[iPMT]->GetBinContent(2)/hHitDist[iPMT]->GetBinContent(5)-0.476455;
00847     cout << "ratioNN[" << iPMT << "] = " << ratioNN[iPMT] << endl;
00848   }
00849   for (Int_t iPMT=17; iPMT<=24; iPMT++) {
00850     ratioNN[iPMT]=hHitDist[iPMT]->GetBinContent(2)/hHitDist[iPMT]->GetEntries()-0.3281; // normalisation to all bin entries, H8500: 0.1025, R11265: 0.3281
00851     //ratioNN[iPMT]=hHitDist[iPMT]->GetBinContent(2)/(hHitDist[iPMT]->GetEntries()-hHitDist[iPMT]->GetBinContent(2))-0.11421;  // normalisation to all bin entries w/o bin2
00852     //ratioNN[iPMT]=hHitDist[iPMT]->GetBinContent(2)/hHitDist[iPMT]->GetBinContent(5)-0.476455;
00853     cout << "ratioNN[" << iPMT << "] = " << ratioNN[iPMT] << endl;
00854   }
00855   
00856   cout << "Number of events with 2 hits: " << nof2hitevents << endl;
00857 }
00858 
00859 // ########################## crosstalk results end ############################## 
00860 
00861 
00862 
00863 
00864 void TRICHProc::SaveParameters()
00865 {
00866         ofstream outfile;
00867         outfile.open ("rich_parameters_130404.txt", ios_base::app);
00868 
00869         double o2 = -1.0;
00870         double h2o = -1.0;
00871         double ptb = -1.0;
00872         double tt1 = -1.0;
00873         double pt4 = -1.0;
00874         double refrIndex = -1.0;
00875         double press0 = -1.0;
00876     double temp0 = -1.0;
00877 
00878         if (fEpicsEvent && fEpicsEvent->IsValid()) {
00879         o2 = fEpicsEvent->GetDouble("CBM:RICH:Gas:O2");
00880                 h2o = fEpicsEvent->GetDouble("CBM:RICH:Gas:H2O");
00881         ptb = fEpicsEvent->GetDouble("CBM:RICH:Gas:PTB");
00882                 tt1 = fEpicsEvent->GetDouble("CBM:RICH:Gas:TT-1");
00883         pt4 = fEpicsEvent->GetDouble("CBM:RICH:Gas:PT-4");
00884         refrIndex = fEpicsEvent->GetDouble("CBM:RICH:Gas:RefrIndex");
00885                 press0 = fEpicsEvent->GetDouble("CBM:RICH:press_0");
00886         temp0 = fEpicsEvent->GetDouble("CBM:RICH:temp_0");
00887         };
00888 
00889         if (o2 == -1. && h2o == -1. && ptb == -1. && tt1 == -1. && pt4 == -1. && refrIndex == -1. && press0 == -1 && temp0 == -1) return;
00890 
00891         TGo4AnalysisStep* firststep = TGo4Analysis::Instance()->GetAnalysisStepNum(0);
00892         if (dynamic_cast<TGo4MbsFileParameter*>(firststep->GetEventSource())) {
00893                 outfile << firststep->GetEventSourceName() << "\t" <<
00894                                 fEventTime << "\t" <<
00895                                 o2 << "\t" <<
00896                                 h2o << "\t" <<
00897                                 ptb << "\t" <<
00898                                 tt1 << "\t" <<
00899                                 pt4 << "\t" <<
00900                                 refrIndex << "\t" << 
00901                                 press0 << "\t" << 
00902                                 temp0 << endl;
00903         }
00904         outfile.close();
00905         time_t rawtime = fEventTime;
00906         //time ( &rawtime );
00907         printf ( "The current local time is: %s", ctime (&rawtime) );
00908 }
00909 
00910 void TRICHProc::ProcessLedPulser()
00911 {
00912 
00913    UInt_t hitsPerPMT[24] = {0, 0, 0, 0, 0, 0,
00914                             0, 0, 0, 0, 0, 0,
00915                             0, 0, 0, 0, 0, 0,
00916                             0, 0, 0, 0, 0, 0};  // hardcoded number of mamps
00917                             
00918    Int_t hitArrayPerPMT[25][9][9] = {0}; // array for each pmt and for each pixel, filled with number of hits (=0,1,...)
00919    Int_t hitArrayPerPMT_ADC[25][9][9] = {0};  // array for each pmt and for each pixel, filled with ADC value of hits                        
00920 
00921    for (UInt_t RICHROC = fPar->RICHFirstROC; RICHROC <= fPar->RICHLastROC; RICHROC++)
00922    {
00923       // cout << "  Processing ROC "<< RICHROC << endl;
00924       // Loop over all ROCs which deliver RICH data
00925       if (RICHROC>MAX_ROC) continue;
00926 
00927       TRocData* RocData=dynamic_cast<TRocData*>(fRocInputEvent->getEventElement(RICHROC));
00928       if (RocData==0) continue;  //no data from this ROC available;
00929 
00930       for (UInt_t MessageNr = 0; MessageNr < RocData->fExtMessages.size(); MessageNr++)
00931       {
00932          // cout << "  Processing message " << MessageNr << "from ROC " << RICHROC << endl
00933          // Loop over all messages from particular RICH-ROC
00934 
00935          TRocMessageExtended msg = RocData->fExtMessages.at(MessageNr);
00936          if (msg.GetMessageType() != roc::MSG_HIT) continue;  // only interrested in HIT messages
00937 
00938          // extraxt message information
00939          UInt_t dataROC = msg.GetRocNumber();
00940          UInt_t dataFEBNr = msg.GetNxNumber();  // 0 or 2
00941          UInt_t dataFEBCh = msg.GetNxChNum();   // 0..127
00942          Double_t dataHitTime = msg.GetTriggerDeltaT();
00943          Int_t dataHitADCraw = msg.GetNxADC();
00944          Int_t dataHitADCcor = msg.GetCorrectedNxADC();
00945 
00946          // calculate ncell from ROC and FEB number
00947          int ncell = (dataROC - fPar->RICHFirstROC)*256 + dataFEBCh + (dataFEBNr==2 ? 128 : 0);
00948 
00949          // get mapping from geometry file
00950          double xbin = fMAPMTCells[ncell].xbin;
00951          double ybin = fMAPMTCells[ncell].ybin;
00952          double xcoord = fMAPMTCells[ncell].xcoord;
00953          double ycoord = fMAPMTCells[ncell].ycoord;
00954          Int_t mapmt = fMAPMTCells[ncell].mapmt;
00955          Int_t pixel = fMAPMTCells[ncell].pixel;
00956          Double_t dataHitNofPh = 1.0*dataHitADCcor/gain_per_pixel[mapmt][pixel];
00957 
00958          if (xbin == -1. && ybin == -1.)
00959          {
00960             // printf ("not connected channel: ncell=%d, dataROC=%d, dataFEBCh=%d, dataFEBNr=%d\n",
00961             //    ncell, dataROC, dataFEBCh, dataFEBNr);
00962             continue;
00963          }
00964 
00965          hLedPulserEntriesPerChannel->Fill(ncell);
00966          hLedPulserEntries2d->Fill(xbin, ybin);
00967          hLedPulserADCvsPixel[mapmt-1]->Fill(pixel, dataHitADCcor);
00968                  hLedPulserADCvsPixel[mapmt-1]->Fill(0., dataHitADCcor);
00969 
00970                  hSinglePhotonSpectra[mapmt]->Fill(dataHitADCcor);
00971                  hLEDPulserADC_nofph[mapmt]->Fill(dataHitNofPh);
00972                  hLEDPulserADC_nofphPP[mapmt]->Fill(pixel, dataHitNofPh);
00973                  
00974                 // if(mapmt == 12) {
00975                 //      cout << "PMT: " << mapmt << " , Pixel: " << pixel << "\t Korrektur: " << gain_per_pixel[mapmt][pixel] << endl;
00976                 // }    
00977 
00978          hitsPerPMT[mapmt-1]++;
00979          
00980                  Int_t pixel_x = (pixel-1)%8+1;
00981                  Int_t pixel_y = (pixel-1)/8+1;
00982          hitArrayPerPMT[mapmt][pixel_x][pixel_y] += 1;
00983       //   hitArrayPerPMT_ADC[mapmt][pixel_x][pixel_y] = dataHitADCcor;
00984          hitArrayPerPMT_ADC[mapmt][pixel_x][pixel_y] = dataHitNofPh*gain_mean[mapmt]; //ADC value corrected with gain per pixel
00985          
00986          
00987          
00988          
00989          
00990 
00991       } // MessageNr
00992 
00993    } // RICHROC
00994 
00995 
00996 
00997         // ###### BEGIN - extended event analysis ################################################
00998         for (Int_t iPMT=1; iPMT<=24; iPMT++) {
00999                 if (hitsPerPMT[iPMT-1] >= 1 && fPar->ExtendedAnalysis) {
01000                         for (UInt_t RICHROC = fPar->RICHFirstROC; RICHROC <= fPar->RICHLastROC; RICHROC++) {
01001                         if (RICHROC>MAX_ROC) continue;
01002 
01003                         TRocData* RocData=dynamic_cast<TRocData*>(fRocInputEvent->getEventElement(RICHROC));
01004                         if (RocData==0) continue;  //no data from this ROC available;
01005 
01006                         for (UInt_t MessageNr = 0; MessageNr < RocData->fExtMessages.size(); MessageNr++) {
01007 
01008                                 TRocMessageExtended msg = RocData->fExtMessages.at(MessageNr);
01009                                 if (msg.GetMessageType() != roc::MSG_HIT) continue;  // only interrested in HIT messages
01010 
01011                                 // extraxt message information
01012                                 UInt_t dataROC = msg.GetRocNumber();
01013                                 UInt_t dataFEBNr = msg.GetNxNumber();  // 0 or 2
01014                                 UInt_t dataFEBCh = msg.GetNxChNum();   // 0..127
01015                                 Double_t dataHitTime = msg.GetTriggerDeltaT();
01016                                 Int_t dataHitADCraw = msg.GetNxADC();
01017                                 Int_t dataHitADCcor = msg.GetCorrectedNxADC();
01018 
01019                                 // calculate ncell from ROC and FEB number
01020                                 int ncell = (dataROC - fPar->RICHFirstROC)*256 + dataFEBCh + (dataFEBNr==2 ? 128 : 0);
01021 
01022                                 // get mapping from geometry file
01023                                 double xbin = fMAPMTCells[ncell].xbin;
01024                                 double ybin = fMAPMTCells[ncell].ybin;
01025                                 double xcoord = fMAPMTCells[ncell].xcoord;
01026                                 double ycoord = fMAPMTCells[ncell].ycoord;
01027                                 Int_t mapmt = fMAPMTCells[ncell].mapmt;
01028                                 Int_t pixel = fMAPMTCells[ncell].pixel;
01029                                 Double_t dataHitNofPh = 1.0*dataHitADCcor/gain_per_pixel[mapmt][pixel];
01030 
01031                                 if (xbin == -1. && ybin == -1.) {
01032                                     continue;
01033                                 }       
01034                                         if(mapmt == iPMT) {             
01035                                                 Int_t hitSum = 0;               
01036                                                 Int_t pixel_x = (pixel-1)%8+1;
01037                                                 Int_t pixel_y = (pixel-1)/8+1;
01038                                                 Bool_t ctHit = kFALSE;
01039                                                 for(int x=pixel_x-1; x<=pixel_x+1; x++) {
01040                                                         for(int y=pixel_y-1; y<=pixel_y+1; y++) {
01041                                                                 hitSum += hitArrayPerPMT[mapmt][x][y];
01042                                                         //      if(dataHitADCcor < 0.5 * hitArrayPerPMT_ADC[mapmt][x][y]) { // old version
01043                                                         //      if(dataHitNofPh < 0.1 * gain_per_pixel[mapmt][pixel]) { // < 10% of single-photon peak
01044                                                         //              ctHit = kTRUE;
01045                                                         //      }
01046                                                         }
01047                                                 }
01048                                                 
01049                                                 
01050                                                 if(dataHitADCcor < 0.1 * gain_per_pixel[mapmt][pixel] && hitSum > 1) { // < 10% of single-photon peak
01051                                                         ctHit = kTRUE;
01052                                                 }
01053                                         
01054                                                 if(hitSum == 1) { // only pixel with no hit in the 8 neighbouring pixel
01055                                                                 hLEDPulserADC_1hit[mapmt]->Fill(dataHitADCcor);
01056                                                                 hLEDPulserADC_nofph_mct[mapmt]->Fill(1.0*dataHitADCcor/gain_per_pixel[mapmt][pixel]);
01057                                                 //      if(pixel_x != 1 && pixel_x !=8 && pixel_y != 1 && pixel_y !=8) { // not a pixel from the edge
01058                                                                 hLEDPulserADC_1hit_perPixel[mapmt]->Fill(pixel,dataHitADCcor);
01059                                                         //      hLEDPulserADCcor_1hit[mapmt]->Fill(dataHitADCcor/gain_per_pixel[mapmt][pixel]*gain_mean[mapmt]);
01060                                                                 hLEDPulserADCcor_1hit[mapmt]->Fill(dataHitADCcor);
01061                                                 //      }
01062                                                 }
01063                                                 if(ctHit) {     
01064                                                         hLEDPulserADC_ctHits[mapmt]->Fill(dataHitADCcor);
01065                                                 }       
01066                                                 if(hitsPerPMT[mapmt-1] == 1) {
01067                                                         hLEDPulserADC_nofph_1hitpMA[mapmt]->Fill(1.0*dataHitADCcor/gain_per_pixel[mapmt][pixel]);
01068                                                 }
01069                                                         
01070                                                         
01071                                                 // test histograms:
01072                                         SP_ADC_all[mapmt]->Fill(dataHitADCcor); 
01073                                         SP_norm_all[mapmt]->Fill(dataHitNofPh); 
01074                                         if(hitSum == 1) {
01075                                                 SP_ADC_1hit[mapmt]->Fill(dataHitADCcor);
01076                                                 SP_norm_1hit[mapmt]->Fill(dataHitNofPh);
01077                                         }
01078                                         if(hitSum > 1 && dataHitADCcor < 0.1 * gain_per_pixel[mapmt][pixel]) {
01079                                                 SP_ADC_ct[mapmt]->Fill(dataHitADCcor);
01080                                                 SP_norm_ct[mapmt]->Fill(dataHitNofPh);
01081                                                 crosstalk[mapmt][1]++;
01082                                         }
01083                                         if(hitSum > 1 && dataHitADCcor > 0.1 * gain_per_pixel[mapmt][pixel]) {
01084                                                 SP_ADC_ctmain[mapmt]->Fill(dataHitADCcor);
01085                                                 SP_norm_ctmain[mapmt]->Fill(dataHitNofPh);
01086                                         }
01087                                         if(dataHitADCcor > 0.1 * gain_per_pixel[mapmt][pixel]) {
01088                                                 SP_ADC_noct[mapmt]->Fill(dataHitADCcor);
01089                                                 SP_norm_noct[mapmt]->Fill(dataHitNofPh);
01090                                         }
01091                                         
01092                                         crosstalk[mapmt][0]++; // number of all hits
01093                         
01094                                                                 
01095                                         }
01096                                         
01097                                 } // MessageNr
01098                    } // RICHROC
01099                 }
01100         }
01101                 
01102         
01103         // ###### END - extended event analysis ################################################        
01104                 
01105                 
01106                 
01107                 
01108    
01109         // ############################################ BEGIN CROSSTALK ANALYSIS (only for LED hits) ##############################
01110         if(fPar->AnalyseCrosstalk) {
01111                 Double_t hitDist;
01112                 Double_t hitDist_py;
01113                 Double_t xbinCrTalk[2];
01114                 Double_t ybinCrTalk[2];
01115                 Int_t NrHit=0;
01116                 xbinCrTalk[0]=-1;
01117                 xbinCrTalk[1]=-1;
01118                 ybinCrTalk[0]=-1;
01119                 ybinCrTalk[1]=-1;
01120                 int n_edge=0;
01121                 int n_corner=0;
01122                 int n_center=0;
01123                 
01124         
01125                 for (Int_t iPMT=1; iPMT<=24; iPMT++) {
01126                         hitDist=0;
01127                         hitDist_py=0;
01128                         NrHit=0;
01129                         xbinCrTalk[0]=-1;
01130                         xbinCrTalk[1]=-1;
01131                         ybinCrTalk[0]=-1;
01132                         ybinCrTalk[1]=-1;
01133                         
01134                         // select only events with 2 hits on 1 PMT
01135                         if (hitsPerPMT[iPMT-1]==2) {
01136                                 nof2hitevents++;
01137                                 hitPMT->Fill(iPMT);
01138                                 NrHit = 0;
01139                                 fEventNr_2ndloop++;
01140                         
01141                                 for (UInt_t RICHROC = fPar->RICHFirstROC; RICHROC <= fPar->RICHLastROC; RICHROC++) {
01142                                         if (RICHROC>MAX_ROC) continue;
01143 
01144                                         TRocData* RocData=dynamic_cast<TRocData*>(fRocInputEvent->getEventElement(RICHROC));
01145                                         if (RocData==0) continue;  //no data from this ROC available;
01146         
01147                                         for (UInt_t MessageNr = 0; MessageNr < RocData->fExtMessages.size(); MessageNr++) {
01148                                                 TRocMessageExtended msg = RocData->fExtMessages.at(MessageNr);
01149                                                 if (msg.GetMessageType() != roc::MSG_HIT) continue;  // only interrested in HIT messages
01150 
01151                                                 // extraxt message information
01152                                                 UInt_t dataROC = msg.GetRocNumber();
01153                                                 UInt_t dataFEBNr = msg.GetNxNumber();  // 0 or 2
01154                                                 UInt_t dataFEBCh = msg.GetNxChNum();   // 0..127
01155                                                 Double_t dataHitTime = msg.GetTriggerDeltaT();
01156                                                 Int_t dataHitADCraw = msg.GetNxADC();
01157                                                 Int_t dataHitADCcor = msg.GetCorrectedNxADC();
01158 
01159                                                 // calculate ncell from ROC and FEB number
01160                                                 int ncell = (dataROC - fPar->RICHFirstROC)*256 + dataFEBCh + (dataFEBNr==2 ? 128 : 0);
01161 
01162                                                 // get mapping from geometry file
01163                                                 double xbin = fMAPMTCells[ncell].xbin;
01164                                                 double ybin = fMAPMTCells[ncell].ybin;
01165                                                 Int_t mapmt = fMAPMTCells[ncell].mapmt;
01166                                                 Int_t pixel = fMAPMTCells[ncell].pixel;
01167 
01168                                         
01169                                                 if (mapmt != iPMT) continue;//iPMT
01170                                                 if (xbin == -1. && ybin == -1.) continue;
01171                                          
01172                                                 //if (NrHit>1) continue;
01173                                                 // if (fMAPMT_cond->Test(dataHitADCcor) && fMAPMT_tcond->Test(dataHitTime)) {
01174                                                 
01175                                                 xbinCrTalk[NrHit] = fMAPMTCells[ncell].xbin;
01176                                                 ybinCrTalk[NrHit] = fMAPMTCells[ncell].ybin;
01177                                                 
01178                                                 
01179                                                 if(iPMT == 18 || iPMT == 20 || iPMT == 22 || iPMT == 24) {
01180                                                         xbinCrTalk[NrHit] = fMAPMTCells[ncell].xcoord;
01181                                                         ybinCrTalk[NrHit] = fMAPMTCells[ncell].ycoord;
01182                                                 }
01183                                                 
01184                                         double xcoord = fMAPMTCells[ncell].xcoord;
01185                                         double ycoord = fMAPMTCells[ncell].ycoord;
01186                                                 
01187                                                 if(NrHit >= 1) continue;
01188                                                 NrHit++;
01189                                         } // MessageNr
01190                                 }
01191                                 
01192                         //######################### START - Optional simulation of crosstalk data #######################
01193                                 if(fPar->SimulateCrosstalk) {
01194                                   Int_t NofPixel; // 64 for H8500 and XP85012 (64 pixel), 16 for R11265 (16 pixel)
01195                                   if(iPMT <= 16) NofPixel = 64;
01196                                   if(iPMT > 16) NofPixel = 16;
01197                                   Int_t sqp = TMath::Sqrt(NofPixel);
01198                                 
01199                                   for(int i=1; i<=10; i++) { //simulate 10x as much events as in data
01200                                         Int_t pixel1 = (rand()%NofPixel)+1;
01201                                         Int_t pixel2 = (rand()%NofPixel)+1;
01202                         
01203                                         Int_t pixel1_x = (pixel1-1)%sqp;
01204                                         Int_t pixel1_y = (pixel1-1)/sqp;
01205                                         Int_t pixel2_x = (pixel2-1)%sqp;
01206                                         Int_t pixel2_y = (pixel2-1)/sqp;
01207                                         
01208                                         Int_t distance_x_sim = TMath::Abs(pixel1_x - pixel2_x);
01209                                         Int_t distance_y_sim = TMath::Abs(pixel1_y - pixel2_y);
01210                                                                 
01211                                         Int_t hitDist_sim = 12;
01212                                         Double_t hitDist_sim_py = 0;
01213                                         
01214                                         // Calculation of distance with separate areas                          
01215                                         if(distance_x_sim == 0 && distance_y_sim == 0) hitDist_sim = 0;
01216                                         if((distance_x_sim == 1 && distance_y_sim == 0) || (distance_x_sim == 0 && distance_y_sim == 1)) hitDist_sim = 1;
01217                                         if(distance_x_sim == 1 && distance_y_sim == 1) hitDist_sim = 2;
01218                                         if((distance_x_sim == 2 && distance_y_sim == 0) || (distance_x_sim == 0 && distance_y_sim == 2)) hitDist_sim = 3;
01219                                         if((distance_x_sim == 2 && distance_y_sim == 1) || (distance_x_sim == 1 && distance_y_sim == 2)) hitDist_sim = 4;
01220                                         if(distance_x_sim == 2 && distance_y_sim == 2) hitDist_sim = 5;
01221                                         if((distance_x_sim == 3 && distance_y_sim == 0) || (distance_x_sim == 0 && distance_y_sim == 3)) hitDist_sim = 6;
01222                                         if((distance_x_sim == 3 && distance_y_sim == 1) || (distance_x_sim == 1 && distance_y_sim == 3)) hitDist_sim = 7;
01223                                         if((distance_x_sim == 3 && distance_y_sim == 2) || (distance_x_sim == 2 && distance_y_sim == 3)) hitDist_sim = 8;
01224                                         if(distance_x_sim == 3 && distance_y_sim == 3) hitDist_sim = 9;
01225                                         
01226                                         // Calculation of distance with Pythagoras
01227                                         hitDist_sim_py = TMath::Sqrt(distance_x_sim*distance_x_sim + distance_y_sim*distance_y_sim);
01228                                         // Filling histograms   
01229                                         hHitDist_sim[iPMT]->Fill(hitDist_sim, 0.1);
01230                                         hHitDist_sim_py[iPMT]->Fill(hitDist_sim_py, 0.1);
01231                                   }
01232                                 }
01233                         //######################### END - Optional simulation of crosstalk data #######################
01234                                 
01235                         
01236                                 if (xbinCrTalk[0]!=-1 && xbinCrTalk[1]!=-1 && ybinCrTalk[0]!=-1 && ybinCrTalk[1]!=-1) {
01237                                         hEntries2dperPMT[iPMT]->Fill(xbinCrTalk[0],ybinCrTalk[0]);//iPMT
01238                                         hEntries2dperPMT[iPMT]->Fill(xbinCrTalk[1],ybinCrTalk[1]);//iPMT
01239                                 
01240                                         Int_t distance_x = TMath::Abs(xbinCrTalk[1]-xbinCrTalk[0]);
01241                                         Int_t distance_y = TMath::Abs(ybinCrTalk[1]-ybinCrTalk[0]);
01242 
01243                                         // Calculation of distance with separate areas
01244                                         if(distance_x == 0 && distance_y == 0) hitDist = 0;
01245                                         if((distance_x == 1 && distance_y == 0) || (distance_x == 0 && distance_y == 1)) hitDist = 1;
01246                                         if(distance_x == 1 && distance_y == 1) hitDist = 2;
01247                                         if((distance_x == 2 && distance_y == 0) || (distance_x == 0 && distance_y == 2)) hitDist = 3;
01248                                         if((distance_x == 2 && distance_y == 1) || (distance_x == 1 && distance_y == 2)) hitDist = 4;
01249                                         if(distance_x == 2 && distance_y == 2) hitDist = 5;
01250                                         if((distance_x == 3 && distance_y == 0) || (distance_x == 0 && distance_y == 3)) hitDist = 6;
01251                                         if((distance_x == 3 && distance_y == 1) || (distance_x == 1 && distance_y == 3)) hitDist = 7;
01252                                         if((distance_x == 3 && distance_y == 2) || (distance_x == 2 && distance_y == 3)) hitDist = 8;
01253                                         if(distance_x == 3 && distance_y == 3) hitDist = 9;
01254                                         
01255                                         
01256                                         if(iPMT == 18 || iPMT == 20 || iPMT == 22 || iPMT == 24) {
01257                                                 Double_t hitDist_mm = TMath::Sqrt((xbinCrTalk[1]-xbinCrTalk[0])*(xbinCrTalk[1]-xbinCrTalk[0]) + (ybinCrTalk[1]-ybinCrTalk[0])*(ybinCrTalk[1]-ybinCrTalk[0]));
01258                                                 if(hitDist_mm < 25) hitDist = 9;
01259                                                 if(hitDist_mm < 22) hitDist = 8;
01260                                                 if(hitDist_mm < 19) hitDist = 7;
01261                                                 if(hitDist_mm < 18) hitDist = 6;
01262                                                 if(hitDist_mm < 17) hitDist = 5;
01263                                                 if(hitDist_mm < 14) hitDist = 4;
01264                                                 if(hitDist_mm < 12) hitDist = 3;
01265                                                 if(hitDist_mm < 9) hitDist = 2;
01266                                                 if(hitDist_mm < 6) hitDist = 1;
01267                                                 if(hitDist_mm < 1) hitDist = 0; 
01268                                                 
01269                                                 hitDist_py = hitDist_mm / 6.125;                                                
01270                                                 
01271                                                 //cout << (xbinCrTalk[1]-xbinCrTalk[0]) << "\t" << (ybinCrTalk[1]-ybinCrTalk[0]) << "\t" << hitDist << endl;
01272                                         }
01273                                         else {                                  
01274                                                 // Calculation of distance with Pythagoras      
01275                                                 hitDist_py = TMath::Sqrt(distance_x*distance_x + distance_y*distance_y); // doesn't work for PMTs 18,20,22,24
01276                                         }
01277                                         
01278                                         hHitDist[iPMT]->Fill(hitDist); //iPMT
01279                                         hHitDist_py[iPMT]->Fill(hitDist_py); //iPMT
01280                                         hHitMultperPMT_CrTalk[iPMT]->Fill(hitsPerPMT[iPMT-1]); //iPMT
01281                                         hEntries2dAll->Fill(xbinCrTalk[0],ybinCrTalk[0]);
01282                                         hEntries2dAll->Fill(xbinCrTalk[1],ybinCrTalk[1]);
01283                                 }                               
01284                         }                       
01285                         
01286                         if(hitsPerPMT[iPMT-1] < 2 || hitsPerPMT[iPMT-1] > 2) {
01287                                 hitPMTmore->Fill(iPMT);
01288                         }
01289                 }       
01290         }
01291         // ############################################ END CROSSTALK ANALYSIS   ##################################################
01292    
01293    
01294             
01295          
01296          
01297          
01298 
01299    // process TRB lightpulser events
01300    if (fTrbInputEvent && fTrbInputEvent->IsValid()) {
01301      for (UInt_t RICHTRB = fPar->RICHFirstTRB; RICHTRB <= fPar->RICHLastTRB; RICHTRB++) {
01302          // Loop over all TRBs which deliver RICH data
01303        if (RICHTRB > TRB_TDC3_NUMBOARDS)
01304          continue;
01305        
01306        TTrbData* theTRB = dynamic_cast<TTrbData*>(fTrbInputEvent->getEventElement(RICHTRB));
01307        if (theTRB == 0)
01308          continue; //no data from this TRB available;
01309        
01310        for (unsigned int i = 0; i<theTRB->fHits.size(); i++) {
01311          UShort_t tdc = theTRB->fHits[i].GetTDC();
01312          UShort_t ch = theTRB->fHits[i].GetChannel();
01313          UChar_t edge = theTRB->fHits[i].GetEdgeType();
01314          Double_t dtm = theTRB->fHits[i].GetDeltaT();
01315          Double_t xbin=-1;
01316          Double_t ybin=-1;
01317          Int_t ncell=-1;
01318 
01319          if (ch==0) continue;  // skip tdc entries from reference channel 0 !
01320          if (fTRBindex[RICHTRB*64+ch] != 0) {
01321            xbin=fTRBindex[RICHTRB*64+ch]->xbin;
01322            ybin=fTRBindex[RICHTRB*64+ch]->ybin;
01323            ncell=fTRBindex[RICHTRB*64+ch]->ncell;
01324            Int_t mapmt = fTRBindex[RICHTRB*64+ch]->mapmt;
01325 
01326            hLedPulserEntries2d->Fill(xbin, ybin);
01327            hLedPulserEntriesPerChannel->Fill(ncell);
01328            hitsPerPMT[mapmt-1]++;
01329          }
01330        }
01331 
01332      } // trb
01333    }
01335    
01336    for (Int_t mapmt = 1; mapmt<=24; mapmt++)
01337      hLedPulserHitVsMapmt->Fill(mapmt-1, hitsPerPMT[mapmt-1]);
01338    
01339 }
01340 
01341 
01342 
01343 void TRICHProc::DoAnalysis(
01344                 const vector<CbmRichHitLight>& hits,
01345                 bool isElectron,
01346                 bool isPion,
01347                 bool isMuon)
01348 {
01349         //cout << "hits.size() = " << hits.size() << endl;
01350         if (hits.size() < 1) return;
01351 
01352     //[0] is electron, [1] is pion, [2] is muon, [3] is all
01353         Int_t hIndex = -1;
01354         if (isElectron) hIndex = 0;
01355         else if (isPion) hIndex = 1;
01356         else if (isMuon) hIndex = 2;
01357         if (hIndex < 0) hIndex = 2;
01358 
01359         hNofHitsEvAll[hIndex]->Fill(hits.size());
01360         hNofHitsEvAll[3]->Fill(hits.size());
01361 
01362         fRingFinder->DoFind(hits);
01363         vector<CbmRichRingLight*> rings = fRingFinder->GetFoundRings();
01364         hNofFoundRings->Fill(rings.size());
01365         if (rings.size() <= 0) return;
01366         CbmRichRingLight* ring = rings[0];
01367         //if (rings.size() == 2) cout << "2 rings were found" << endl;
01368 
01369         Int_t nofHits = ring->GetNofHits();
01370         for (Int_t iH = 0; iH < hits.size(); iH++){
01371                 Bool_t hitFound = false;
01372                 for (Int_t i = 0; i < nofHits; i++) {
01373                    if (ring != NULL && ring->GetHit(i).fX == hits[iH].fX && ring->GetHit(i).fY == hits[iH].fY){
01374                            hitFound = true;
01375                    }
01376                 }
01377                 if (hitFound == false) hNotAssignedHitsXY->Fill(hits[iH].fX, hits[iH].fY);
01378         }
01379 
01380         hNofHitsRingFinder[hIndex]->Fill(ring->GetNofHits());
01381         hNofHitsRingFinder[3]->Fill(ring->GetNofHits());
01382         //hNofHitsRingFinder_cor->Fill(1.0*ring->GetNofHits()/fPress0*1000.0*(fTemp0+273.15)/273.15); //using temp sensor of camera
01383         hNofHitsRingFinder_cor->Fill(1.0*ring->GetNofHits()/fPTB*1000.0*fTT1/273.15); // using temp-sensor of gas-system
01384         //hNofHitsRingFinder_cor2->Fill(1.0*ring->GetNofHits()/fPress0*1000.0/fPress0*1000.0*fTT1/273.15*fTT1/273.15); // using temp-sensor of gas-system
01385         //hNofHitsRingFinder_cor3->Fill(1.0*ring->GetNofHits()/(1-fTT1*fTT1/(fPress0*fPress0))*(1-273.15*273.15/1000000)); // using temp-sensor of gas-system
01386 
01387         hNofHitsEvRingFinder[hIndex]->Fill(hits.size());
01388         hNofHitsEvRingFinder[3]->Fill(hits.size());
01389         hNofHitsEvRingFinderEff[hIndex]->Divide(hNofHitsRingFinder[hIndex], hNofHitsEvAll[hIndex]);
01390         hNofHitsEvRingFinderEff[hIndex]->Divide(hNofHitsRingFinder[3], hNofHitsEvAll[3]);
01391         //cout << "Ring finder eff:" << 100.* (double)hNofHitsRingFinder[3]->GetEntries()/ hNofHitsAll[3]->GetEntries() << endl;
01392 
01393         //fit COP
01394         for (int iR = 0; iR < rings.size(); iR++){
01395                 fRingFitter->DoFit(rings[iR]);
01396         }
01397         FillCircleFitHistograms(hIndex, ring, hits.size());// fill histograms for a specific particle
01398         FillCircleFitHistograms(3, ring, hits.size());// fill histograms for all particle. [3] is all
01399         DrawEventHits(hits, fNrSED);
01400         DrawRingHits(ring, fNrSED);
01401         DrawCircles(rings, fNrSED);
01402 
01403         for (int iR = 0; iR < rings.size(); iR++){
01404                 fEllipseFitter->DoFit(rings[iR]);
01405         }
01406         if (hIndex>=0) FillEllipseFitHistograms(hIndex, ring, hits.size());// fill histograms for a specific particle
01407         FillEllipseFitHistograms(3, ring, hits.size());// fill histograms for a specific particle
01408         DrawEllipses(rings, fNrSED);
01409 
01410         // rotate single event printing to next histogram ONLY if last one had a non-empty event
01411         if ( rings.size() == 2 ){// && ring->GetBaxis()/ring->GetAaxis() < 0.7) {
01412                 fNrSED++;
01413                 if (fNrSED >= hSED.size()) fNrSED = 0;
01414         }
01415         //if (ring != NULL) delete ring;
01416 }
01417 
01418 void TRICHProc::FillCircleFitHistograms(
01419                 Int_t hIndex,
01420                 CbmRichRingLight* ring,
01421                 int nofEventHits)
01422 {
01423    if (hIndex != -1){
01424       if (ring->GetRadius() > 0. && ring->GetCenterX() > 0. && ring->GetCenterY() > 0.) {
01425          hCircleFitRadius[hIndex]->Fill(ring->GetRadius());
01426          hCircleFitRadius_cor[hIndex]->Fill(ring->GetRadius()/(TMath::Pi()/2-fTT1/fPTB)*(TMath::Pi()/2-273.15/1000));
01427          hCircleFitCenter[hIndex]->Fill(ring->GetCenterX(), ring->GetCenterY());
01428          hCircleFitChi2[hIndex]->Fill(ring->GetChi2() / ring->GetNofHits());
01429          for (int iH = 0; iH < ring->GetNofHits(); iH++){
01430                  double xc = ring->GetCenterX();
01431                  double yc = ring->GetCenterY();
01432                  CbmRichHitLight h = ring->GetHit(iH);
01433                  double dr = sqrt((xc - h.fX)*(xc - h.fX) + (yc - h.fY)*(yc - h.fY)) - ring->GetRadius();
01434                  hCircleFitDR[hIndex]->Fill(dr);
01435          }
01436       }
01437       // good fitted ring
01438       if (ring->GetRadius() > 20 && ring->GetRadius() < 70) hNofHitsEvCircleFit[hIndex]->Fill(nofEventHits);
01439       hNofHitsEvCircleFitEff[hIndex]->Divide(hNofHitsEvCircleFit[hIndex], hNofHitsEvAll[hIndex]);
01440    }
01441 }
01442 
01443 void TRICHProc::FillEllipseFitHistograms(
01444                 Int_t hIndex,
01445                 CbmRichRingLight* ring,
01446                 int nofEventHits)
01447 {
01448    if (ring->GetAaxis() > 0. && ring->GetBaxis() > 0. && ring->GetPhi() != -1. && ring->GetCenterX() != -1. && ring->GetCenterY()!=-1.) {
01449       hEllipseFitAaxis[hIndex]->Fill(ring->GetAaxis());
01450       hEllipseFitBaxis[hIndex]->Fill(ring->GetBaxis()); 
01451       hEllipseFitBoverA[hIndex]->Fill(ring->GetBaxis()/ring->GetAaxis());  
01452       hEllipseFitPhi[hIndex]->Fill(ring->GetPhi());    
01453       hEllipseFitCenter[hIndex]->Fill(ring->GetCenterX(),ring->GetCenterY());
01454       hEllipseFitChi2[hIndex]->Fill(ring->GetChi2() / ring->GetNofHits());
01455    }
01456    //good fitted ring with ellipse fitter
01457    if (ring->GetAaxis() > 20 && ring->GetAaxis() < 80 && ring->GetBaxis() > 20
01458                    && ring->GetBaxis() < 80) hNofHitsEvEllipseFit[hIndex]->Fill(nofEventHits);
01459    hNofHitsEvEllipseFitEff[hIndex]->Divide(hNofHitsEvEllipseFit[hIndex], hNofHitsEvAll[hIndex]);
01460 }
01461 
01462 void TRICHProc::DrawRingHits(
01463                 CbmRichRingLight* ring,
01464                 Int_t sedNum)
01465 {
01466         if (ring->GetNofHits() > fMaxNofHitsInEventDraw) return;
01467 
01468    Int_t nofHits = ring->GetNofHits();
01469    for (Int_t i = 0; i < nofHits; i++) {
01470           fRingHits[sedNum][i]->SetX1(ring->GetHit(i).fX);
01471           fRingHits[sedNum][i]->SetY1(ring->GetHit(i).fY);
01472           fRingHits[sedNum][i]->SetR1(3);
01473           fRingHits[sedNum][i]->SetR2(3);
01474    }
01475    for (Int_t i = nofHits; i < fMaxNofHitsInEventDraw; i++){
01476            fRingHits[sedNum][i]->SetX1(-999999);
01477    }
01478 }
01479 
01480 void TRICHProc::DrawEventHits(
01481                 const vector<CbmRichHitLight>& hits,
01482                 Int_t sedNum)
01483 {
01484    if (hits.size() > fMaxNofHitsInEventDraw) return;
01485    Int_t nofHits = hits.size();
01486    for (Int_t i = 0; i < nofHits; i++) {
01487            fEventHits[sedNum][i]->SetX1(hits[i].fX);
01488            fEventHits[sedNum][i]->SetY1(hits[i].fY);
01489            fEventHits[sedNum][i]->SetR1(4);
01490            fEventHits[sedNum][i]->SetR2(4);
01491    }
01492    for (Int_t i = nofHits; i < fMaxNofHitsInEventDraw; i++){
01493            fEventHits[sedNum][i]->SetX1(-999999);
01494    }
01495 }
01496 
01497 void TRICHProc::DrawCircles(
01498                 const vector<CbmRichRingLight*>& rings,
01499                 Int_t sedNum)
01500 {
01501         if (rings.size() <= 0) return;
01502         for (int iR = 0; iR < rings.size(); iR++){
01503            fCircle[sedNum][iR]->SetX1(rings[iR]->GetCenterX());
01504            fCircle[sedNum][iR]->SetY1(rings[iR]->GetCenterY());
01505            fCircle[sedNum][iR]->SetR1(rings[iR]->GetRadius());
01506            fCircle[sedNum][iR]->SetR2(rings[iR]->GetRadius());
01507    }
01508    stringstream ss;
01509    ss.precision(3);
01510    ss << "(x,y,r,n):(" << rings[0]->GetCenterX() << ", " << rings[0]->GetCenterY() << ","
01511                    << rings[0]->GetRadius()<< ", " << rings[0]->GetNofHits() << ")";
01512    fCircleText[sedNum]->SetText(5, 200, ss.str().c_str());
01513 }
01514 
01515 
01516 void TRICHProc::DrawEllipses(
01517                 const vector<CbmRichRingLight*>& rings,
01518                 Int_t sedNum)
01519 {
01520         for (int iR = 0; iR < rings.size(); iR++){
01521            fEllipse[sedNum][iR]->SetX1(rings[iR]->GetCenterX());
01522            fEllipse[sedNum][iR]->SetY1(rings[iR]->GetCenterY());
01523            fEllipse[sedNum][iR]->SetR1(rings[iR]->GetAaxis());
01524            fEllipse[sedNum][iR]->SetR2(rings[iR]->GetBaxis());
01525            fEllipse[sedNum][iR]->SetTheta(rings[iR]->GetPhi() * 180. / 3.1415);
01526         }
01527    stringstream ss;
01528    ss.precision(3);
01529    ss << "(x,y,a,b,p,n):(" << rings[0]->GetCenterX() << ", " << rings[0]->GetCenterY() << ","
01530       << rings[0]->GetAaxis()<< ", "<< rings[0]->GetBaxis()<< ", " << rings[0]->GetPhi()<< ", "
01531       << rings[0]->GetNofHits() << ")";
01532    fEllipseText[sedNum]->SetText(5, 185, ss.str().c_str());
01533 }