X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=T0%2FAliT0CalibWalk.cxx;h=59007069f5302a3b39b0155b2eb26deff40d64fd;hb=841039b29d8dd211c987b0061fb7a13ad578fe8b;hp=130b3ca35f2a5d3271049db7d8b7871086240276;hpb=700e6b362299f3fd6dc29cb8b773e032e9821a7c;p=u%2Fmrichter%2FAliRoot.git

diff --git a/T0/AliT0CalibWalk.cxx b/T0/AliT0CalibWalk.cxx
index 130b3ca35f2..59007069f53 100644
--- a/T0/AliT0CalibWalk.cxx
+++ b/T0/AliT0CalibWalk.cxx
@@ -12,7 +12,7 @@
  * about the suitability of this software for any purpose. It is          *
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
-
+ 
 /* $Id$ */
 
 ///////////////////////////////////////////////////////////////////////////////
@@ -33,14 +33,18 @@
 #include <Riostream.h>
 #include <TSpectrum.h>
 #include <TProfile.h>
+#include <TF1.h>
 
-
+using std::cout;
+using std::endl;
 ClassImp(AliT0CalibWalk)
 
 //________________________________________________________________
   AliT0CalibWalk::AliT0CalibWalk():   TNamed(),
 				      fWalk(0),
-				      fAmpLEDRec(0)
+				      fAmpLEDRec(0),
+				      fQTC(0),
+				      fAmpLED(0)
 {
   //
 }
@@ -48,7 +52,9 @@ ClassImp(AliT0CalibWalk)
 //________________________________________________________________
 AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
 				      fWalk(0),
