[u/mrichter/AliRoot.git] / T0 / MakeTrendT0.C

#define NPMTs 24

int MakeTrendT0( char *infile, int run, char* ocdbStorage="raw://") {

  gSystem->Load("libANALYSIS");
  gSystem->Load("libANALYSISalice");
  gSystem->Load("libCORRFW");
  gSystem->Load("libTENDER");
  gSystem->Load("libPWGPP.so");

  char *outfile = "trending.root";
 
  if(!infile) return -1;
  if(!outfile) return -1;

  TFile *f = TFile::Open(infile,"read");
  if (!f) {
    printf("File %s not available\n", infile);
    return -1;
  }

  //            LOAD HISTOGRAMS FROM QAresults.root   
  TObjArray   *fTzeroObject = (TObjArray*) f->Get("T0_Performance/QAT0chists"); 

  TH1F* fTzeroORAplusORC =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORAplusORC"))->Clone("A"); 
  TH1F* fResolution      =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fResolution"))->Clone("B");
  TH1F* fTzeroORA        =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORA"))->Clone("C");
  TH1F* fTzeroORC        =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORC"))->Clone("D");

  TH2F *fTimeVSAmplitude[NPMTs];//counting PMTs from 0
  TH1D *fAmplitude[NPMTs];
  TH1D *fTime[NPMTs];
  //-----> add new histogram here  


  // sigma fit of resolution, mean T0A, T0C and T0C&A,  mean amplitude for each PMT
  double resolutionSigma = -9999; //dummy init
  double tzeroOrAPlusOrC = -9999; //dummy init
  double tzeroOrA        = -9999; //dummy init
  double tzeroOrC        = -9999; //dummy init
  double meanAmplitude[NPMTs]; //dummy init
  double meanTime[NPMTs]; //dummy init
  double timeDelayOCDB[NPMTs]; //dummy init
 
  //.................................................................................

  for(int ipmt=1; ipmt<=NPMTs; ipmt++){ //loop over all PMTs
    fTimeVSAmplitude[ipmt-1] =(TH2F*) ((TH2F*) fTzeroObject->FindObject(Form("fTimeVSAmplitude%d",ipmt)))->Clone(Form("E%d",ipmt));
    int nbinsX = fTimeVSAmplitude[ipmt-1]->GetNbinsX();
    int nbinsY = fTimeVSAmplitude[ipmt-1]->GetNbinsY();
    //Amplitude
    fAmplitude[ipmt-1] = (TH1D*) fTimeVSAmplitude[ipmt-1]->ProjectionX(Form("fAmplitude%d",ipmt), 1, nbinsY); 
    meanAmplitude[ipmt-1] = -9999; //dummy init 
    if(fAmplitude[ipmt-1]->GetEntries()>0){
      meanAmplitude[ipmt-1] = fAmplitude[ipmt-1]->GetMean();
    }

    //Time
    fTime[ipmt-1] = (TH1D*) fTimeVSAmplitude[ipmt-1]->ProjectionY(Form("fTime%d",ipmt), 1, nbinsX); 
    meanTime[ipmt-1] = -9999; //dummy init 
    if(fTime[ipmt-1]->GetEntries()>0){
      if(fTime[ipmt-1]->GetEntries()>20){
        meanTime[ipmt-1] = GetParameterGaus((TH1F*) fTime[ipmt-1], 1); // Mean Time
      }else if(fTime[ipmt-1]->GetEntries()>0){
        meanTime[ipmt-1] = fTime[ipmt-1]->GetMean();
      } 
    }
  }


  if(fResolution->GetEntries()>20){
    resolutionSigma = GetParameterGaus(fResolution, 2); //gaussian sigma 
  }else if(fResolution->GetEntries()>0){
    resolutionSigma = fResolution->GetRMS(); //gaussian sigma 
  }

  if(fTzeroORAplusORC->GetEntries()>20){
    tzeroOrAPlusOrC = GetParameterGaus(fTzeroORAplusORC, 1); //gaussian mean 
  }else if(fTzeroORAplusORC->GetEntries()>0){
    tzeroOrAPlusOrC = fTzeroORAplusORC->GetMean();
  }

  if(fTzeroORA->GetEntries()>20){ 
    tzeroOrA = GetParameterGaus(fTzeroORA, 1); //gaussian mean
  }else if(fTzeroORA->GetEntries()>0){
    tzeroOrA = fTzeroORA->GetMean();
  }

  if(fTzeroORC->GetEntries()>20){
    tzeroOrC = GetParameterGaus(fTzeroORC, 1); //gaussian mean 
  }else if(fTzeroORC->GetEntries()>0){
    tzeroOrC = fTzeroORC->GetMean(); //gaussian mean 
  }

  //-----> analyze the new histogram here and set mean/sigma


  f->Close();

  //-------------------- READ OCDB TIME DELAYS ---------------------------
  // Arguments:
  AliCDBManager* man = AliCDBManager::Instance();
  man->SetDefaultStorage(ocdbStorage);
  man->SetRun(run);
  AliCDBEntry *entry = AliCDBManager::Instance()->Get("T0/Calib/TimeDelay");
  AliT0CalibTimeEq *clb = (AliT0CalibTimeEq*)entry->GetObject(); 
  for (Int_t i=0; i<NPMTs; i++){
    timeDelayOCDB[i] = 0;
    if(clb)
      timeDelayOCDB[i] = clb->GetCFDvalue(i,0);
  }


  //--------------- write walues to the output ------------ 
  TTreeSRedirector* pcstream = NULL;
  pcstream = new TTreeSRedirector(outfile);
  if (!pcstream) return;


  TFile *x =  pcstream->GetFile();
  x->cd();

  TObjString runType;
  Int_t startTimeGRP=0;
  Int_t stopTimeGRP=0;
  Int_t time=0;
  Int_t duration=0;

  time     = (startTimeGRP+stopTimeGRP)/2;
  duration = (stopTimeGRP-startTimeGRP);

  (*pcstream)<<"t0QA"<<
                  "run="<<run;//<<
                 //"time="<<time<<
            // "startTimeGRP="<<startTimeGRP<<
            //  "stopTimeGRP="<<stopTimeGRP<<
              ///    "duration="<<duration<<
               //   "runType.="<<&runType;
               

 (*pcstream)<<"t0QA"<<
             "resolution="<< resolutionSigma<<
             "tzeroOrAPlusOrC="<< tzeroOrAPlusOrC<<
             "tzeroOrA="<< tzeroOrA<<
             "tzeroOrC="<< tzeroOrC<<
             "amplPMT1="<<meanAmplitude[0]<<
             "amplPMT2="<<meanAmplitude[1]<<
             "amplPMT3="<<meanAmplitude[2]<<
             "amplPMT4="<<meanAmplitude[3]<<
             "amplPMT5="<<meanAmplitude[4]<<
             "amplPMT6="<<meanAmplitude[5]<<
             "amplPMT7="<<meanAmplitude[6]<<
             "amplPMT8="<<meanAmplitude[7]<<
             "amplPMT9="<<meanAmplitude[8]<<
             "amplPMT10="<<meanAmplitude[9]<<
             "amplPMT11="<<meanAmplitude[10]<<
             "amplPMT12="<<meanAmplitude[11]<<
             "amplPMT13="<<meanAmplitude[12]<<
             "amplPMT14="<<meanAmplitude[13]<<
             "amplPMT15="<<meanAmplitude[14]<<
             "amplPMT16="<<meanAmplitude[15]<<
             "amplPMT17="<<meanAmplitude[16]<<
             "amplPMT18="<<meanAmplitude[17]<<
             "amplPMT19="<<meanAmplitude[18]<<
             "amplPMT20="<<meanAmplitude[19]<<
             "amplPMT21="<<meanAmplitude[20]<<
             "amplPMT22="<<meanAmplitude[21]<<
             "amplPMT23="<<meanAmplitude[22]<<
             "amplPMT24="<<meanAmplitude[23]<<
             "timePMT1="<<meanTime[0]<<
             "timePMT2="<<meanTime[1]<<
             "timePMT3="<<meanTime[2]<<
             "timePMT4="<<meanTime[3]<<
             "timePMT5="<<meanTime[4]<<
             "timePMT6="<<meanTime[5]<<
             "timePMT7="<<meanTime[6]<<
             "timePMT8="<<meanTime[7]<<
             "timePMT9="<<meanTime[8]<<
             "timePMT10="<<meanTime[9]<<
             "timePMT11="<<meanTime[10]<<
             "timePMT12="<<meanTime[11]<<
             "timePMT13="<<meanTime[12]<<
             "timePMT14="<<meanTime[13]<<
             "timePMT15="<<meanTime[14]<<
             "timePMT16="<<meanTime[15]<<
             "timePMT17="<<meanTime[16]<<
             "timePMT18="<<meanTime[17]<<
             "timePMT19="<<meanTime[18]<<
             "timePMT20="<<meanTime[19]<<
             "timePMT21="<<meanTime[20]<<
             "timePMT22="<<meanTime[21]<<
             "timePMT23="<<meanTime[22]<<
             "timePMT24="<<meanTime[23];
(*pcstream)<<"t0QA"<<
             "timeDelayPMT1="<<timeDelayOCDB[0]<<
             "timeDelayPMT2="<<timeDelayOCDB[1]<<
             "timeDelayPMT3="<<timeDelayOCDB[2]<<
             "timeDelayPMT4="<<timeDelayOCDB[3]<<
             "timeDelayPMT5="<<timeDelayOCDB[4]<<
             "timeDelayPMT6="<<timeDelayOCDB[5]<<
             "timeDelayPMT7="<<timeDelayOCDB[6]<<
             "timeDelayPMT8="<<timeDelayOCDB[7]<<
             "timeDelayPMT9="<<timeDelayOCDB[8]<<
             "timeDelayPMT10="<<timeDelayOCDB[9]<<
             "timeDelayPMT11="<<timeDelayOCDB[10]<<
             "timeDelayPMT12="<<timeDelayOCDB[11]<<
             "timeDelayPMT13="<<timeDelayOCDB[12]<<
             "timeDelayPMT14="<<timeDelayOCDB[13]<<
             "timeDelayPMT15="<<timeDelayOCDB[14]<<
             "timeDelayPMT16="<<timeDelayOCDB[15]<<
             "timeDelayPMT17="<<timeDelayOCDB[16]<<
             "timeDelayPMT18="<<timeDelayOCDB[17]<<
             "timeDelayPMT19="<<timeDelayOCDB[18]<<
             "timeDelayPMT20="<<timeDelayOCDB[19]<<
             "timeDelayPMT21="<<timeDelayOCDB[20]<<
             "timeDelayPMT22="<<timeDelayOCDB[21]<<
             "timeDelayPMT23="<<timeDelayOCDB[22]<<
             "timeDelayPMT24="<<timeDelayOCDB[23];

            //-----> add the mean/sigma of the new histogram here      
 
 (*pcstream)<<"t0QA"<<"\n";
  
  pcstream->Close(); 
 
  delete pcstream; 

  return;

}

