BmnRoot
Loading...
Searching...
No Matches
utilites.h
Go to the documentation of this file.
1// utilities for ROOT macro analysis
2//Author: Dmitry Finogeev, dmitry.finogeev@cern.ch
3
4
5#ifndef UTILITES_H
6#define UTILITES_H
7//header what include most popular root classes
8
9//#include "CanvasShow/canvasshow.h"
10//#include "CanvasShow/canvasshow.cxx"
11
12#include <TObjArray.h>
13#include <TFile.h>
14#include <TChain.h>
15#include <TObjString.h>
16#include <TH1.h>
17#include <TH2.h>
18#include <TF1.h>
19
20#include<TROOT.h>
21#include<TSystem.h>
22#include<TSystemDirectory.h>
23#include<TSystemFile.h>
24#include<TList.h>
25#include<TCollection.h>
26#include<TMath.h>
27#include<TDirectory.h>
28#include<TStopwatch.h>
29#include<TGraph.h>
30#include<TBrowser.h>
31#include<TCanvas.h>
32#include<TProfile.h>
33#include<TStyle.h>
34#include<TVectorD.h>
35
36
37#include <iostream>
38#include <fstream>
39
40#define OUTINFILE 0
41
42#if OUTINFILE
43std::ofstream current_out ("CurrentLog.txt");
44#else
45std::ostream &current_out = std::cout;
46#endif
47
48#define uVERBOSE 1 // 1 - errors; 2 - results; 3 - details; 4 - debugging; 5 - flood
49#define fit_funk_name "fitFunction"
50
51using namespace std;
52
53Int_t FilesListFromDir(const TString filespath, TObjArray *const collection, Int_t max_num, TString content);
54Int_t SaveArrayStructInFile(TCollection *array, TFile *file, TString start_path = "");
55Int_t SaveArrayOfObjectsArraysInFile(TObjArray *ObjectsArray = 0, const TString filename = "data", const TString path = "./", const Int_t deleteArrays = 0);
56
57Int_t FitHistogramm(TH1 *histogramm, Double_t &mean, Double_t &sigma, Double_t RangeXmin = 0., Double_t RangeXmax = 0., TString fit_formula = "gaus");
58Int_t FitHistogrammInFirstMaxPeak(TH1* histogramm, Double_t &Mean, Double_t &Sigma, Double_t peaksRangeMin,
59 Double_t peaksRangeMax, Double_t peakRange, Double_t peakRMSRange = 1.5);
60Int_t FitHistogrammInRMSRange(TH1 *histogramm, Double_t &mean, Double_t &sigma, Double_t RMSrange, Double_t RangeXmin = 0, Double_t RangeXmax = 0);
61Int_t FitHistogrammInMaxPeak(TH1 *histogramm, Double_t &mean, Double_t &sigma, Double_t PeakWidth, Double_t RangeXmin = 0., Double_t RangeXmax = 0.);
62
63Int_t ArrayStat(const Int_t nElements, const Double_t* Array, Double_t &maxElement, Double_t &minElement, Double_t &Amplitude, Double_t &Mean, Double_t &RMS);
64void ResetArray(TObjArray **array, const TString option);
65
66
67Int_t ConstructFilesListFromDirectory(const TString &Directory, TCollection *const collection)
68{
69 TString filespath = Directory;
70 if(!filespath.EndsWith("/")) filespath += "/";
71
72 if(!collection)
73 {
74 if(uVERBOSE >= 1) current_out << Form("utilites.h/ConstructFilesListInDirectory: collection ptr is NULL")<<std::endl;
75 return 0;
76 }
77
78 TSystemFile *currFile;
79 TString currFileName;
80 TSystemDirectory SDirectory(filespath, filespath);
81 TIter nextfile(SDirectory.GetListOfFiles());
82
83 Int_t nfiles = 0;
84 while ((currFile=(TSystemFile*)nextfile()))
85 {
86 currFileName = filespath + (currFile->GetName());
87 if (!currFile->IsDirectory() && currFileName.EndsWith("root"))
88 {
89 collection->Add( new TObjString( currFileName.Data() ) );
90
91 if(uVERBOSE >= 3) current_out << Form("utilites.h/ConstructFilesListInDirectory: found root file \"%s\"",
92 currFileName.Data() )<<std::endl;
93 nfiles++;
94 }
95 }
96
97 return nfiles;
98}
99
101Int_t RUNnumberInString(const TString &str)
102{
103 const Int_t lenght = 6;
104 Int_t value;
105 TString sub_str;
106 for(Int_t start=0; start <= str.Sizeof() - lenght; start++)
107 {
108 sub_str = str(start, lenght);
109 value = sub_str.Atoi();
110 //printf("debug:start = %i; strpath = %s; value = %i\n",start, sub_str.Data(), value);
111 if( value < 999999 && 100000 < value) return value;
112 }
113
114 return 0;
115}
116
118TObject* CreateXprofile(const TH1 *sourcehist)
119{
120 TProfile *profile = ((TH2*)sourcehist)->ProfileX( Form("%s_prof",sourcehist->GetName()) , 1, -1, "o");
121 profile->SetTitle( sourcehist->GetTitle() );
122 if(uVERBOSE >= 2) current_out << Form("utilites.h/CreateXprofile: profile \"%s\" was created",sourcehist->GetName())<<std::endl;
123
124 return (TObject*)profile;
125}
126
128Int_t CollectRAWfilesInChain(TChain* rawTree, const TString dirname, const Int_t maxNfiles = 50)
129{
130
131 TString filespath = dirname;
132 if(!filespath.EndsWith("/")) filespath += "/";
133
134 if(uVERBOSE >= 2) current_out << Form("utilites.h/CollectRAWfilesInChain: collecting files from \"%s\"",filespath.Data())<<std::endl;
135
136 TSystemFile *currFile;
137 TString currFileName;
138 TSystemDirectory SDirectory(filespath, filespath);
139 TIter nextfile(SDirectory.GetListOfFiles());
140
141 Int_t NfilesLoaded = 0;
142 while ((currFile=(TSystemFile*)nextfile()))
143 {
144 currFileName = filespath + (currFile->GetName());
145 if (!currFile->IsDirectory() && currFileName.EndsWith("root"))
146 {
147 rawTree->AddFile( currFileName.Data(), 0, "t0tree" );
148 if(uVERBOSE >= 4) current_out << Form("utilites.h/CollectRAWfilesInChain: %i file \"%s\" was added in chain",
149 NfilesLoaded,currFileName.Data())<<std::endl;
150 NfilesLoaded++;
151 }
152
153 if(maxNfiles <= NfilesLoaded) break;
154 }
155
156 if(uVERBOSE >= 2) current_out << Form("utilites.h/CollectRAWfilesInChain: collected %d files, total events: %lld",NfilesLoaded, rawTree->GetEntries())<<std::endl;
157
158 return NfilesLoaded;
159
160}
161
163Int_t FilesListFromDir(const TString Directory, TObjArray *const collection, Int_t max_num, TString content)
164{
165 TString filespath = Directory;
166
167 if(!filespath.EndsWith("/")) filespath += "/";
168
169 if(!collection)
170 {
171 current_out << Form("ERROR: Collection ptr is NULL") << std::endl;
172 return 0;
173 }
174
175 TSystemFile *currFile;
176 TString currFileName;
177 TSystemDirectory SDirectory(filespath, filespath);
178 TIter nextfile(SDirectory.GetListOfFiles());
179
180 Int_t nfiles = 0;
181 while ((currFile=(TSystemFile*)nextfile()))
182 {
183 currFileName = currFile->GetName();
184 if (!currFile->IsDirectory() && currFileName.Contains(content))
185 {
186 collection->Add( new TObjString( currFileName.Data() ) );
187 nfiles++;
188 if(max_num > 0)
189 if(nfiles == max_num) break;
190 }
191 }
192
193 return nfiles;
194}
195
196
198Int_t SaveArrayStructInFile(TCollection *array, TFile *file, TString start_path)
199{
200 //if(start_path = "") start_path = "/";
201 //file->cd(start_path.Data());
202
203 //file->mkdir(array->GetName());
204
205 if(array == 0) return 0;
206
207 if(start_path == "")
208 {
209 start_path = (TString)array->GetName();
210 }
211 else
212 {
213 start_path += "/" + (TString)array->GetName();
214 }
215
216 file->mkdir(start_path.Data());
217 file->cd(start_path.Data());
218
219 TIter nx_iter((TCollection*)(array));
220 TObject* obj_ptr;
221 while ( (obj_ptr=(TObjArray*)nx_iter()) )
222 {
223 if(obj_ptr == 0) continue;
224
225 if(TCollection *new_coll = dynamic_cast<TCollection*>(obj_ptr))
226 {
227 SaveArrayStructInFile(new_coll, file, start_path);
228 file->cd(start_path.Data());
229 }
230 else
231 {
232 if(uVERBOSE >= 2) current_out << Form("Writing %s object", obj_ptr->GetName()) << std::endl;
233 obj_ptr->Write();
234 }
235 }
236
237 return 1;
238}
239
240
242Int_t SaveArrayOfObjectsArraysInFile(TObjArray *ObjectsArray,const TString filename, const TString path, const Int_t deleteArrays)
243{
244 // TList *ObjectsArray = (TList*)ObjectsArray_;
245 if( !ObjectsArray )
246 {
247 if(uVERBOSE >= 1) current_out <<Form("utilites.h/SaveArrayObjectsArrayInFile::ERROR: Array does NOT exist")<<std::endl;
248 return 0;
249 }
250
251 if(ObjectsArray->GetLast() < 0)
252 {
253 if(uVERBOSE >= 1) current_out <<Form("utilites.h/SaveArrayObjectsArrayInFile: WARING: Objects array is empty")<<std::endl;
254 return 0;
255 }
256
257 TString resultfilename = filename;
258 if(!resultfilename.EndsWith(".root")) resultfilename += ".root";
259
260 TObject* pobjCurrObject = NULL;
261 TString objectname;
262
263 TFile *file = new TFile(resultfilename.Data(), "UPDATE");
264
265 if(file->IsOpen())
266 {
267 if(uVERBOSE >= 2) current_out <<Form("utilites.h/SaveArrayObjectsArrayInFile: Writing arrays in file \"%s\" ... ", resultfilename.Data())<<std::endl;
268 }
269 else
270 {
271 if(uVERBOSE >= 1) current_out <<Form("utilites.h/SaveArrayObjectsArrayInFile::ERROR: File \"%s\" was NOT opened for UPDATE",resultfilename.Data())<<std::endl;
272 return 0;
273 }
274
275
276 for(Int_t n=0; n<ObjectsArray->GetEntries(); n++)
277 {
278 pobjCurrObject = ObjectsArray->At(n);
279 if(pobjCurrObject == NULL) continue;
280
281 objectname = pobjCurrObject->GetName();
282 if(uVERBOSE >= 3) current_out <<Form("utilites.h/SaveArrayObjectsArrayInFile::Saving \"%s\" ...",objectname.Data())<<std::endl;
283 file->rmdir(objectname.Data()); //delete old data to rewrite
284 TDirectory *CurrDir = file->mkdir(objectname.Data()); CurrDir->cd();
285 pobjCurrObject->Write();
286
287 if(uVERBOSE >= 2) current_out <<Form("utilites.h/SaveArrayObjectsArrayInFile::Object \"%s\" was successfully writed",objectname.Data())<<std::endl;
288 }
289
290 file->Close();
291 delete file;
292
293 gDirectory->cd("Rint:/");
294
295 TObjArray *currArray = NULL;
296 if(deleteArrays) //
297 {
298 for(Int_t n=0; n<ObjectsArray->GetEntries(); n++)
299 {
300 currArray = (TObjArray*)ObjectsArray->At(n);
301 if(currArray == NULL) continue;
302 objectname = currArray->GetName();
303 if(uVERBOSE >= 2) current_out <<Form("utilites.h/SaveArrayObjectsArrayInFile::Deleting \"%s\"",objectname.Data())<<std::endl;
304 ResetArray(&currArray, "delete");
305 }
306 ObjectsArray->Clear();
307 }
308
309 return 1;
310}
311
313Int_t FitHistogrammInFirstMaxPeak(TH1* histogramm, Double_t &Mean, Double_t &Sigma, Double_t peaksRangeMin,
314 Double_t peaksRangeMax, Double_t peakRange, Double_t peakRMSRange)
315{
316 const Double_t peakLeftMinValueRatio = 0.5;
317
318 histogramm->SetAxisRange(peaksRangeMin, peaksRangeMax); //cutting saturation peak
319 Double_t PeakMax = histogramm->GetBinCenter( histogramm->GetMaximumBin() );
320 Double_t PeakMaxValue = histogramm->GetBinContent(histogramm->GetMaximumBin());
321 Double_t PeakLeftEdge = histogramm->GetBinCenter(histogramm->FindFirstBinAbove(PeakMaxValue*peakLeftMinValueRatio)-2);
322 histogramm->SetAxisRange(PeakLeftEdge, PeakMax+peakRange); //cutting saturation peak
323 Double_t PeakRMS = histogramm->GetRMS();
324
325 FitHistogramm(histogramm, Mean, Sigma, PeakLeftEdge-PeakRMS*0.5, PeakMax+PeakRMS*peakRMSRange);
326
327 histogramm->SetAxisRange(peaksRangeMin, peaksRangeMax); //cutting saturation peak
328
329 return 1;
330}
331
333Int_t FitHistogrammInRMSRange(TH1 *histogramm, Double_t &mean, Double_t &sigma, Double_t RMSrange, Double_t RangeXmin, Double_t RangeXmax)
334{
335 if( !histogramm ){if(uVERBOSE >= 1)
336 current_out <<Form("utilites.h/FitHistogrammInRMSRange::ERROR: histogramm pointer is NULL")<<std::endl; return 0;}
337
338 if(RangeXmax <= RangeXmin)
339 {
340 RangeXmin = histogramm->GetXaxis()->GetXmin();
341 RangeXmax = histogramm->GetXaxis()->GetXmax();
342 }
343
344 Double_t hist_start_Xmin = histogramm->GetXaxis()->GetXmin();
345 Double_t hist_start_Xmax = histogramm->GetXaxis()->GetXmax();
346
347 histogramm->SetAxisRange(RangeXmin, RangeXmax);
348 Double_t histMean = histogramm->GetMean();
349 Double_t histRMS = histogramm->GetRMS();
350 histogramm->SetAxisRange(hist_start_Xmin, hist_start_Xmax);
351
352 Double_t fit_range_min = histMean-RMSrange*histRMS;
353 Double_t fit_range_max = histMean+RMSrange*histRMS;
354
355 if(fit_range_min < RangeXmin) fit_range_min = RangeXmin;
356 if(fit_range_max > RangeXmax) fit_range_max = RangeXmax;
357
358 return FitHistogramm(histogramm, mean, sigma, fit_range_min, fit_range_max);
359}
361Int_t FitHistogrammInMaxPeak(TH1 *histogramm, Double_t &mean, Double_t &sigma, Double_t PeakWidth, Double_t RangeXmin, Double_t RangeXmax)
362{
363 if( !histogramm ){if(uVERBOSE >= 1)
364 current_out <<Form("utilites.h/FitHistogrammInMaxPeak::ERROR: histogramm pointer is NULL")<<std::endl; return 0;}
365
366
367 if(RangeXmax <= RangeXmin)
368 {
369 RangeXmin = histogramm->GetXaxis()->GetXmin();
370 RangeXmax = histogramm->GetXaxis()->GetXmax();
371 }
372
373 Double_t hist_start_Xmin = histogramm->GetXaxis()->GetXmin();
374 Double_t hist_start_Xmax = histogramm->GetXaxis()->GetXmax();
375
376 histogramm->SetAxisRange(RangeXmin, RangeXmax);
377 Double_t PeakMax = histogramm->GetBinCenter( histogramm->GetMaximumBin() );
378 histogramm->SetAxisRange(hist_start_Xmin, hist_start_Xmax);
379
380 Double_t FitRangeMin = PeakMax-PeakWidth*0.5;
381 Double_t FitRangeMax = PeakMax+PeakWidth*0.5;
382
383 FitRangeMin = (FitRangeMin < RangeXmin)? RangeXmin : FitRangeMin;
384 FitRangeMax = (FitRangeMax > RangeXmax)? RangeXmax : FitRangeMax;
385
386
387 if(uVERBOSE >= 4)current_out<<Form("utilites.h/FitHistogrammInMaxPeak: Hsit: \"%s\" Maximum: %.2f, Peak Width: %.2f, Range:[%.2f, %.2f]",
388 histogramm->GetName(), PeakMax, PeakWidth, FitRangeMin, FitRangeMax)<<std::endl;
389
390
391 return FitHistogramm(histogramm, mean, sigma,FitRangeMin,FitRangeMax);
392}
393
395Int_t FitHistogrammSmartInMaxPeak(TH1 *histogramm, Double_t &mean, Double_t &sigma,Double_t nRMSfit = 2., Double_t RangeXmin = 0., Double_t RangeXmax = 0.)
396{
397 static const Int_t nSmooth = 1;
398 static const Double_t nSmoothRMScut = 1.;
399 //static const Double_t nRMSfit = 2;
400
401 if( !histogramm ){if(uVERBOSE >= 1)
402 current_out <<Form("utilites.h/FitHistogrammSmartInMaxPeak::ERROR: histogramm pointer is NULL")<<std::endl; return 0;}
403
404 if(RangeXmax <= RangeXmin)
405 {
406 RangeXmin = histogramm->GetXaxis()->GetXmin();
407 RangeXmax = histogramm->GetXaxis()->GetXmax();
408 }
409
410 histogramm->SetAxisRange(RangeXmin, RangeXmax);
411
412 //searching fit ranges
413 TH1 *tempSmoothProjection = (TH1*)histogramm->Clone("tempSmoothFitHist");
414 tempSmoothProjection->Smooth(nSmooth);
415 tempSmoothProjection->SetAxisRange(RangeXmin, RangeXmax);
416 Double_t SmoothPeakMax = tempSmoothProjection->GetBinCenter( tempSmoothProjection->GetMaximumBin() );
417 Double_t SmoothRMS = tempSmoothProjection->GetRMS();
418
419 //maximum on hist edge:
420 if( TMath::Abs(SmoothPeakMax - RangeXmax) < SmoothRMS*0.5 )
421 {
422 RangeXmax += SmoothRMS*0.5;
423 tempSmoothProjection->SetAxisRange(RangeXmin, RangeXmax);
424 SmoothPeakMax = tempSmoothProjection->GetBinCenter( tempSmoothProjection->GetMaximumBin() );
425 SmoothRMS = tempSmoothProjection->GetRMS();
426 } else if( TMath::Abs(SmoothPeakMax - RangeXmin) < SmoothRMS*0.5 )
427 {
428 RangeXmin -= SmoothRMS*0.5;
429 tempSmoothProjection->SetAxisRange(RangeXmin, RangeXmax);
430 SmoothPeakMax = tempSmoothProjection->GetBinCenter( tempSmoothProjection->GetMaximumBin() );
431 SmoothRMS = tempSmoothProjection->GetRMS();
432 }
433
434 tempSmoothProjection->SetAxisRange(SmoothPeakMax-SmoothRMS*nSmoothRMScut, SmoothPeakMax+SmoothRMS*nSmoothRMScut);
435 Double_t SmoothRMSreal = tempSmoothProjection->GetRMS();
436 delete tempSmoothProjection;
437
438 Double_t FitRangeMin = SmoothPeakMax-SmoothRMSreal*nRMSfit;
439 Double_t FitRangeMax = SmoothPeakMax+SmoothRMSreal*nRMSfit;
440
441 FitRangeMin = (FitRangeMin < RangeXmin)? RangeXmin : FitRangeMin;
442 FitRangeMax = (FitRangeMax > RangeXmax)? RangeXmax : FitRangeMax;
443
444 if(uVERBOSE >= 3)current_out<<Form("utilites.h/FitHistogrammSmartInMaxPeak: Hsit: \"%s\" Maximum: %.2f, SmoothRMS: %.2f, SmoothRMSreal: %.2f, Range:[%.2f, %.2f]",
445 histogramm->GetName(), SmoothPeakMax, SmoothRMS, SmoothRMSreal, FitRangeMin, FitRangeMax)<<std::endl;
446 //return NULL;
447 return FitHistogramm(histogramm, mean, sigma, FitRangeMin, FitRangeMax);
448}
449
451Int_t FitHistogramm(TH1 *histogramm, Double_t &mean, Double_t &sigma, Double_t RangeXmin, Double_t RangeXmax, TString fit_formula)
452{
453 static const Int_t minFitEvents = 1; //minimum events in range for fitting
454 static const Int_t EventsNoFit = 15;
455
456 TString fitoptions[] = {"QR", "QRL", "QRI", "END"};
457
458 if( !histogramm ){if(uVERBOSE >= 1)
459 current_out <<Form("utilites.h/FitHistogramm::ERROR: histogramm pointer is NULL")<<std::endl; return 0;}
460
461
462 if(RangeXmax <= RangeXmin)
463 {
464 RangeXmin = histogramm->GetXaxis()->GetXmin();
465 RangeXmax = histogramm->GetXaxis()->GetXmax();
466 }
467
468
469 //if too low events for fitting returning 0
470 Int_t intBinMin = histogramm->FindBin(RangeXmin);
471 Int_t intBinMax = histogramm->FindBin(RangeXmax);
472 Int_t integral = histogramm->Integral(intBinMin, intBinMax);
473
474 if( integral < minFitEvents){if(uVERBOSE >= 1)
475 current_out <<Form("utilites.h/FitHistogramm::ERROR: Too low events in \"%s\" histogram: %i",histogramm->GetName(), integral)<<std::endl; return 0;}
476
477 if(uVERBOSE >= 2)
478 current_out <<Form("utilites.h/FitHistogramm::Fiting histogramm \"%s\"",histogramm->GetName())<<std::endl;
479
480
481 Double_t hist_start_Xmin = histogramm->GetXaxis()->GetXmin();
482 Double_t hist_start_Xmax = histogramm->GetXaxis()->GetXmax();
483
484 histogramm->SetAxisRange(RangeXmin, RangeXmax);
485 Double_t histMean = histogramm->GetMean();
486 Double_t histRMS = histogramm->GetRMS();
487 histogramm->SetAxisRange(hist_start_Xmin, hist_start_Xmax);
488
489
490 //if low events, only RMS and mean;
491 if( integral <= EventsNoFit )
492 {
493 mean = histMean;
494 sigma = histRMS;
495 return 2;
496 }
497
498 TF1 *currfitfunction = NULL;
499 TF1 *bestfitfunction = NULL;
500 Int_t bestOption = -111;
501
502 Double_t bestChisquare = -111., currChisquare = 0;
503
504 for(Int_t opt = 0; !fitoptions[opt].Contains("END"); opt++){
505
506 currfitfunction = new TF1(fit_funk_name, fit_formula.Data(), RangeXmin, RangeXmax);
507
508 histogramm->Fit(currfitfunction, fitoptions[opt], "", RangeXmin, RangeXmax);
509 currChisquare = currfitfunction->GetChisquare();
510
511 if(uVERBOSE >= 3)current_out<<Form("utilites.h/FitHistogramm: current fiting option: \"%s\", chi-square: %.1f",fitoptions[opt].Data(), currChisquare)<<std::endl;
512 if( ( (currChisquare < bestChisquare)&&(currChisquare > 0.0001) )||( (bestChisquare < 0.)&&(currChisquare > 0.0001)) )
513 {
514 bestChisquare = currChisquare;
515 if(bestfitfunction) delete bestfitfunction;
516 bestfitfunction = currfitfunction;
517 bestOption = opt;
518 }
519 else
520 {
521 if(currfitfunction) delete currfitfunction;
522 }
523 }
524
525
526 if( 0. < bestChisquare)
527 {
528 histogramm->Fit(bestfitfunction, fitoptions[bestOption], "", RangeXmin, RangeXmax);
529 if(uVERBOSE >= 2)current_out<<Form("utilites.h/FitHistogramm: best fit [%.2f, %.2f]: %.2f, option %i",
530 RangeXmin, RangeXmax, bestChisquare, bestOption)<<std::endl;
531 histogramm->GetFunction(fit_funk_name)->SetRange(RangeXmin, RangeXmax); //to see selected range on func at hist
532
533
534 mean = bestfitfunction->GetParameter(1);
535 sigma = bestfitfunction->GetParameter(2);
536 }
537 else
538 {
539 if(uVERBOSE >= 1)current_out<<Form("WARING: fitting failure for all options")<<std::endl;
540 mean = histMean;
541 sigma = histRMS;
542 }
543
544 if(bestfitfunction) delete bestfitfunction;
545 return 1;
546
547}
549void ResetArray(TObjArray **array, TString option)
550{
551 Int_t ioption;
552 TObject *currObject = NULL;
553
554 //if option is incorrect
555 if(option.Contains("delete")) ioption = 0;
556 else if(option.Contains("recreate")) ioption = 1;
557 else {if(uVERBOSE >= 1) current_out <<Form("utilites.h/ResetArray::ERROR: array \"%s\" was NOT correctly deleted, wrong option \"%s\"",(*array)->GetName(), option.Data())<<std::endl; return;}
558
559 //if array is empty
560 if(!(*array))
561 {
562 if(ioption == 1){*array = new TObjArray(); return;}
563 if(ioption == 0)return;
564 }
565
566 for(Int_t element = 0; element <= (*array)->GetLast(); element++)
567 {
568 currObject = (*array)->At(element);
569 if(!currObject)continue;
570 if(uVERBOSE >= 4) current_out <<Form("utilites.h/ResetArray::Deleting object \"%s\"",currObject->GetName())<<std::endl;
571 delete currObject;
572 }
573
574 (*array)->Clear();
575
576 delete *array;
577
578 if(ioption == 1){*array = new TObjArray(); return;}
579 if(ioption == 0){*array = NULL; return;}
580
581}
583Int_t ArrayStat(Int_t nElements, Double_t* Array, Double_t &maxElement, Double_t &minElement, Double_t &Amplitude, Double_t &Mean, Double_t &RMS)
584{
585 maxElement = Array[0];
586 minElement = Array[0];
587 Amplitude = Mean = RMS = 0;
588 Double_t currElement = 0, SqSumm = 0;
589
590 for(Int_t i = 0; i < nElements; i++)
591 {
592 currElement = Array[i];
593 if(maxElement < currElement)maxElement = currElement;
594 if(currElement < minElement)minElement = currElement;
595 Mean += currElement;
596 }
597
598 Mean /= (Double_t)nElements;
599 Amplitude = maxElement - minElement;
600
601 for(Int_t i = 0; i < nElements; i++)
602 {
603 currElement = Array[i];
604 SqSumm += (currElement-Mean)*(currElement-Mean);
605 }
606
607 RMS = TMath::Sqrt( SqSumm / (Double_t)nElements );
608
609 return 1;
610}
612TH1 *CreateSmartPojection(TH2 *sourcehist, TString option = "n", const Float_t minbinevX = 1,
613 const Float_t minbinevedge = 0.001, const Int_t rebinX = 1, const Int_t rebinY = 1, Int_t ProjectionRebin = 5, Float_t projWidth = 1.5, TObjArray *monitor_slices_array = NULL)
614{
615 //options: "SIGMA" - draw sigma (RMS) else mean; "FIT" - fit projection, else mean and RMS; "s" - smooth; "n" - no rebin ; "r" - rebin
616
617 if(!sourcehist)
618 {
619 if(uVERBOSE >= 1)cout <<Form("utilites.h/CreateSmartPojection::ERROR: hist pointer is NULL !!!\n")<<endl;
620 return NULL;
621 }
622
623 static const Int_t minBinEvOnEdge = 5;
624 //static const Int_t ProjectionRebin = 5;
625
626 const Bool_t IsDrawSIGMA = option.Contains("SIGMA");
627 const Bool_t IsFitProjection = option.Contains("FIT");
628
629 if(uVERBOSE >= 2)cout <<Form("utilites.h/CreateSmartPojection: Creating %s Projection%s for \"%s\" ... ",
630 (IsDrawSIGMA)? "sigma":"mean", (IsFitProjection)?" by fit":"", sourcehist->GetName())<<endl;
631
632 Int_t NumTotalEv = sourcehist->GetEntries();
633 if(NumTotalEv < 1000)
634 {
635 cout <<Form("utilites.h/CreateSmartPojection::ERROR: Too low events in \"%s\" hist: %i",sourcehist->GetName(), NumTotalEv)<<endl;
636 return NULL;
637 }
638
639 TH2 *rebinedhist = sourcehist->Rebin2D(rebinX, rebinY, Form("%s_rebin",sourcehist->GetName()) );
640
641 if(option.Contains("s"))rebinedhist->Smooth(3);
642
643 Int_t NumBinsY = rebinedhist->GetNbinsY(); //number of bin along Y
644 Int_t NumBinsX = rebinedhist->GetNbinsX(); //number of bin along X
645
646 //average events in bin, does not calculated empty bins
647 Int_t AverageEvInBinX = 0, bincalculated = 0, integral; //events in bin (X axis) on the averages
648 for(Int_t xbin = 1; xbin <= NumBinsX; xbin++)
649 {
650 integral = rebinedhist->Integral(xbin,xbin,1,NumBinsY);
651 if(0 < integral){ AverageEvInBinX += integral; bincalculated++;}
652 }
653 AverageEvInBinX = ceil((Float_t)AverageEvInBinX / (Float_t)bincalculated);
654
655
656 if(uVERBOSE >= 4)cout <<Form("utilites.h/CreateSmartPojection: Average events in bin: %i", AverageEvInBinX)<<endl;
657
658 //fist and last bin with minbinevedge*AverageEvInBinX
659 Int_t firstBinX = -1, lastBinX = -1;
660 for(Int_t xbin = 1; xbin <= NumBinsX; xbin++)
661 {
662 if(firstBinX < 0)
663 if( (ceil(minbinevedge*AverageEvInBinX) <= rebinedhist->Integral(xbin,xbin,1,NumBinsY)) &&
664 ( minBinEvOnEdge <= rebinedhist->Integral(xbin,xbin,1,NumBinsY)) ) firstBinX = xbin;
665
666 if(lastBinX < 0)
667 if( (ceil(minbinevedge*AverageEvInBinX) <= rebinedhist->Integral(NumBinsX-xbin+1,NumBinsX-xbin+1,1,NumBinsY)) &&
668 ( minBinEvOnEdge <= rebinedhist->Integral(NumBinsX-xbin+1,NumBinsX-xbin+1,1,NumBinsY))) lastBinX = NumBinsX-xbin+1;
669
670 if( (lastBinX > 0) && (firstBinX > 0) ) break;
671 }
672
673 if(uVERBOSE >= 4)cout <<Form("utilites.h/CreateSmartPojection: FIRST BIN: %i, %.1f (%.1f ev); LAST BIN: %i, %.1f (%.1f ev)",
674 firstBinX, rebinedhist->GetXaxis()->GetBinCenter(firstBinX), rebinedhist->Integral(firstBinX,firstBinX,1,NumBinsY),
675 lastBinX, rebinedhist->GetXaxis()->GetBinCenter( lastBinX), rebinedhist->Integral(lastBinX, lastBinX, 1,NumBinsY))<<endl;
676 if(option.Contains("r"))if(uVERBOSE >= 4)cout <<Form("Bin before rebinning: %i", lastBinX-firstBinX)<<endl;
677
678 //finding bins: in each bin minimum minQTCCFDbinEv*QTCCFDnumEvInBinXAverage events
679 Float_t *NewBinning = new Float_t[lastBinX-firstBinX+2];
680 Int_t eventsINbin = 0, numberOfBins = 0;
681
682 if(uVERBOSE >= 5)cout <<Form("BINNING: ");
683 NewBinning[0] = rebinedhist->GetXaxis()->GetBinLowEdge(firstBinX);
684 for(Int_t xbin = firstBinX+1; xbin <= lastBinX; xbin++)
685 {
686 eventsINbin += rebinedhist->Integral(xbin,xbin,1,NumBinsY);
687
688 if( (eventsINbin >= AverageEvInBinX*minbinevX) || (xbin == lastBinX) || (!option.Contains("r")))
689 {
690 numberOfBins++;
691 NewBinning[numberOfBins] = rebinedhist->GetXaxis()->GetBinUpEdge(xbin);
692 if(uVERBOSE >= 5)cout <<Form("%.1f ", NewBinning[numberOfBins]);
693 eventsINbin = 0;
694 }
695 } //for sbin
696 if(uVERBOSE >= 4)cout <<""<<endl;
697 if(option.Contains("r"))if(uVERBOSE >= 4)cout <<Form("Bin after rebining: %i", numberOfBins)<<endl;
698
699 if(numberOfBins <= 0)
700 {
701 if(uVERBOSE >= 1)cout <<Form("utilites.h/CreateSmartPojection: Binning failure, number of bins: %i", numberOfBins)<<endl;
702 return NULL;
703 }
704
705 //filling 1D histograms with mean and errors by gauss fit
706 TH1 *ptrProjection = new TH1F(Form("%s_smProj_%s", sourcehist->GetName(),(IsDrawSIGMA)?"sigma":"mean"),"title",numberOfBins,0,1);
707 ptrProjection->GetXaxis()->Set(numberOfBins, NewBinning);
708 ptrProjection->SetOption("E1");
709 ptrProjection->SetTitle(Form("[%s projection%s] %s",
710 (IsDrawSIGMA)? "SIGMA":"MEAN", (IsFitProjection)?" by fit":"", sourcehist->GetTitle()));
711 if(IsFitProjection && IsDrawSIGMA)
712 {
713 ptrProjection->SetMarkerStyle(22);
714
715 }
716
717 Int_t firstbin, lastbin;
718 Double_t mean = 0., sigma = 0., sigmaERROR = 0., meanERROR = 0.;
719 TH1 *ptrProjectionY;
720 static Int_t monitor_slices_number = 0;
721
722 for(Int_t bin = 0; bin < numberOfBins; bin++)
723 {
724 ptrProjectionY = NULL;
725 mean = sigma = sigmaERROR = 0.0001;
726 firstbin = rebinedhist->GetXaxis()->FindBin( NewBinning[bin] );
727 lastbin = rebinedhist->GetXaxis()->FindBin( NewBinning[bin+1] );
728
729 ptrProjectionY = rebinedhist->ProjectionY(Form("%i_%s_slice_bin%i_binmin%i_binmax%i",
730 monitor_slices_number++,(IsDrawSIGMA)? "SIGMA":"MEAN",bin, firstbin,lastbin),firstbin,lastbin);
731 ptrProjectionY->SetTitle( Form("[proj slice: minX %.2f, maxX %.2f] %s", NewBinning[bin], NewBinning[bin+1], ptrProjection->GetTitle()));
732
733 if(ptrProjectionY)
734 {
735 ptrProjectionY->Rebin(ProjectionRebin);
736
737 if(IsFitProjection)
738 {
739 //FitHistogrammInMaxPeak(ptrProjectionY,mean,sigma, projWidth); //width = 1.5
740 FitHistogrammSmartInMaxPeak(ptrProjectionY,mean,sigma); //width = 1.5
741
742 if(ptrProjectionY->GetFunction(fit_funk_name))
743 {
744 meanERROR = ptrProjectionY->GetFunction(fit_funk_name)->GetParError(1);
745 sigmaERROR = ptrProjectionY->GetFunction(fit_funk_name)->GetParError(2);
746 }
747 }
748 else
749 {
750 mean = ptrProjectionY->GetMean();
751 sigma = ptrProjectionY->GetRMS();
752 meanERROR = ptrProjectionY->GetMeanError();
753 sigmaERROR = ptrProjectionY->GetRMSError();
754 }
755
756 if(monitor_slices_array) monitor_slices_array->Add(ptrProjectionY);
757 else delete ptrProjectionY;
758
759 } else { cout <<Form("utilites.h/CreateSmartPojection: WARING: ProjectionY \"%s\" was NOT created", Form("QTCCFDslice_bin_%i_%i",firstbin,lastbin) )<<endl; }
760
761 ptrProjection->SetBinContent(bin+1,(IsDrawSIGMA) ? sigma : mean);
762 ptrProjection->SetBinError(bin+1, (IsDrawSIGMA) ? sigmaERROR*0.5 : meanERROR*0.5);
763
764 }
765
766 if(monitor_slices_array) monitor_slices_array->Add(rebinedhist);
767 else delete rebinedhist;
768
769 delete[] NewBinning;
770 return ptrProjection;
771}
773TObjArray* HDraw(const TObjArray *Histogramms_arr, TString name = "",
774 int Hr = 4, int Vt = 4, const int Pair = 1,
775 int min = 0, int max = 100, int Log = 0, int Leg = 0, const TString param = "", int Norm = 0)
776{
777 gStyle->SetPalette(1);
778
779 //name, HxV,number paired hits, min, max,log (1 - y, 2 - z), legend (0-no, 1-all, 2-first pad)
780 int H_Last = Histogramms_arr->GetLast();
781 if(H_Last <= 0) return NULL;
782
783 if(min < 0) min = 0;
784 if(min > H_Last) min = H_Last;
785 if(max < 0) max = 0;
786 if(max > H_Last) max = H_Last;
787 if(Hr < 1) Hr = 1;
788 if(Vt < 1) Vt = 1;
789 if((Leg == 2)&&(Hr*Vt == 1))Leg = 0;
790
791 if(name == "") name = Histogramms_arr->GetName();
792
793 TObjArray *Canvas_arr = new TObjArray();
794 Canvas_arr->SetOwner();
795
796 int flagD = Hr*Vt;
797
798 for(int t=min; t <= max; t++){
799 //((TH1F*)Histogramms_arr->At(t)) -> SetLineWidth(3);
800 // ((TH1F*)Histogramms_arr->At(t)) -> SetLineColor(t%Pair +1);
801 //((TH1F*)Histogramms_arr->At(t)) -> SetFillColor(t%Pair +1);
802 //((TH1F*)Histogramms_arr->At(t)) -> SetFillStyle(3001);
803 //((TH1F*)Histogramms_arr->At(t)) -> SetMarkerColor(t%Pair +1);
804 }
805
806 for(int t=min; t <= max; t++){
807
808 if(++flagD > Hr*Vt){
809 char n_Canv[30]; sprintf(n_Canv, "%s#%d",name.Data(),Canvas_arr->GetLast()+1);
810 Canvas_arr->Add( new TCanvas(n_Canv, n_Canv, 2000, 2000) );
811 ((TCanvas*)Canvas_arr->Last())->Divide(Hr,Vt,0.0001,0.0001);
812 flagD = 1;}
813
814 int t_0 = t;
815 ((TCanvas*)Canvas_arr->Last())->cd(flagD);
816 //if(Histogramms_arr->At(t) != NULL)
817 if(!Norm){ ((TH1F*)Histogramms_arr->At(t))->Draw(param.Data()); }
818 else
819 //if(Histogramms_arr->At(t) != NULL)
820 ((TH1F*)Histogramms_arr->At(t))->DrawNormalized(param.Data(), 1);
821 for(int c_Pair = 1; c_Pair < Pair; c_Pair++){
822 t++;
823 //if(Histogramms_arr->At(t) != NULL)
824 if(!Norm){ ((TH1F*)Histogramms_arr->At(t))->Draw("sames");}
825 else
826 //if(Histogramms_arr->At(t) != NULL)
827 ((TH1F*)Histogramms_arr->At(t))->DrawNormalized("sames", 1);
828 }
829
830 //to build big legend on 1st pad drawing hist set twice
831 if((Leg == 2)&&(flagD == 1)){
832 gPad -> BuildLegend(0,0,1,1);
833 flagD++;
834 t = t_0;
835 ((TCanvas*)Canvas_arr->Last())->cd(flagD);
836 //if(Histogramms_arr->At(t) != NULL)
837 if(!Norm) ((TH1F*)Histogramms_arr->At(t))->Draw(param.Data());
838 else
839 //if(Histogramms_arr->At(t) != NULL)
840 ((TH1F*)Histogramms_arr->At(t))->DrawNormalized(param.Data(),1);
841 for(int c_Pair = 1; c_Pair < Pair; c_Pair++){
842 t++;
843 //if(Histogramms_arr->At(t) != NULL)
844 if(!Norm) ((TH1F*)Histogramms_arr->At(t))->Draw("sames");
845 else
846 //if(Histogramms_arr->At(t) != NULL)
847 ((TH1F*)Histogramms_arr->At(t))->DrawNormalized("sames", 1);
848 }
849 }
850
851
852 if(Log == 1) gPad -> SetLogy(1);
853 if(Log == 2) gPad -> SetLogz(1);
854
855 //if(Leg == 1) ;//gPad -> BuildLegend();
856
857 }
858
859 return Canvas_arr;
860}
862#endif // UTILITES_H
Sigma
Double_t Sigma(vector< Double_t > dist, vector< Double_t > w)
Definition BmnMath.cxx:890
Draw
int Draw(TCanvas *c, vector< TH1 * > hv, vector< string > optV)
Definition BmnMonProf.cxx:332
i
int i
Definition P4_F32vec4.h:22
min
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:30
max
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition P4_F32vec4.h:31
SaveArrayOfObjectsArraysInFile
Int_t SaveArrayOfObjectsArraysInFile(TObjArray *ObjectsArray=0, const TString filename="data", const TString path="./", const Int_t deleteArrays=0)
Definition utilites.h:242
SaveArrayStructInFile
Int_t SaveArrayStructInFile(TCollection *array, TFile *file, TString start_path="")
Definition utilites.h:198
current_out
std::ostream & current_out
Definition utilites.h:45
ArrayStat
Int_t ArrayStat(const Int_t nElements, const Double_t *Array, Double_t &maxElement, Double_t &minElement, Double_t &Amplitude, Double_t &Mean, Double_t &RMS)
uVERBOSE
#define uVERBOSE
Definition utilites.h:48
CollectRAWfilesInChain
Int_t CollectRAWfilesInChain(TChain *rawTree, const TString dirname, const Int_t maxNfiles=50)
Definition utilites.h:128
FilesListFromDir
Int_t FilesListFromDir(const TString filespath, TObjArray *const collection, Int_t max_num, TString content)
Definition utilites.h:163
CreateXprofile
TObject * CreateXprofile(const TH1 *sourcehist)
Definition utilites.h:118
RUNnumberInString
Int_t RUNnumberInString(const TString &str)
Definition utilites.h:101
FitHistogrammSmartInMaxPeak
Int_t FitHistogrammSmartInMaxPeak(TH1 *histogramm, Double_t &mean, Double_t &sigma, Double_t nRMSfit=2., Double_t RangeXmin=0., Double_t RangeXmax=0.)
Definition utilites.h:395
ResetArray
void ResetArray(TObjArray **array, const TString option)
Definition utilites.h:549
ConstructFilesListFromDirectory
Int_t ConstructFilesListFromDirectory(const TString &Directory, TCollection *const collection)
Definition utilites.h:67
FitHistogramm
Int_t FitHistogramm(TH1 *histogramm, Double_t &mean, Double_t &sigma, Double_t RangeXmin=0., Double_t RangeXmax=0., TString fit_formula="gaus")
Definition utilites.h:451
FitHistogrammInMaxPeak
Int_t FitHistogrammInMaxPeak(TH1 *histogramm, Double_t &mean, Double_t &sigma, Double_t PeakWidth, Double_t RangeXmin=0., Double_t RangeXmax=0.)
Definition utilites.h:361
HDraw
TObjArray * HDraw(const TObjArray *Histogramms_arr, TString name="", int Hr=4, int Vt=4, const int Pair=1, int min=0, int max=100, int Log=0, int Leg=0, const TString param="", int Norm=0)
Definition utilites.h:773
FitHistogrammInFirstMaxPeak
Int_t FitHistogrammInFirstMaxPeak(TH1 *histogramm, Double_t &Mean, Double_t &Sigma, Double_t peaksRangeMin, Double_t peaksRangeMax, Double_t peakRange, Double_t peakRMSRange=1.5)
Definition utilites.h:313
FitHistogrammInRMSRange
Int_t FitHistogrammInRMSRange(TH1 *histogramm, Double_t &mean, Double_t &sigma, Double_t RMSrange, Double_t RangeXmin=0, Double_t RangeXmax=0)
Definition utilites.h:333
CreateSmartPojection
TH1 * CreateSmartPojection(TH2 *sourcehist, TString option="n", const Float_t minbinevX=1, const Float_t minbinevedge=0.001, const Int_t rebinX=1, const Int_t rebinY=1, Int_t ProjectionRebin=5, Float_t projWidth=1.5, TObjArray *monitor_slices_array=NULL)
Definition utilites.h:612
fit_funk_name
#define fit_funk_name
Definition utilites.h:49
std
STL namespace.