-				      fAmpLEDRec(0)				      
+				      fAmpLEDRec(0),				  				      fQTC(0),
+				      fAmpLED(0)
+    
 {
   TString namst = "Calib_";
   namst += name;
@@ -61,8 +67,10 @@ AliT0CalibWalk::AliT0CalibWalk(const char* name):TNamed(),
 AliT0CalibWalk::AliT0CalibWalk(const AliT0CalibWalk& calibda) :
   TNamed(calibda),		
   fWalk(0),
-  fAmpLEDRec(0)
-
+  fAmpLEDRec(0),
+  fQTC(0),
+  fAmpLED(0)
+ 
 {
 // copy constructor
   SetName(calibda.GetName());
@@ -86,215 +94,167 @@ AliT0CalibWalk::~AliT0CalibWalk()
 {
   //
 }
-//________________________________________________________________
-void AliT0CalibWalk::SetWalk(Int_t ipmt)
-{
-  //read QTC walk graph from external file
+//________________________________________________________________  
 
-  Int_t mv, ps; 
-  Int_t x[70000], y[70000], index[70000];
-  Float_t time[10000],amplitude[10000];
-  string buffer;
-  Bool_t down=false;
-  
-  const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-Amp.txt");
-  ifstream inFile(filename);
-  if(!inFile) AliError(Form("Cannot open file %s !",filename));
+Bool_t AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
+{
+  //make walk corerction for preprocessor
+  Float_t sigma,cfdmean, qtmean, ledmean;
+  Bool_t ok=true;
+  Float_t   mips[50];
   
-  Int_t i=0;
-  while(getline(inFile,buffer)){
-    inFile >> ps >> mv;
-
-    x[i]=ps; y[i]=mv;
-    i++;
+  gFile = TFile::Open(laserFile);
+  if(!gFile) {
+    AliError("No input laser data found ");
   }
-  inFile.close();
-  cout<<" number of data "<<i<<endl;
- 
-  TMath::Sort(i, y, index,down);
-  Int_t amp=0, iin=0, isum=0, sum=0;
-  Int_t ind=0;
-  for (Int_t ii=0; ii<i; ii++)
+  else
     {
-      ind=index[ii];
-      if(y[ind] == amp)
-	{
-	  sum +=x[ind];
-	  iin++;
-	}
-      else
+      //      gFile->ls();
+      TH1F* hAmp = (TH1F*) gFile->Get("hAmpLaser");
+      Int_t nmips=0;
+      for (Int_t ibin=0; ibin<2000; ibin++) {
+	Float_t bincont = hAmp->GetBinContent(ibin);
+	if(bincont>0){ 
+	  mips[nmips] = hAmp->GetXaxis()->GetBinCenter(ibin);
+	  cout<<ibin<<" bincont "<<bincont<<" amp "<< mips[nmips]<<endl;
+	  nmips++;
+	}	
+      }    
+      
+      if (nmips<17) {
+	ok=false;
+	return ok;
+      } 
+       
+      Float_t x1[50], y1[50]; 
+      Float_t x2[50], xx2[50],y2[50];
+      Float_t xx1[50],yy1[50], xx[50];
+      
+      Float_t cfd0 = 0;
+      
+      for (Int_t ii=0; ii<nmips; ii++)
+	x1[ii] = y1[ii] = x2[ii] = y2[ii] = 0; 
+      
+      for (Int_t i=0; i<24; i++)
 	{
-	  if(iin>0)
-	    time[isum] = Float_t (sum/(iin));
-	  else
-	    time[isum] =Float_t (x[ind]);
-	  amplitude[isum] = Float_t (amp);
-	  amp=y[ind];
-	  iin=0;
-	  isum++;
-	  sum=0;
+	  cfd0 = 0;
+	  for (Int_t im=0; im<nmips; im++)
+	    {	      
+	      TString cfd = Form("hCFD%i_%i",i+1,im+1);
+	      TString qtc = Form("hQTC%i_%i",i+1,im+1);
+	      TString led = Form("hLED%i_%i",i+1,im+1);
+	      
+	      TH1F *hCFD = (TH1F*) gFile->Get(cfd.Data()) ;
+	      TH1F *hLED = (TH1F*) gFile->Get(led.Data());
+	      TH1F *hQTC = (TH1F*) gFile->Get(qtc.Data()) ;
+	      //	      hCFD->SetDirectory(0);
+	      //	      hQTC->SetDirectory(0);
+	      //	      hLED->SetDirectory(0);
+	      if(!hCFD )
+ 	      	AliWarning(Form(" no CFD data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
+	      if(!hQTC )
+ 	      	AliWarning(Form(" no QTC correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
+	      if(!hLED)	      
+	      	AliWarning(Form(" no LED correction data in LASER DA for channel %i for amplitude %f MIPs",i,mips[im]));
+	      if( hCFD && hCFD->GetEntries()<500 ) {
+		ok=false;
+		printf("no peak in CFD spectrum for PMT %i amplitude %i\n",i,im);
+		return ok;
+	      }
+	      if(hCFD && hCFD->GetEntries()>500 ) {
+		if( hCFD->GetRMS() >= 1.5) 
+		  GetMeanAndSigma(hCFD, cfdmean, sigma);
+		else
+		  cfdmean = hCFD->GetMean();
+		
+		Int_t   maxBin = hCFD->GetMaximumBin(); 
+		Double_t  meanEstimate = hCFD->GetBinCenter( maxBin); 
+		if(TMath::Abs(meanEstimate - cfdmean) > 20 ) cfdmean = meanEstimate; 
+		if (im == 0) cfd0 = cfdmean;
+		y1[im] =  cfdmean - cfd0;
+	      }	
+	      if(hQTC && hQTC->GetEntries()>500) {
+		GetMeanAndSigma(hQTC, qtmean, sigma);
+		
+		x1[im] = qtmean;
+		if( x1[im] == 0) {
+		  ok=false;
+		  printf("no peak in QTC signal for PMT %i amplitude %i\n",i,im);
+		  return ok;
+		}
+	      }
+		
+	      if( hLED && hLED->GetEntries()>500) {
+		GetMeanAndSigma(hLED, ledmean, sigma);
+	      }				   
+	      else
+		{ 
+		  ok=false;
+		  printf("no peak in LED spectrum for PMT %i amplitude %i\n",i,im);
+		  return ok;
+		}
+	      x2[im] = ledmean;
+	      xx2[im] = x2[nmips-im-1];
+		
+		xx[im]=mips[im];
+		if (hQTC) delete  hQTC;
+		if (hCFD) delete  hCFD;
+		if (hLED) delete  hLED;
+		
+	    }
+	  
+	  for (Int_t imi=0; imi<nmips; imi++)
+	    {
+	      yy1[imi] = Float_t (mips[nmips-imi-1]);
+	      xx1[imi] = x2[nmips-imi-1]; 
+	      //	      cout<<" LED rev "<<i<<" "<<imi<<" "<<xx1[imi]<<" "<< yy1[imi]<<"nmips-imi-1 = "<< nmips-imi-1<<endl;
+	    }
+	  
+	  if(i==0) cout<<"Making graphs..."<<endl;
+	  
+      TGraph *grwalkqtc = new TGraph (nmips,x1,y1);
+      grwalkqtc->SetTitle(Form("PMT%i",i));
+      TGraph *grwalkled = new TGraph (nmips,x2,y1);
+      grwalkled->SetTitle(Form("PMT%i",i));
+      fWalk.AddAtAndExpand(grwalkqtc,i);
+      fAmpLEDRec.AddAtAndExpand(grwalkled,i);
+      //	  cout<<" add walk "<<i<<endl;
+      
+      //fit amplitude graphs to make comparison wth new one	  
+      TGraph *grampled = new TGraph (nmips,xx1,yy1);
+      TGraph *grqtc = new TGraph (nmips,x1,xx);
+      fQTC.AddAtAndExpand(grqtc,i);	 
+      fAmpLED.AddAtAndExpand(grampled,i);
+      //	  cout<<" add amp "<<i<<endl;
+      
+      if(i==23)
+	cout<<"Graphs created..."<<endl;   
 	}
-    }
-
-  inFile.close();
+    } //if gFile exits
 
-  TGraph* gr = new TGraph(isum, amplitude, time);
-  fWalk.AddAtAndExpand(gr,ipmt);
+  return ok;
 }
 
-//________________________________________________________________
-
-void AliT0CalibWalk::SetAmpLEDRec(Int_t ipmt)
-{
-  // read LED walk from external file
-  Float_t mv, ps; 
-  Float_t x[100], y[100];
-  string buffer;
-  
-  const char * filename = gSystem->ExpandPathName("$ALICE_ROOT/T0/data/CFD-LED.txt");
-   ifstream inFile(filename);
-  if(!inFile) {AliError(Form("Cannot open file %s !",filename));}
-  
-  inFile >> mv>>ps;
-  Int_t i=0;
-  
-  while(getline(inFile,buffer)){
-    x[i]=mv; y[i]=ps;	
-    inFile >> mv >> ps;
-    i++;
-  }
-  inFile.close();
-  Float_t y1[100], x1[100];
-  for (Int_t ir=0; ir<i; ir++){
-    y1[ir]=y[i-ir]; x1[ir]=x[i-ir];}
-  TGraph* gr = new TGraph(i,y1,x1);
-  fAmpLEDRec.AddAtAndExpand(gr,ipmt);
-  
-}
-
-//________________________________________________________________
-
-void AliT0CalibWalk::MakeWalkCorrGraph(const char *laserFile)
+void AliT0CalibWalk::GetMeanAndSigma(TH1F* hist, Float_t &mean, Float_t &sigma) 
 {
-  //make walk corerction for preprocessor
 
-  TFile *gFile = TFile::Open(laserFile);
-  //gSystem->Load("libSpectrum");
+  const double window = 2.;  //fit window 
  
-  Int_t npeaks = 20;
-  Int_t sigma=3.;
-  Bool_t down = false;
+  double meanEstimate, sigmaEstimate; 
+  int maxBin;
+  maxBin        =  hist->GetMaximumBin(); //position of maximum
+  meanEstimate  =  hist->GetBinCenter( maxBin); // mean of gaussian sitting in maximum
+  sigmaEstimate = hist->GetRMS();
+  TF1* fit= new TF1("fit","gaus", meanEstimate - window*sigmaEstimate, meanEstimate + window*sigmaEstimate);
+  fit->SetParameters(hist->GetBinContent(maxBin), meanEstimate, sigmaEstimate);
+  hist->Fit("fit","RQ","Q");
 
-  Int_t index[20];
-  Char_t buf1[10], buf2[10],buf3[20];
-  
-  
-  for (Int_t i=0; i<24; i++)
-    {
-      if(i>11)	{  sprintf(buf1,"T0_A_%i_CFD",i+1-12); }
-
-      if(i<12)	{ sprintf(buf1,"T0_C_%i_CFD",i+1); }
-
-      sprintf(buf2,"CFD_QTC%i",i+1);
-      sprintf(buf3,"CFD_LED%i",i+1);
-      cout<<buf1<<" "<<buf2<<" "<<buf3<<endl;
-      TH2F *qtc_cfd = (TH2F*) gFile->Get(buf2);
-      TH2F *led_cfd = (TH2F*) gFile->Get(buf3);
-      TH1F *cfd = (TH1F*) gFile->Get(buf1);
-      //       cfd->Draw();
-      TSpectrum *s = new TSpectrum(2*npeaks,1);
-      Int_t nfound = s->Search(cfd,sigma," ",0.05);
-      cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
-      if(nfound!=0){
-	Float_t *xpeak = s->GetPositionX();
-  	TMath::Sort(nfound, xpeak, index,down);
-  	Float_t xp = xpeak[index[0]];
-	Int_t xbin = cfd->GetXaxis()->FindBin(xp);
-	Float_t yp = cfd->GetBinContent(xbin);
-	cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
-	Float_t hmax = xp+10*sigma;
-	Float_t hmin = xp-10*sigma;
-	cout<<hmin<< " "<<hmax<<endl;
-	//	    cfd->GetXaxis()->SetRange(hmin,hmax);
-	//	    TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
-	// cfd->Fit("g1","R");
-	//	Int_t hminbin=  qtc_cfd->GetXaxis()->GetFirst();
-	//	Int_t hmaxbin=  qtc_cfd->GetXaxis()->GetLast();
-	Int_t nbins= qtc_cfd->GetXaxis()->GetNbins();
-	cout<<" nbins "<<nbins<<endl;
-	//	cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
-	//  qtc_cfd->Draw();
-	TProfile *pr_y = qtc_cfd->ProfileX();
-	//	Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
-	
-	pr_y->SetMaximum(hmax);
-	pr_y->SetMinimum(hmin);
-	Int_t np=nbins/20;
-	Double_t *xx = new Double_t[np];
-	Double_t *yy = new Double_t[np];
-	Int_t ng=0;
-	Double_t yg=0;
-	for (Int_t ip=1; ip<nbins; ip++)
-	  {
-	    if(ip%20 != 0 ) {
-	      if (pr_y->GetBinContent(ip) !=0)
-		yg +=pr_y->GetBinContent(ip);
-	      //	cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
-	      ng++;}
-	    else {
-	      xx[ip/20] = Float_t (pr_y->GetBinCenter(ip));
-	      yy[ip/20] = yg/ng;
-	      yg=0;
-	      ng=0;
-	      //	     	      cout<<ip<<" "<<ip/20<<" "<< xx[ip/20]<<" "<< yy[ip/20]<<endl;
-	    }
-	  }
-	TGraph *gr = new TGraph(np,xx,yy);
-	gr->SetMinimum(hmin);
-	gr->SetMaximum(hmax);
-	fWalk.AddAtAndExpand(gr,i);	  
-	//	fWalk.AddAtAndExpand(gr,i+12);
+  mean  = (Float_t) fit->GetParameter(1);
+  sigma = (Float_t) fit->GetParameter(2);
+  printf(" mean %f max %f sigma %f \n",mean, meanEstimate , sigma);  
 
-	Int_t nbinsled= led_cfd->GetXaxis()->GetNbins();
-	cout<<" nbins led "<<nbinsled<<endl;
-	//	cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
-	//  qtc_cfd->Draw();
-	TProfile *prled_y = led_cfd->ProfileX();
-	//	Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
-	
-	prled_y->SetMaximum(hmax);
-	prled_y->SetMinimum(hmin);
-	Int_t npled=nbinsled/20;
-	Double_t *xxled = new Double_t[np];
-	Double_t *yyled = new Double_t[np];
-	Int_t ngled=0;
-	Double_t ygled=0;
-	for (Int_t ip=1; ip<nbinsled; ip++)
-	  {
-	    if(ip%20 != 0 ) {
-	      if (prled_y->GetBinContent(ip) !=0)
-		ygled +=prled_y->GetBinContent(ip);
-	      //	cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
-	      ngled++;}
-	    else {
-	      xxled[ip/20] = Float_t (prled_y->GetBinCenter(ip));
-	      yyled[ip/20] = ygled/ngled;
-	      ygled=0;
-	      ngled=0;
-	    }
-	  }
-	TGraph *grled = new TGraph(npled,xxled,yyled);
-	grled->SetMinimum(hmin);
-	grled->SetMaximum(hmax);
-	fAmpLEDRec.AddAtAndExpand(grled,i);	  
- 	//	delete [] xx;
-	//	delete [] yy;
-	
-      }
-      
-    }
-  
+  delete fit;
+}
 
 
-}