Changeset 4429 for trunk/MagicSoft/Mars/mtemp/mifae/macros
- Timestamp:
- 07/28/04 11:01:15 (20 years ago)
- Location:
- trunk/MagicSoft/Mars/mtemp/mifae/macros
- Files:
-
- 7 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mtemp/mifae/macros/alpha_plot.C
r3947 r4429 25 25 } 26 26 27 //void alpha_plot(TString f_on_name = "../HillasFiles/20040319_Mrk421_30_15.KDummy*.root", 28 void alpha_plot(TString f_on_name = "../HillasFiles/20040319_OffMrk421_30_15.KDummy*.root", 29 TString f_off_name = "../HillasFiles/20040319_OffMrk421_30_15.KDummy*.root", 27 void alpha_plot(TString f_on_name = "../HillasFiles/Mrk421/*_H.root", 28 TString f_off_name = "../HillasFiles/OffMrk421/*_H.root", 30 29 TString f_src_name = "../HillasFiles/20040319_OffMrk421.fake.pos") 31 /*void alpha_plot(TString f_on_name = "../HillasFiles/mrk421OnMisp3015*.root",32 TString f_off_name = "../HillasFiles/mrk421OffMisp3015*.root",33 TString f_src_name = "../HillasFiles/20040215_Mrk421.pos")34 */35 30 { 36 31 … … 38 33 39 34 //cuts 40 Float_t sizemin = 500.; //[ADC]35 Float_t sizemin = 2000.; //[ADC] 41 36 Float_t sizemax = 10000000000.; //[ADC] 42 Float_t widthmin = 0. ;43 Float_t widthmax = 0. 8;44 Float_t lengthmin = 0. ;45 Float_t lengthmax = 0. 8;46 Float_t distmin = 0. ;47 Float_t distmax = 2.;37 Float_t widthmin = 0.06; 38 Float_t widthmax = 0.12; 39 Float_t lengthmin = 0.10; 40 Float_t lengthmax = 0.26; 41 Float_t distmin = 0.3; 42 Float_t distmax = 1.2; 48 43 Float_t alphamin = 0.; 49 44 Float_t alphamax = 90.; … … 52 47 Float_t sigexccmin = 0.; 53 48 Float_t sigexccmax = 30.; 54 Float_t bkgnormmin = 0.;49 Float_t bkgnormmin = 30.; 55 50 Float_t bkgnormmax = 90.; 56 51 … … 173 168 174 169 //cuts 175 TString sizestr = " MHillas.fSize < ";170 TString sizestr = "(MHillas.fSize < "; 176 171 sizestr += sizemin; 177 sizestr += " ||";172 sizestr += ") || ("; 178 173 sizestr += "MHillas.fSize > "; 179 174 sizestr += sizemax; 175 sizestr += ")"; 180 176 MF sizefilter(sizestr); 181 177 182 TString widthstr = " {MHillas.fWidth/315.} < ";178 TString widthstr = "({MHillas.fWidth/315.} < "; 183 179 widthstr += widthmin; 184 widthstr += " ||";180 widthstr += ") || ("; 185 181 widthstr += "{MHillas.fWidth/315.} > "; 186 182 widthstr += widthmax; 183 widthstr += ")"; 187 184 MF widthfilter(widthstr); 188 185 189 TString lengthstr = " {MHillas.fLength/315.} < ";186 TString lengthstr = "({MHillas.fLength/315.} < "; 190 187 lengthstr += lengthmin; 191 lengthstr += " ||";188 lengthstr += ") || ("; 192 189 lengthstr += "{MHillas.fLength/315.} > "; 193 190 lengthstr += lengthmax; 191 lengthstr += ")"; 194 192 MF lengthfilter(lengthstr); 195 193 196 TString diststr = " {MHillasSrc.fDist/315.} < ";194 TString diststr = "({MHillasSrc.fDist/315.} < "; 197 195 diststr += distmin; 198 diststr += " ||";196 diststr += ") || ("; 199 197 diststr += "{MHillasSrc.fDist/315.} > "; 200 198 diststr += distmax; 199 diststr += ")"; 201 200 MF distfilter(diststr); 202 201 203 TString alphastr = " {abs(MHillasSrc.fAlpha)} < ";202 TString alphastr = "({abs(MHillasSrc.fAlpha)} < "; 204 203 alphastr += alphamin; 205 alphastr += " ||";204 alphastr += ") || ("; 206 205 alphastr += "{abs(MHillasSrc.fAlpha)} > "; 207 206 alphastr += alphamax; 207 alphastr += ")"; 208 208 MF alphafilter(alphastr); 209 209 … … 219 219 MContinue cont_odd(&oddfilter); 220 220 221 MSrcPosFromFile srccalc(f_src_name); 221 // MSrcPosFromFile srccalc(f_src_name); 222 MSrcPlace srccalc; 222 223 223 224 MHillasSrcCalc csrc_on; … … 245 246 tlist_on.AddToList(&csrc_on); 246 247 tlist_on.AddToList(&fsrcpos_on); 247 tlist_on.AddToList(&cont_odd);248 // tlist_on.AddToList(&cont_odd); 248 249 tlist_on.AddToList(&cont_size); 249 250 tlist_on.AddToList(&cont_width); … … 357 358 read_off.DisableAutoScheme(); 358 359 359 srccalc.SetMode(MSrcP osFromFile::kOff);360 srccalc.SetMode(MSrcPlace::kOff); 360 361 361 362 MHillasSrcCalc csrc_off; … … 376 377 tlist_off.AddToList(&csrc_off); 377 378 tlist_off.AddToList(&fsrcpos_off); 378 tlist_off.AddToList(&cont_even);379 // tlist_off.AddToList(&cont_even); 379 380 tlist_off.AddToList(&cont_size); 380 381 tlist_off.AddToList(&cont_width); -
trunk/MagicSoft/Mars/mtemp/mifae/macros/findstars.C
r4294 r4429 47 47 48 48 49 void findstars(const TString filename=" dc_2004_03_17_01_16_51_20440_Mrk421.root", const TString directory="/nfs/magic/CaCodata/rootdata/Mrk421/Period015/2004_03_17/", const UInt_t numEvents = 0)49 void findstars(const TString filename="20040422_*_D_Mrk421_E.root", const TString directory="/local_disk/Data/rootdata/Mrk421/Period016/2004_04_22/", const UInt_t numEvents = 0) 50 50 { 51 51 … … 81 81 MGeomApply geomapl; 82 82 TString continuoslightfile = 83 // "/home/Javi/mnt_magic_data/CaCo/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root"; 84 "/nfs/magic/CaCodata/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root"; 83 "/local_disk/CaCoData/rootdata/Miscellaneous/Period016/2004_04_16/dc_2004_04_16_04_46_18_22368_Off3c279-2CL100.root"; 85 84 86 85 Float_t mindc = 0.9; //[uA] … … 93 92 const TArrayS blindpixels(numblind,(Short_t*)x); 94 93 Float_t ringinterest = 100; //[mm] 95 Float_t tailcut = 2.5;96 UInt_t integratedevents = 1;94 Float_t tailcut = 3.5; 95 UInt_t integratedevents = 4; 97 96 98 97 MFindStars findstars; -
trunk/MagicSoft/Mars/mtemp/mifae/macros/hillasONOFF.C
r4035 r4429 27 27 28 28 // define file, tree, canvas 29 TFile *fdata = new TFile(onname); 29 TChain *fdata = new TChain; 30 fdata->AddFile(onname); 30 31 TTree *t = Parameters; 31 32 32 33 // define file, tree, canvas 33 TFile *fdata2 = new TFile(offname); 34 TChain *fdata2 = new TChain; 35 fdata2->AddFile(offname); 34 36 TTree *t2 = Parameters; 35 37 -
trunk/MagicSoft/Mars/mtemp/mifae/macros/hvnotnominal.C
r3977 r4429 47 47 48 48 49 void hvnotnominal(const TString filename="20040319_20821_D_Mrk421_S.root", const TString directory="/nfs/magic/CaCodata/2004_03_19/" )49 void hvnotnominal(const TString filename="20040319_20821_D_Mrk421_S.root", const TString directory="/nfs/magic/CaCodata/2004_03_19/", Float_t percent = 0.01) 50 50 { 51 51 … … 86 86 MGeomApply geomapl; 87 87 MFHVNotNominal fHVNominal; 88 TString hvnominal = "/ nfs/ifae.es/user/j/jlopez/Mars/ConfigFiles/HVSettings_FF35q.conf";88 TString hvnominal = "/mnt/users/jlopez/Mars/Files4Mars/Config/HVSettings_FF35q.conf"; 89 89 fHVNominal.SetHVNominalValues(hvnominal); 90 90 // fHVNominal.SetMaxNumPixelsDeviated(10); … … 105 105 // 106 106 107 if (!evtloop.Eventloop()) 107 MHCamera hvnotnominal(geomcam); 108 MHCamera lowhvnotnominal(geomcam); 109 MHCamera uphvnotnominal(geomcam); 110 111 if (!evtloop.PreProcess()) 108 112 return; 109 113 114 TArrayD nominal = fHVNominal.GetHVNominal(); 115 TArrayD lownominal = nominal; 116 TArrayD upnominal = nominal; 117 118 for (UInt_t pix=0; pix<nominal.GetSize(); pix++) 119 { 120 lownominal[pix] *= (1-percent); 121 upnominal[pix] *= (1+percent); 122 } 123 124 125 while(tlist.Process()) 126 { 127 128 hvnotnominal.CntCamContent(hvcam,lownominal,0,kFALSE); 129 hvnotnominal.CntCamContent(hvcam,upnominal,0,kTRUE); 130 131 lowhvnotnominal.CntCamContent(hvcam,lownominal,0,kFALSE); 132 uphvnotnominal.CntCamContent(hvcam,upnominal,0,kTRUE); 133 } 134 135 136 evtloop.PostProcess(); 137 110 138 tlist.PrintStatistics(); 139 140 TCanvas c1; 141 c1.Divide(2,1); 142 c1.cd(1); 143 lowhvnotnominal.SetPrettyPalette(); 144 lowhvnotnominal.Draw(); 145 gPad->cd(1); 146 gPad->Modified(); 147 gPad->Update(); 148 c1.cd(2); 149 uphvnotnominal.SetPrettyPalette(); 150 uphvnotnominal.Draw(); 151 gPad->cd(2); 152 gPad->Modified(); 153 gPad->Update(); 154 155 TCanvas c2; 156 hvnotnominal.SetPrettyPalette(); 157 hvnotnominal.Draw(); 158 gPad->cd(1); 159 gPad->Modified(); 160 gPad->Update(); 161 162 if (!HandleInput()) 163 break; 111 164 112 165 } -
trunk/MagicSoft/Mars/mtemp/mifae/macros/lightcurve.C
r4408 r4429 1 Double_t ChiSquareNDof(TH1D *h1, TH1D *h2) 2 { 3 Double_t chiq = 0.; 4 Double_t chi; 5 Double_t error; 6 Int_t nbinsnozero = 0; 7 8 Int_t nbins = h1->GetNbinsX(); 9 if (nbins != h2->GetNbinsX() || nbins == 0) 10 return -1; 11 12 for (UInt_t bin=1; bin<=nbins; bin++) 13 { 14 error = sqrt(h1->GetBinError(bin)*h1->GetBinError(bin) + 15 h2->GetBinError(bin)*h2->GetBinError(bin)); 16 if (error != 0) 17 { 18 chi = (h1->GetBinContent(bin)-h2->GetBinContent(bin))/error; 19 chiq += chi*chi; 20 nbinsnozero++; 21 } 22 } 23 24 return (nbinsnozero>0?chiq/nbinsnozero:0); 25 } 26 27 Int_t GetBin(Double_t size, Int_t numberSizeBins, Double_t sizeBins[]) 28 { 29 30 Int_t result = -1; 31 32 Int_t lowerbin = 0; 33 Int_t upperbin = numberSizeBins; 34 Int_t bin; 35 36 Int_t count = 0; 37 38 if (size >= sizeBins[0]) 39 { 40 while (upperbin - lowerbin > 1 && count++<=numberSizeBins) 41 { 42 bin = (upperbin+lowerbin)/2; 43 if (size >= sizeBins[bin]) 44 lowerbin = bin; 45 else 46 upperbin = bin; 47 } 48 result = count<=numberSizeBins?lowerbin:-1; 49 } 50 51 return result; 52 53 } 54 1 2 Double_t ChiSquareNDof(TH1D *h1, TH1D *h2); 3 Int_t GetBin(Double_t size, Int_t numberSizeBins, Double_t sizeBins[]); 55 4 56 5 void lightcurve(TString f_on_name = "../HillasFiles/Mrk421/*_H.root", … … 67 16 // Configuration 68 17 const Bool_t debug = kFALSE; 69 const Double_t timebin = 600.; //[sec]70 TString psname = "./2004042 0_Mrk421.600s.ps";18 const Double_t timebin = 1380.; //[sec] 19 TString psname = "./20040421_Mrk421.1380s.ps"; 71 20 72 21 //Constanst … … 94 43 TObjArray timediffHistoOn; 95 44 96 const Int_t numberSizeBins = 4;97 Double_t sizeBins[numberSizeBins] = {1 292., 1897., 2785., 4087.}; //[Photons]45 const Int_t numberSizeBins = 3; 46 Double_t sizeBins[numberSizeBins] = {1897., 2785., 4087.}; //[Photons] 98 47 99 48 //cuts 100 49 101 Double_t widthmin[numberSizeBins] = { 0.06, 0.06, 0.06, 0.06 }; 102 Double_t widthmax[numberSizeBins] = { 0.10, 0.10, 0.12, 0.12 }; 103 Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10, 0.10 }; 104 Double_t lengthmax[numberSizeBins] = { 0.20, 0.25, 0.26, 0.36 }; 105 Double_t distmin[numberSizeBins] = { 0.30, 0.30, 0.30, 0.30 }; 106 Double_t distmax[numberSizeBins] = { 1.20, 1.20, 1.20, 1.20 }; 50 Double_t widthmin[numberSizeBins] = { 0.06, 0.06, 0.06 }; 51 Double_t widthmax[numberSizeBins] = { 0.10, 0.12, 0.12 }; 52 Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10 }; 53 Double_t lengthmax[numberSizeBins] = { 0.25, 0.26, 0.36 }; 54 Double_t distmin[numberSizeBins] = { 0.30, 0.30, 0.30 }; 55 Double_t distmax[numberSizeBins] = { 1.20, 1.20, 1.20 }; 56 57 // const Int_t numberSizeBins = 4; 58 // Double_t sizeBins[numberSizeBins] = {1292., 1897., 2785., 4087.}; //[Photons] 59 60 // //cuts 61 62 // Double_t widthmin[numberSizeBins] = { 0.06, 0.06, 0.06, 0.06 }; 63 // Double_t widthmax[numberSizeBins] = { 0.10, 0.10, 0.12, 0.12 }; 64 // Double_t lengthmin[numberSizeBins] = { 0.10, 0.10, 0.10, 0.10 }; 65 // Double_t lengthmax[numberSizeBins] = { 0.20, 0.25, 0.26, 0.36 }; 66 // Double_t distmin[numberSizeBins] = { 0.30, 0.30, 0.30, 0.30 }; 67 // Double_t distmax[numberSizeBins] = { 1.20, 1.20, 1.20, 1.20 }; 107 68 108 69 //alpha plot integration for gammas 109 70 Double_t sigexccmin = 0.; 110 Double_t sigexccmax = 20.;71 Double_t sigexccmax = 15.; 111 72 Double_t bkgnormmin = 40.; 112 73 Double_t bkgnormmax = 80.; … … 128 89 MRawEvtHeader evtheader_on; 129 90 MTime time_on; 91 MLiveTime livetime_on; 130 92 MHillas hillas_on; 131 93 MHillasSrc hillassrc_on; … … 137 99 plist_on.AddToList(&evtheader_on); 138 100 plist_on.AddToList(&time_on); 101 plist_on.AddToList(&livetime_on); 139 102 plist_on.AddToList(&hillas_on); 140 103 plist_on.AddToList(&hillassrc_on); … … 186 149 MSrcPlace srcplace; 187 150 MHillasSrcCalc csrc_on; 188 151 152 MLiveTimeCalc livetimecalc_on; 153 livetimecalc_on.SetRealTimeBinSize(timebin); 154 189 155 // // prints 190 156 // MPrint pevent("MRawEvtHeader"); … … 195 161 //tasklist 196 162 tlist_on.AddToList(&read_on); 163 tlist_on.AddToList(&livetimecalc_on); 197 164 tlist_on.AddToList(&srcplace); 198 165 tlist_on.AddToList(&csrc_on); … … 212 179 return; 213 180 214 Double_t mjdFirstEventofBin=0;215 Double_t mjdFirstEvent=0;216 217 MTime mtimeFirstEventofBin;218 MTime mtimeFirstEvent;219 220 Double_t mjdLastEvent=0;221 MTime mtimeLastEvent=0;222 UInt_t evtLastEvent;223 UInt_t runLastEvent=0;224 MTime mtimeLastEvent;225 226 Bool_t flag = kFALSE;227 228 181 numberOnEventsAfterCleaning[numberTimeBins-1] = 0; 229 182 zenithMinimumOn[numberTimeBins-1] = 100.; … … 252 205 TF1 *f1 = new TF1("exp","expo",mintimediff,maxtimediff); 253 206 207 Double_t lastSrcPosX = 0; 208 Double_t lastSrcPosY = 0; 209 254 210 while(tlist_on.Process()) 255 211 { 256 // Compute live time 257 UInt_t run = runheader_on.GetRunNumber(); 258 UInt_t evt = evtheader_on.GetDAQEvtNumber(); 259 Double_t mjd = time_on.GetMjd(); 212 numberOnEventsAfterCleaning[numberTimeBins-1]++; 260 213 261 numberOnEventsAfterCleaning[numberTimeBins-1]++; 262 263 if (mjd == 0) 264 { 265 266 if (debug) 267 { 268 cout << "MTime::GetTime() == 0" << endl; 269 cout.precision(15); 270 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl; 271 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl; 272 mtimeLastEvent.Print(); 273 time_on.Print(); 274 } 275 276 flag = kTRUE; 277 } 278 else if (mjd-mjdLastEvent < 0 && flag) 279 { 280 281 if (debug) 282 { 283 cout << "mjd-mjdLastEvent < 0" << endl; 284 cout.precision(15); 285 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl; 286 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl; 287 mtimeLastEvent.Print(); 288 time_on.Print(); 289 } 290 291 flag = kTRUE; 292 } 293 else if (mjd-mjdLastEvent == 0 && flag) 294 { 295 296 if (debug) 297 { 298 cout << "mjd-mjdLastEvent == 0" << endl; 299 cout.precision(15); 300 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl; 301 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl; 302 mtimeLastEvent.Print(); 303 time_on.Print(); 304 } 305 306 flag = kTRUE; 307 } 308 else if (evt-evtLastEvent <= 0) 309 { 310 311 if (debug) 312 { 313 cout << "evt-evtLastEvent <= 0" << endl; 314 cout.precision(15); 315 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl; 316 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl; 317 mtimeLastEvent.Print(); 318 time_on.Print(); 319 } 320 321 flag = kTRUE; 322 } 323 324 else if ((Int_t)mjd == mjd) 325 { 326 327 if (debug) 328 { 329 cout << "(Int_t)mjd == mjd" << endl; 330 cout.precision(15); 331 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl; 332 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl; 333 mtimeLastEvent.Print(); 334 time_on.Print(); 335 } 336 337 flag = kTRUE; 338 } 214 if (srcpos_on.GetX() == 0. && srcpos_on.GetY() == 0.) 215 srcpos_on.SetXY(lastSrcPosX,lastSrcPosY); 339 216 else 340 217 { 341 342 if (flag && debug) 218 lastSrcPosX = srcpos_on.GetX(); 219 lastSrcPosY = srcpos_on.GetY(); 220 } 221 srcplace.CallProcess(); 222 csrc_on.CallProcess(); 223 224 225 Double_t size = hillas_on.GetSize(); 226 Double_t width = hillas_on.GetWidth()*geomcam.GetConvMm2Deg(); 227 Double_t length = hillas_on.GetLength()*geomcam.GetConvMm2Deg(); 228 Double_t meanx = hillas_on.GetMeanX()*geomcam.GetConvMm2Deg(); 229 Double_t meany = hillas_on.GetMeanY()*geomcam.GetConvMm2Deg(); 230 Double_t centerdist = TMath::Sqrt(meanx*meanx + meany*meany); 231 Double_t dist = hillassrc_on.GetDist()*geomcam.GetConvMm2Deg(); 232 Double_t alpha = hillassrc_on.GetAlpha(); 233 Double_t srcposx = srcpos_on.GetX(); 234 Double_t srcposy = srcpos_on.GetY(); 235 Double_t zenith = drive_on.GetNominalZd(); 236 237 238 Int_t sizebin = GetBin(size,numberSizeBins,sizeBins); 239 240 if (sizebin >= 0) 241 { 242 if (width > widthmin[sizebin] && width < widthmax[sizebin]) 343 243 { 344 345 cout << "flag = TRUE" << endl; 346 cout.precision(15); 347 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl; 348 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl; 349 mtimeLastEvent.Print(); 350 time_on.Print(); 351 352 flag = kFALSE; 353 354 } 355 356 if ( run != runLastEvent ) 357 { 358 if ( mjdLastEvent != 0 ) 244 if (length > lengthmin[sizebin] && length < lengthmax[sizebin]) 359 245 { 360 timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay; 361 362 cout.precision(10); 363 cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl; 364 if (flag && debug) 365 { 366 cout << "last event:\t run "<< runLastEvent << " event " << evtLastEvent << " mjd " << mjdLastEvent << endl; 367 cout << "present event:\t run " << run << " event " << evt << " mjd " << mjd << endl; 368 mtimeLastEvent.Print(); 369 time_on.Print(); 370 } 371 372 } 373 374 mjdFirstEvent = mjd; 375 } 376 else if (mjdFirstEventofBin != 0) 377 hTimeDiff_on_timebin->Fill((mjd-mjdLastEvent)*kDay*kSec); 378 379 if (mjdFirstEventofBin == 0) 380 { 381 mjdFirstEvent = mjd; 382 mjdFirstEventofBin = mjd; 383 } 384 385 evtLastEvent = evt; 386 runLastEvent = run; 387 mjdLastEvent = mjd; 388 mtimeLastEvent = time_on; 389 390 Double_t size = hillas_on.GetSize(); 391 Double_t width = hillas_on.GetWidth()*geomcam.GetConvMm2Deg(); 392 Double_t length = hillas_on.GetLength()*geomcam.GetConvMm2Deg(); 393 Double_t meanx = hillas_on.GetMeanX()*geomcam.GetConvMm2Deg(); 394 Double_t meany = hillas_on.GetMeanY()*geomcam.GetConvMm2Deg(); 395 Double_t centerdist = TMath::Sqrt(meanx*meanx + meany*meany); 396 Double_t dist = hillassrc_on.GetDist()*geomcam.GetConvMm2Deg(); 397 Double_t alpha = hillassrc_on.GetAlpha(); 398 Double_t srcposx = srcpos_on.GetX(); 399 Double_t srcposy = srcpos_on.GetY(); 400 Double_t zenith = drive_on.GetNominalZd(); 401 402 403 Int_t sizebin = GetBin(size,numberSizeBins,sizeBins); 404 405 if (sizebin >= 0) 406 { 407 if (width > widthmin[sizebin] && width < widthmax[sizebin]) 408 { 409 if (length > lengthmin[sizebin] && length < lengthmax[sizebin]) 410 { 411 if (dist > distmin[sizebin] && centerdist < distmax[sizebin]) 412 { 413 414 //General histos 415 hSize_on->Fill(log10(size)); 416 hWidth_on->Fill(width); 417 hLength_on->Fill(length); 418 hDist_on->Fill(dist); 419 hAlpha_on_abs->Fill(TMath::Abs(alpha)); 420 hAlpha_on->Fill(alpha); 421 hSrcPos_on->Fill(srcposx,srcposy); 422 423 // Time bin histos 424 hAlpha_on_abs_timebin->Fill(TMath::Abs(alpha)); 425 hSrcPos_on_timebin->Fill(srcposx,srcposy); 426 hCosZenith_on_timebin->Fill(TMath::Cos(zenith*kConvDegToRad)); 427 428 if (zenith != 0 && zenith < zenithMinimumOn[numberTimeBins-1]) 429 zenithMinimumOn[numberTimeBins-1] = zenith; 430 if (zenith>zenithMaximumOn[numberTimeBins-1]) 431 zenithMaximumOn[numberTimeBins-1] = zenith; 432 433 } 246 if (dist > distmin[sizebin] && centerdist < distmax[sizebin]) 247 { 248 249 //General histos 250 hSize_on->Fill(log10(size)); 251 hWidth_on->Fill(width); 252 hLength_on->Fill(length); 253 hDist_on->Fill(dist); 254 hAlpha_on_abs->Fill(TMath::Abs(alpha)); 255 hAlpha_on->Fill(alpha); 256 hSrcPos_on->Fill(srcposx,srcposy); 257 258 // Time bin histos 259 hAlpha_on_abs_timebin->Fill(TMath::Abs(alpha)); 260 hSrcPos_on_timebin->Fill(srcposx,srcposy); 261 hCosZenith_on_timebin->Fill(TMath::Cos(zenith*kConvDegToRad)); 262 263 if (zenith != 0 && zenith < zenithMinimumOn[numberTimeBins-1]) 264 zenithMinimumOn[numberTimeBins-1] = zenith; 265 if (zenith>zenithMaximumOn[numberTimeBins-1]) 266 zenithMaximumOn[numberTimeBins-1] = zenith; 267 434 268 } 435 269 } 436 270 } 271 } 272 273 274 // Check if you are overload the maxim time bin 275 if (numberTimeBins != livetime_on.GetNumberTimeBins()) 276 { 277 timeOn[numberTimeBins-1] = livetime_on.GetLiveTimeTArray().At(numberTimeBins-1); 278 meanTimeBinOn[numberTimeBins-1] = livetime_on.GetMeanRealTimeTArray().At(numberTimeBins-1); 279 widthTimeBinOn[numberTimeBins-1] = livetime_on.GetWidthRealTimeTArray().At(numberTimeBins-1); 280 281 //Compute the number of on events 282 numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90); 283 numberBkgEventsToNormOn[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90); 437 284 438 439 // Check if you are overload the maxim time bin 440 if ((mjd-mjdFirstEventofBin)*kDay >= timebin) 441 { 442 //Compute the time on 443 timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay; 444 widthTimeBinOn[numberTimeBins-1] = (mjdLastEvent-mjdFirstEventofBin)/2; 445 meanTimeBinOn[numberTimeBins-1] = mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1]; 446 447 cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl; 448 cout << "mjd " << mjd << " mjdFirstEventofBin " << mjdFirstEventofBin << " widthTimeBinOn " << widthTimeBinOn[numberTimeBins-1] << " meanTimeBinOn " << meanTimeBinOn[numberTimeBins-1] << ' ' << mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1] << endl; 449 450 //Compute the number of on events 451 numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1); 452 numberBkgEventsToNormOn[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1); 453 //Initialize variables 454 mjdFirstEvent = mjd; 455 mjdFirstEventofBin = mjd; 456 457 hTimeDiff_on_timebin->Fit("exp","RQ0"); 458 deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries()); 459 cout << "1-> Exp(" << f1->GetParameter(0) << " + " << f1->GetParameter(1) << "*x) +- [" << f1->GetParError(0) << ' ' << f1->GetParError(1) << "]" << endl; 460 cout << "Calc entries " << (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))) << " +- Nt*(" << TMath::Sqrt( f1->GetParameter(0)*f1->GetParError(0)*f1->GetParameter(0)*f1->GetParError(0) + (f1->GetParError(1)*f1->GetParError(1))/(f1->GetParameter(1)*f1->GetParameter(1))) << ") meas entries " << hTimeDiff_on_timebin->GetEntries() << " dead fraction " << deadFractionOn[numberTimeBins-1] << endl; 461 462 alphaHistoOn.AddLast(hAlpha_on_abs_timebin); 463 srcposHistoOn.AddLast(hSrcPos_on_timebin); 464 coszenithHistoOn.AddLast(hCosZenith_on_timebin); 465 timediffHistoOn.AddLast(hTimeDiff_on_timebin); 466 467 // cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin << " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl; 468 469 //Increase the number of time bins and all needed arrays 470 numberTimeBins++; 471 472 timeOn.Set(numberTimeBins); 473 numberOnEvents.Set(numberTimeBins); 474 numberBkgEventsToNormOn.Set(numberTimeBins); 475 widthTimeBinOn.Set(numberTimeBins); 476 meanTimeBinOn.Set(numberTimeBins); 477 zenithMinimumOn.Set(numberTimeBins); 478 zenithMaximumOn.Set(numberTimeBins); 479 deadFractionOn.Set(numberTimeBins); 480 numberOnEventsAfterCleaning.Set(numberTimeBins); 481 482 timeOn[numberTimeBins-1] = 0.; 483 zenithMinimumOn[numberTimeBins-1] = 100.; 484 zenithMaximumOn[numberTimeBins-1] = 0.; 485 numberOnEventsAfterCleaning[numberTimeBins-1] = 0; 486 487 alphaTitle = Form("%s%02i","hAlphaOn",numberTimeBins-1); 488 hAlpha_on_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs); 489 490 srcposTitle = Form("%s%02i","hSrcPosOn",numberTimeBins-1); 491 hSrcPos_on_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos); 492 493 coszenithTitle = Form("%s%02i","hCosZenithOn",numberTimeBins-1); 494 hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith); 495 496 timediffTitle = Form("%s%02i","hTimeDiffOn",numberTimeBins-1); 497 hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff); 498 499 flag = kTRUE; 500 } 285 // hTimeDiff_on_timebin->Fit("exp","RQ0"); 286 // deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries()); 287 // cout << "1-> Exp(" << f1->GetParameter(0) << " + " << f1->GetParameter(1) << "*x) +- [" << f1->GetParError(0) << ' ' << f1->GetParError(1) << "]" << endl; 288 // cout << "Calc entries " << (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))) << " +- Nt*(" << TMath::Sqrt( f1->GetParameter(0)*f1->GetParError(0)*f1->GetParameter(0)*f1->GetParError(0) + (f1->GetParError(1)*f1->GetParError(1))/(f1->GetParameter(1)*f1->GetParameter(1))) << ") meas entries " << hTimeDiff_on_timebin->GetEntries() << " dead fraction " << deadFractionOn[numberTimeBins-1] << endl; 289 deadFractionOn[numberTimeBins-1] = 1.; 290 291 alphaHistoOn.AddLast(hAlpha_on_abs_timebin); 292 srcposHistoOn.AddLast(hSrcPos_on_timebin); 293 coszenithHistoOn.AddLast(hCosZenith_on_timebin); 294 timediffHistoOn.AddLast(hTimeDiff_on_timebin); 295 296 numberTimeBins = livetime_on.GetNumberTimeBins(); 297 298 timeOn.Set(numberTimeBins); 299 numberOnEvents.Set(numberTimeBins); 300 numberBkgEventsToNormOn.Set(numberTimeBins); 301 widthTimeBinOn.Set(numberTimeBins); 302 meanTimeBinOn.Set(numberTimeBins); 303 zenithMinimumOn.Set(numberTimeBins); 304 zenithMaximumOn.Set(numberTimeBins); 305 deadFractionOn.Set(numberTimeBins); 306 numberOnEventsAfterCleaning.Set(numberTimeBins); 307 308 timeOn[numberTimeBins-1] = 0.; 309 zenithMinimumOn[numberTimeBins-1] = 100.; 310 zenithMaximumOn[numberTimeBins-1] = 0.; 311 numberOnEventsAfterCleaning[numberTimeBins-1] = 0; 312 313 alphaTitle = Form("%s%02i","hAlphaOn",numberTimeBins-1); 314 hAlpha_on_abs_timebin = new TH1F (alphaTitle,"",nbins_abs,minalpha_abs,maxalpha_abs); 315 316 srcposTitle = Form("%s%02i","hSrcPosOn",numberTimeBins-1); 317 hSrcPos_on_timebin = new TH2F (srcposTitle,"",nbins_srcpos,minsrcpos,maxsrcpos,nbins_srcpos,minsrcpos,maxsrcpos); 318 319 coszenithTitle = Form("%s%02i","hCosZenithOn",numberTimeBins-1); 320 hCosZenith_on_timebin = new TH1F (coszenithTitle,"",nbins_coszenith,mincoszenith,maxcoszenith); 321 322 timediffTitle = Form("%s%02i","hTimeDiffOn",numberTimeBins-1); 323 hTimeDiff_on_timebin = new TH1F (timediffTitle,"",nbins_timediff,mintimediff,maxtimediff); 324 325 } 326 327 } 501 328 502 } 503 504 } 505 506 507 //Compute the time on for last time bin 508 timeOn[numberTimeBins-1] += (mjdLastEvent-mjdFirstEvent)*kDay; 509 widthTimeBinOn[numberTimeBins-1] = (mjd-mjdFirstEventofBin)/2; 510 meanTimeBinOn[numberTimeBins-1] = mjdFirstEventofBin + widthTimeBinOn[numberTimeBins-1]; 329 loop_on.PostProcess(); 330 331 tlist_on.PrintStatistics(); 332 333 timeOn[numberTimeBins-1] = livetime_on.GetLiveTimeTArray().At(numberTimeBins-1); 334 meanTimeBinOn[numberTimeBins-1] = livetime_on.GetMeanRealTimeTArray().At(numberTimeBins-1); 335 widthTimeBinOn[numberTimeBins-1] = livetime_on.GetWidthRealTimeTArray().At(numberTimeBins-1); 336 511 337 //Compute the number of on events for the last time bin 512 numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90+1); 513 numberBkgEventsToNormOn[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90+1); 514 515 hTimeDiff_on_timebin->Fit("exp","RQ0"); 516 deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries()); 517 518 cout.precision(5); 519 cout << "2-> Exp(" << f1->GetParameter(0) << " + " << f1->GetParameter(1) << "*x) +- [" << f1->GetParError(0) << ' ' << f1->GetParError(1) << "]" << endl; 520 cout << "Calc entries " << (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/TMath::Abs(f1->GetParameter(1)) << " +- Nt*(" << TMath::Sqrt( f1->GetParameter(0)*f1->GetParError(0)*f1->GetParameter(0)*f1->GetParError(0) + (f1->GetParError(1)*f1->GetParError(1))/(f1->GetParameter(1)*f1->GetParameter(1))) << ") meas entries " << hTimeDiff_on_timebin->GetEntries() << " dead fraction " << deadFractionOn[numberTimeBins-1] << endl; 521 522 cout.precision(10); 523 cout << "timeOn[" << numberTimeBins-1 << "] " << timeOn[numberTimeBins-1] << " mjdLastEvent " << mjdLastEvent << " mjdFirstEvent " << mjdFirstEvent << endl; 524 cout << "mjd " << mjd << " mjdFirstEventofBin " << mjdFirstEventofBin << " meanTimeBinOn " << meanTimeBinOn[numberTimeBins-1] << " widthTimeBinOn " << widthTimeBinOn[numberTimeBins-1] << endl; 338 numberOnEvents[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90); 339 numberBkgEventsToNormOn[numberTimeBins-1] = (Double_t) hAlpha_on_abs_timebin->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90); 340 341 // hTimeDiff_on_timebin->Fit("exp","RQ0"); 342 // deadFractionOn[numberTimeBins-1] = (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/(TMath::Abs(f1->GetParameter(1))*hTimeDiff_on_timebin->GetEntries()); 343 344 // cout.precision(5); 345 // cout << "2-> Exp(" << f1->GetParameter(0) << " + 346 // cout << "Calc entries " << (nbins_timediff/(maxtimediff-mintimediff))*TMath::Exp(f1->GetParameter(0))/TMath::Abs(f1->GetParameter(1)) << " +- Nt*(" << TMath::Sqrt( f1->GetParameter(0)*f1->GetParError(0)*f1->GetParameter(0)*f1->GetParError(0) + (f1->GetParError(1)*f1->GetParError(1))/(f1->GetParameter(1)*f1->GetParameter(1))) << ") meas entries " << hTimeDiff_on_timebin->GetEntries() << " dead fraction " << deadFractionOn[numberTimeBins-1] << endl; 347 deadFractionOn[numberTimeBins-1] = 1.; 525 348 526 349 alphaHistoOn.AddLast(hAlpha_on_abs_timebin); … … 529 352 timediffHistoOn.AddLast(hTimeDiff_on_timebin); 530 353 531 // cout << "hAlpha_on_abs_timebin " << hAlpha_on_abs_timebin << " alphaHistoOn [" << numberTimeBins-1 << "] " << alphaHistoOn[numberTimeBins-1] << endl;532 533 loop_on.PostProcess();534 535 tlist_on.PrintStatistics();536 537 for (UInt_t bin=0; bin<numberTimeBins; bin++)538 cout << bin << " timeOn " << timeOn[bin] << " min-max zenithOn " << zenithMinimumOn[bin] << "-" << zenithMaximumOn[bin] << " min-max cos(zenithOn) " << TMath::Cos(zenithMinimumOn[bin]*kConvDegToRad) << "-" << TMath::Cos(zenithMaximumOn[bin]*kConvDegToRad) << " numberOnEvents " << numberOnEvents[bin] << " numberOnEventsAfterCleaning " << numberOnEventsAfterCleaning[bin] << " deadFractionOn " << deadFractionOn[bin] << endl;539 540 541 542 543 354 //-----------------------OFF------------------------ 544 355 … … 552 363 TArrayD numberBkgEventsToNormOff(numberTimeBins); 553 364 TArrayD onOffNormFactor(numberTimeBins); 365 TArrayD onOffNormFactorWithEvents(numberTimeBins); 366 TArrayD onOffNormFactorWithLiveTime(numberTimeBins); 554 367 555 368 TArrayD numberExcessEvents(numberTimeBins); … … 590 403 MRawEvtHeader evtheader_off; 591 404 MTime time_off; 405 MLiveTime livetime_off; 592 406 MHillas hillas_off; 593 407 MHillasSrc hillassrc_off; … … 599 413 plist_off.AddToList(&evtheader_off); 600 414 plist_off.AddToList(&time_off); 415 plist_off.AddToList(&livetime_off); 601 416 plist_off.AddToList(&hillas_off); 602 417 plist_off.AddToList(&hillassrc_off); … … 627 442 MHillasSrcCalc csrc_off; 628 443 629 //tasklist 444 MLiveTimeCalc livetimecalc_off; 445 446 //tasklist 630 447 tlist_off.AddToList(&read_off); 448 tlist_off.AddToList(&livetimecalc_off); 631 449 tlist_off.AddToList(&srcplace_timebin); // This is just to run the preprocess correctely 632 450 tlist_off.AddToList(&csrc_off); // This is just to run the preprocess correctely … … 714 532 } 715 533 } 716 717 // if (zenith > zenithMaximumOn[numberTimeBins-1])718 // {719 // cout << "Exit loop because we arrive to zenith = " << zenith << endl;720 // break;721 // }722 534 } 723 535 … … 734 546 for (UInt_t bin=0; bin<numberTimeBins; bin++) 735 547 { 736 548 cout.precision(5); 549 cout << bin << " timeOn " << timeOn[bin] << " mean-width time " << meanTimeBinOn[bin] << "-" << widthTimeBinOn[bin] << endl; 550 } 551 livetime_off.Print("last"); 552 livetime_off.Print("all"); 553 cout << "livetime_off.GetLiveTime() " << livetime_off.GetLiveTime() << endl; 554 555 for (UInt_t bin=0; bin<numberTimeBins; bin++) 556 { 557 timeOff[bin] = livetime_off.GetLiveTime(); 558 559 // 737 560 meanTriggerRateOn[bin] = numberOnEventsAfterCleaning[bin]/timeOn[bin]; 738 561 errorMeanTriggerRateOn[bin] = TMath::Sqrt(numberOnEventsAfterCleaning[bin])/timeOn[bin]; … … 745 568 widthTimeBinOnInSec[bin] = widthTimeBinOn[bin]*kDay; 746 569 747 numberOffEvents[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90 +1);748 numberBkgEventsToNormOff[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90 +1);570 numberOffEvents[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90); 571 numberBkgEventsToNormOff[bin] = (Double_t) ((TH1F*)alphaHistoOff[bin])->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90); 749 572 if (numberOffEvents[bin] != 0 && numberBkgEventsToNormOff[bin] != 0) 750 573 { 751 onOffNormFactor[bin] = numberBkgEventsToNormOn[bin]/numberBkgEventsToNormOff[bin]; 574 onOffNormFactorWithEvents[bin] = numberBkgEventsToNormOn[bin]/numberBkgEventsToNormOff[bin]; 575 onOffNormFactorWithLiveTime[bin] = timeOn[bin]/timeOff[bin]; 576 onOffNormFactor[bin] = onOffNormFactorWithEvents[bin]; 752 577 numberExcessEvents[bin] = numberOnEvents[bin] - onOffNormFactor[bin]*numberOffEvents[bin]; 753 578 errorNumberExcessEvents[bin] = TMath::Sqrt(numberOnEvents[bin] + onOffNormFactor[bin]*onOffNormFactor[bin]*numberOffEvents[bin]); … … 756 581 numberExcessEventsPerMin[bin] = 60.*numberExcessEvents[bin]/timeOn[bin]*(deadFractionOn[bin]>1.?deadFractionOn[bin]:1.); 757 582 errorNumberExcessEventsPerMin[bin] = 60.*errorNumberExcessEvents[bin]/timeOn[bin]; 583 758 584 } 759 585 } … … 762 588 { 763 589 cout.precision(5); 764 cout << bin << " timeOn " << timeOn[bin] << " mean-width time " << meanTimeBinOn [bin] << "-" << widthTimeBinOn[bin] << endl;590 cout << bin << " timeOn " << timeOn[bin] << " mean-width time " << meanTimeBinOnInSec[bin] << "-" << widthTimeBinOnInSec[bin] << endl; 765 591 cout << " numberOnEvents\t " << numberOnEvents[bin] << endl; 766 592 cout << " numberOffEvents\t " << numberOffEvents[bin] << endl; 767 593 cout << " numberBkgEventsToNormOn\t " << numberBkgEventsToNormOn[bin] << endl; 768 594 cout << " numberBkgEventsToNormOff\t " << numberBkgEventsToNormOff[bin] << endl; 769 cout << " onOffNormFactor\t " << onOffNormFactor[bin] << endl; 595 cout << " onOffNormFactorWithEvents\t " << onOffNormFactorWithEvents[bin] << endl; 596 cout << " onOffNormFactorWithLiveTime\t " << onOffNormFactorWithLiveTime[bin] << endl; 770 597 cout << " numberExcessEvents\t " << numberExcessEvents[bin] << " +- " << errorNumberExcessEvents[bin] << endl; 771 598 cout << " deadFraction\t" << deadFractionOn[bin] << endl; … … 796 623 TGraphErrors* cosmicrategraph = new TGraphErrors(numberTimeBins,meanTimeBinOnInSec.GetArray(),meanTriggerRateOn.GetArray(),widthTimeBinOnInSec.GetArray(),errorMeanTriggerRateOn.GetArray()); 797 624 cosmicrategraph->SetTitle("Cosmic Rate"); 798 cosmicrategraph->SetMinimum(0.);799 625 cosmicrategraph->SetMarkerStyle(21); 800 626 cosmicrategraph->SetMarkerSize(0.03); … … 810 636 Float_t maxAlphaHistoHeight = 0; 811 637 812 for (UInt_t bin=0; bin<numberTimeBins -1; bin++)638 for (UInt_t bin=0; bin<numberTimeBins; bin++) 813 639 { 814 640 for (UInt_t i=1; i<=nbins_abs; i++) … … 816 642 maxAlphaHistoHeight = ((TH1F*)alphaHistoOn[bin])->GetBinContent(i); 817 643 } 818 644 645 cout << "maxAlphaHistoHeight " << maxAlphaHistoHeight << endl; 819 646 820 647 for (UInt_t bin=0; bin<numberTimeBins-1; bin++) … … 897 724 // ############################################################################ 898 725 899 Double_t norm_on_abs = (Double_t) hAlpha_on_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90 +1);900 Double_t exces_on_abs = (Double_t) hAlpha_on_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90 +1);901 Double_t norm_off_abs = (Double_t) hAlpha_off_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90 +1);902 Double_t exces_off_abs = (Double_t) hAlpha_off_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90 +1);726 Double_t norm_on_abs = (Double_t) hAlpha_on_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90); 727 Double_t exces_on_abs = (Double_t) hAlpha_on_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90); 728 Double_t norm_off_abs = (Double_t) hAlpha_off_abs->Integral((Int_t)bkgnormmin*nbins_abs/90+1,(Int_t)bkgnormmax*nbins_abs/90); 729 Double_t exces_off_abs = (Double_t) hAlpha_off_abs->Integral((Int_t)sigexccmin*nbins_abs/90+1,(Int_t)sigexccmax*nbins_abs/90); 903 730 Double_t norm = norm_on_abs/norm_off_abs; 904 731 … … 1296 1123 1297 1124 1298 1125 Double_t ChiSquareNDof(TH1D *h1, TH1D *h2) 1126 { 1127 Double_t chiq = 0.; 1128 Double_t chi; 1129 Double_t error; 1130 Int_t nbinsnozero = 0; 1131 1132 Int_t nbins = h1->GetNbinsX(); 1133 if (nbins != h2->GetNbinsX() || nbins == 0) 1134 return -1; 1135 1136 for (UInt_t bin=1; bin<=nbins; bin++) 1137 { 1138 error = sqrt(h1->GetBinError(bin)*h1->GetBinError(bin) + 1139 h2->GetBinError(bin)*h2->GetBinError(bin)); 1140 if (error != 0) 1141 { 1142 chi = (h1->GetBinContent(bin)-h2->GetBinContent(bin))/error; 1143 chiq += chi*chi; 1144 nbinsnozero++; 1145 } 1146 } 1147 1148 return (nbinsnozero>0?chiq/nbinsnozero:0); 1149 } 1150 1151 Int_t GetBin(Double_t size, Int_t numberSizeBins, Double_t sizeBins[]) 1152 { 1153 1154 Int_t result = -1; 1155 1156 Int_t lowerbin = 0; 1157 Int_t upperbin = numberSizeBins; 1158 Int_t bin; 1159 1160 Int_t count = 0; 1161 1162 if (size >= sizeBins[0]) 1163 { 1164 while (upperbin - lowerbin > 1 && count++<=numberSizeBins) 1165 { 1166 bin = (upperbin+lowerbin)/2; 1167 if (size >= sizeBins[bin]) 1168 lowerbin = bin; 1169 else 1170 upperbin = bin; 1171 } 1172 result = count<=numberSizeBins?lowerbin:-1; 1173 } 1174 1175 return result; 1176 1177 } 1178 -
trunk/MagicSoft/Mars/mtemp/mifae/macros/makeHillas.C
r4310 r4429 71 71 72 72 MJPedestal pedloop0; 73 pedloop0.SetOutputPath("./"); 73 74 pedloop0.SetInput(&pedcaliter); 74 75 pedloop0.SetExtractor(&extractor); … … 86 87 MJCalibration calloop; 87 88 89 calloop.SetOutputPath("./"); 88 90 calloop.SetExtractor(&extractor); 89 91 calloop.SetInput(&caliter); … … 101 103 // First Compute the pedestals 102 104 MJPedestal pedloop; 105 pedloop.SetOutputPath("./"); 103 106 pedloop.SetInput(&pediter); 104 107 … … 225 228 read4.AddTree("Drive"); 226 229 static_cast<MRead&>(read4).AddFiles(datiter); 230 read4.AddToBranchList("MReportCurrents.*"); 227 231 228 232 // Task that use Slow control information … … 252 256 MSrcPosFromStars srcposfromstar; 253 257 254 // set bad pixels 255 MBlindPixelCalc blind; 256 MBlindPixelCalc blind2; 257 const Short_t x[16] = {330,395,329,396,389, 258 323,388,322,384,385, 259 386,387,321,320,319, 260 394}; 261 const TArrayS bp(16,(Short_t*)x); 262 blind.SetPixelIndices(bp); 263 blind2.SetPixelIndices(bp); 258 MBadPixelsTreat interpolatebadpixels; 259 interpolatebadpixels.SetUseInterpolation(); 264 260 265 261 MImgCleanStd clean(lcore,ltail); … … 315 311 tlist4.AddToList(&extractor,"Events"); 316 312 tlist4.AddToList(&photcalc,"Events"); 317 //tlist4.AddToList(&blind,"Events"); 313 if(calflag==11) 314 tlist4.AddToList(&interpolatebadpixels); 318 315 tlist4.AddToList(&clean,"Events"); 319 316 … … 332 329 //tlist4.AddToList(&blind2,"Events"); 333 330 tlist4.AddToList(&hcalc,"Events"); 334 // tlist4.AddToList(&srcposcalc,"Events");335 331 tlist4.AddToList(&csrc1,"Events"); 336 332 tlist4.AddToList(&write,"Events"); … … 350 346 351 347 tlist4.PrintStatistics(); 348 352 349 } 353 350 //------------------------------------------------------------------------------- … … 476 473 } 477 474 478 479 480 481 475 pedcaliter.Reset(); 476 pediter.Reset(); 477 caliter.Reset(); 478 datiter.Reset(); 482 479 TString pfile; 483 480 -
trunk/MagicSoft/Mars/mtemp/mifae/macros/observationTime.C
r4259 r4429 6 6 Int_t iniRun=0, 7 7 Int_t finRun=99999, 8 TString filename="/mnt/users/jrico/magic/mars/mars/mtemp/mifae/hillas/AllCrabNebula_3015_kDefault_MispE.root" ,8 TString filename="/mnt/users/jrico/magic/mars/mars/mtemp/mifae/hillas/AllCrabNebula_3015_kDefault_MispE.root" 9 9 ) 10 10 {
Note:
See TracChangeset
for help on using the changeset viewer.