//_____________________________________________________________________________
double GetParameterGaus(TH1F *histo, int whichParameter){

   int    maxBin   =  histo->GetMaximumBin(); 
   double max      =  (histo->GetBinContent(maxBin-1) + histo->GetBinContent(maxBin) + histo->GetBinContent(maxBin+1))/3;
   double mean     =  histo->GetBinCenter(maxBin); //mean
   double lowfwhm  =  histo->GetBinCenter(histo->FindFirstBinAbove(max/2));
   double highfwhm =  histo->GetBinCenter(histo->FindLastBinAbove(max/2));
   double sigma    = (highfwhm - lowfwhm)/2.35482; //estimate fwhm  FWHM = 2.35482*sigma

   TF1 *gaussfit   = new TF1("gaussfit","gaus", mean - 4*sigma, mean + 4*sigma); // fit in +- 4 sigma window
   gaussfit->SetParameters(max, mean, sigma); 

   if(whichParameter==2)
     histo->Fit(gaussfit,"RQNI");
   else
     histo->Fit(gaussfit,"RQN");

   double parValue = gaussfit->GetParameter(whichParameter); 

   delete gaussfit; 
 
   return parValue;
}
Commit	Line	Data
b1a2430c	1	#define NPMTs 24
b1a2430c	2
302db621	3	int MakeTrendT0( char infile, int run, char ocdbStorage="raw://") {
b1a2430c	4
427f2633	5	gSystem->Load("libANALYSIS");
	6	gSystem->Load("libANALYSISalice");
	7	gSystem->Load("libCORRFW");
	8	gSystem->Load("libTENDER");
	9	gSystem->Load("libPWGPP.so");
	10
b1a2430c	11	char *outfile = "trending.root";
	12
	13	if(!infile) return -1;
	14	if(!outfile) return -1;
	15
	16	TFile *f = TFile::Open(infile,"read");
	17	if (!f) {
	18	printf("File %s not available\n", infile);
	19	return -1;
	20	}
	21
	22	// LOAD HISTOGRAMS FROM QAresults.root
	23	TObjArray fTzeroObject = (TObjArray) f->Get("T0_Performance/QAT0chists");
	24
	25	TH1F* fTzeroORAplusORC =(TH1F) ((TH1F) fTzeroObject->FindObject("fTzeroORAplusORC"))->Clone("A");
	26	TH1F* fResolution =(TH1F) ((TH1F) fTzeroObject->FindObject("fResolution"))->Clone("B");
	27	TH1F* fTzeroORA =(TH1F) ((TH1F) fTzeroObject->FindObject("fTzeroORA"))->Clone("C");
	28	TH1F* fTzeroORC =(TH1F) ((TH1F) fTzeroObject->FindObject("fTzeroORC"))->Clone("D");
	29
	30	TH2F *fTimeVSAmplitude[NPMTs];//counting PMTs from 0
	31	TH1D *fAmplitude[NPMTs];
	32	TH1D *fTime[NPMTs];
	33	//-----> add new histogram here
	34
	35
	36
	37	// sigma fit of resolution, mean T0A, T0C and T0C&A, mean amplitude for each PMT
	38	double resolutionSigma = -9999; //dummy init
	39	double tzeroOrAPlusOrC = -9999; //dummy init
	40	double tzeroOrA = -9999; //dummy init
	41	double tzeroOrC = -9999; //dummy init
	42	double meanAmplitude[NPMTs]; //dummy init
	43	double meanTime[NPMTs]; //dummy init
	44	double timeDelayOCDB[NPMTs]; //dummy init
	45
	46	//.................................................................................
	47
	48	for(int ipmt=1; ipmt<=NPMTs; ipmt++){ //loop over all PMTs
	49	fTimeVSAmplitude[ipmt-1] =(TH2F) ((TH2F) fTzeroObject->FindObject(Form("fTimeVSAmplitude%d",ipmt)))->Clone(Form("E%d",ipmt));
	50	int nbinsX = fTimeVSAmplitude[ipmt-1]->GetNbinsX();
	51	int nbinsY = fTimeVSAmplitude[ipmt-1]->GetNbinsY();
	52	//Amplitude
	53	fAmplitude[ipmt-1] = (TH1D*) fTimeVSAmplitude[ipmt-1]->ProjectionX(Form("fAmplitude%d",ipmt), 1, nbinsY);
	54	meanAmplitude[ipmt-1] = -9999; //dummy init
	55	if(fAmplitude[ipmt-1]->GetEntries()>0){
	56	meanAmplitude[ipmt-1] = fAmplitude[ipmt-1]->GetMean();
	57	}
	58
	59	//Time
	60	fTime[ipmt-1] = (TH1D*) fTimeVSAmplitude[ipmt-1]->ProjectionY(Form("fTime%d",ipmt), 1, nbinsX);
	61	meanTime[ipmt-1] = -9999; //dummy init
	62	if(fTime[ipmt-1]->GetEntries()>0){
	63	if(fTime[ipmt-1]->GetEntries()>20){
	64	meanTime[ipmt-1] = GetParameterGaus((TH1F*) fTime[ipmt-1], 1); // Mean Time
	65	}else if(fTime[ipmt-1]->GetEntries()>0){
	66	meanTime[ipmt-1] = fTime[ipmt-1]->GetMean();
	67	}
	68	}
	69	}
	70
	71
	72
	73	if(fResolution->GetEntries()>20){
	74	resolutionSigma = GetParameterGaus(fResolution, 2); //gaussian sigma
75	}else if(fResolution->GetEntries()>0){
76	resolutionSigma = fResolution->GetRMS(); //gaussian sigma
77	}
78
79	if(fTzeroORAplusORC->GetEntries()>20){
80	tzeroOrAPlusOrC = GetParameterGaus(fTzeroORAplusORC, 1); //gaussian mean
81	}else if(fTzeroORAplusORC->GetEntries()>0){
82	tzeroOrAPlusOrC = fTzeroORAplusORC->GetMean();
83	}
84
85	if(fTzeroORA->GetEntries()>20){
86	tzeroOrA = GetParameterGaus(fTzeroORA, 1); //gaussian mean
87	}else if(fTzeroORA->GetEntries()>0){
88	tzeroOrA = fTzeroORA->GetMean();
89	}
90
91	if(fTzeroORC->GetEntries()>20){
92	tzeroOrC = GetParameterGaus(fTzeroORC, 1); //gaussian mean
93	}else if(fTzeroORC->GetEntries()>0){
94	tzeroOrC = fTzeroORC->GetMean(); //gaussian mean
95	}
96
97	//-----> analyze the new histogram here and set mean/sigma
98
99
100	f->Close();
101
102	//-------------------- READ OCDB TIME DELAYS ---------------------------
103	// Arguments:
104	AliCDBManager* man = AliCDBManager::Instance();
302db621	105	man->SetDefaultStorage(ocdbStorage);
b1a2430c	106	man->SetRun(run);
	107	AliCDBEntry *entry = AliCDBManager::Instance()->Get("T0/Calib/TimeDelay");
	108	AliT0CalibTimeEq clb = (AliT0CalibTimeEq)entry->GetObject();
	109	for (Int_t i=0; i<NPMTs; i++){
	110	timeDelayOCDB[i] = 0;
	111	if(clb)
	112	timeDelayOCDB[i] = clb->GetCFDvalue(i,0);
	113	}
	114
	115
	116
	117	//--------------- write walues to the output ------------
	118	TTreeSRedirector* pcstream = NULL;
	119	pcstream = new TTreeSRedirector(outfile);
	120	if (!pcstream) return;
	121
	122
	123	TFile *x = pcstream->GetFile();
	124	x->cd();
	125
	126	TObjString runType;
	127	Int_t startTimeGRP=0;
	128	Int_t stopTimeGRP=0;
	129	Int_t time=0;
	130	Int_t duration=0;
	131
	132	time = (startTimeGRP+stopTimeGRP)/2;
	133	duration = (stopTimeGRP-startTimeGRP);
	134
	135	(*pcstream)<<"t0QA"<<
	136	"run="<<run;//<<
	137	//"time="<<time<<
	138	// "startTimeGRP="<<startTimeGRP<<
	139	// "stopTimeGRP="<<stopTimeGRP<<
	140	/// "duration="<<duration<<
	141	// "runType.="<<&runType;
	142
	143
	144	(*pcstream)<<"t0QA"<<
	145	"resolution="<< resolutionSigma<<
	146	"tzeroOrAPlusOrC="<< tzeroOrAPlusOrC<<
	147	"tzeroOrA="<< tzeroOrA<<
	148	"tzeroOrC="<< tzeroOrC<<
	149	"amplPMT1="<<meanAmplitude[0]<<
	150	"amplPMT2="<<meanAmplitude[1]<<
	151	"amplPMT3="<<meanAmplitude[2]<<
	152	"amplPMT4="<<meanAmplitude[3]<<
	153	"amplPMT5="<<meanAmplitude[4]<<
	154	"amplPMT6="<<meanAmplitude[5]<<
	155	"amplPMT7="<<meanAmplitude[6]<<
	156	"amplPMT8="<<meanAmplitude[7]<<
	157	"amplPMT9="<<meanAmplitude[8]<<
	158	"amplPMT10="<<meanAmplitude[9]<<
	159	"amplPMT11="<<meanAmplitude[10]<<
	160	"amplPMT12="<<meanAmplitude[11]<<
	161	"amplPMT13="<<meanAmplitude[12]<<
	162	"amplPMT14="<<meanAmplitude[13]<<
	163	"amplPMT15="<<meanAmplitude[14]<<
	164	"amplPMT16="<<meanAmplitude[15]<<
	165	"amplPMT17="<<meanAmplitude[16]<<
	166	"amplPMT18="<<meanAmplitude[17]<<
	167	"amplPMT19="<<meanAmplitude[18]<<
	168	"amplPMT20="<<meanAmplitude[19]<<
	169	"amplPMT21="<<meanAmplitude[20]<<
170	"amplPMT22="<<meanAmplitude[21]<<
171	"amplPMT23="<<meanAmplitude[22]<<
172	"amplPMT24="<<meanAmplitude[23]<<
173	"timePMT1="<<meanTime[0]<<
174	"timePMT2="<<meanTime[1]<<
175	"timePMT3="<<meanTime[2]<<
176	"timePMT4="<<meanTime[3]<<
177	"timePMT5="<<meanTime[4]<<
178	"timePMT6="<<meanTime[5]<<
179	"timePMT7="<<meanTime[6]<<
180	"timePMT8="<<meanTime[7]<<
181	"timePMT9="<<meanTime[8]<<
182	"timePMT10="<<meanTime[9]<<
183	"timePMT11="<<meanTime[10]<<
184	"timePMT12="<<meanTime[11]<<
185	"timePMT13="<<meanTime[12]<<
186	"timePMT14="<<meanTime[13]<<
187	"timePMT15="<<meanTime[14]<<
188	"timePMT16="<<meanTime[15]<<
189	"timePMT17="<<meanTime[16]<<
190	"timePMT18="<<meanTime[17]<<
191	"timePMT19="<<meanTime[18]<<
192	"timePMT20="<<meanTime[19]<<
193	"timePMT21="<<meanTime[20]<<
194	"timePMT22="<<meanTime[21]<<
195	"timePMT23="<<meanTime[22]<<
196	"timePMT24="<<meanTime[23];
197	(*pcstream)<<"t0QA"<<
198	"timeDelayPMT1="<<timeDelayOCDB[0]<<
199	"timeDelayPMT2="<<timeDelayOCDB[1]<<
200	"timeDelayPMT3="<<timeDelayOCDB[2]<<
201	"timeDelayPMT4="<<timeDelayOCDB[3]<<
202	"timeDelayPMT5="<<timeDelayOCDB[4]<<
203	"timeDelayPMT6="<<timeDelayOCDB[5]<<
204	"timeDelayPMT7="<<timeDelayOCDB[6]<<
205	"timeDelayPMT8="<<timeDelayOCDB[7]<<
206	"timeDelayPMT9="<<timeDelayOCDB[8]<<
207	"timeDelayPMT10="<<timeDelayOCDB[9]<<
208	"timeDelayPMT11="<<timeDelayOCDB[10]<<
209	"timeDelayPMT12="<<timeDelayOCDB[11]<<
210	"timeDelayPMT13="<<timeDelayOCDB[12]<<
211	"timeDelayPMT14="<<timeDelayOCDB[13]<<
212	"timeDelayPMT15="<<timeDelayOCDB[14]<<
213	"timeDelayPMT16="<<timeDelayOCDB[15]<<
214	"timeDelayPMT17="<<timeDelayOCDB[16]<<
215	"timeDelayPMT18="<<timeDelayOCDB[17]<<
216	"timeDelayPMT19="<<timeDelayOCDB[18]<<
217	"timeDelayPMT20="<<timeDelayOCDB[19]<<
218	"timeDelayPMT21="<<timeDelayOCDB[20]<<
219	"timeDelayPMT22="<<timeDelayOCDB[21]<<
220	"timeDelayPMT23="<<timeDelayOCDB[22]<<
221	"timeDelayPMT24="<<timeDelayOCDB[23];
222
223	//-----> add the mean/sigma of the new histogram here
224
225	(*pcstream)<<"t0QA"<<"\n";
226
227	pcstream->Close();
228
229	delete pcstream;
230
231	return;
232
233	}
234
235	//_____________________________________________________________________________
236	double GetParameterGaus(TH1F *histo, int whichParameter){
237
238	int maxBin = histo->GetMaximumBin();
239	double max = (histo->GetBinContent(maxBin-1) + histo->GetBinContent(maxBin) + histo->GetBinContent(maxBin+1))/3;
240	double mean = histo->GetBinCenter(maxBin); //mean
241	double lowfwhm = histo->GetBinCenter(histo->FindFirstBinAbove(max/2));
242	double highfwhm = histo->GetBinCenter(histo->FindLastBinAbove(max/2));
243	double sigma = (highfwhm - lowfwhm)/2.35482; //estimate fwhm FWHM = 2.35482*sigma
244
245	TF1 gaussfit = new TF1("gaussfit","gaus", mean - 4sigma, mean + 4*sigma); // fit in +- 4 sigma window
246	gaussfit->SetParameters(max, mean, sigma);
247
248	if(whichParameter==2)
249	histo->Fit(gaussfit,"RQNI");
250	else
251	histo->Fit(gaussfit,"RQN");
252
253	double parValue = gaussfit->GetParameter(whichParameter);
254
255	delete gaussfit;
256
257	return parValue;
258	}
259
260