[u/mrichter/AliRoot.git] / TRD / TRDbase / AliTRDCalibraVdriftLinearFit.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

////////////////////////////////////////////////////////////////////////////
//                                                                        //
// AliTRDCalibraVdriftLinearFit                                           //
//                                                                        //
// Does the Vdrift an ExB calibration by applying a linear fit            //
//                                                                        //
// Author:                                                                //
//   R. Bailhache (R.Bailhache@gsi.de)                                    //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

//Root includes
#include <TObjArray.h>
#include <TH2F.h>
#include <TString.h>
#include <TVectorD.h>
#include <TAxis.h>
#include <TLinearFitter.h>
#include <TMath.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TTreeStream.h>
#include <TGraphErrors.h>
#include <TDirectory.h>
#include <TTreeStream.h>
#include <TF1.h>

//header file
#include "AliTRDCalibraVdriftLinearFit.h"

ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/

//_____________________________________________________________________
AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
  TObject(),
  fVersion(0),
  fNameCalibUsed(""),
  fLinearFitterHistoArray(540),
  fLinearFitterPArray(540),
  fLinearFitterEArray(540),
  fRobustFit(kTRUE),
  fMinNpointsFit(11),
  fNbBindx(32),
  fNbBindy(70),
  fRangedx(0.8),
  fRangedy(1.4),
  fDebugStreamer(0x0),
  fDebugLevel(0),
  fSeeDetector(0)
{
  //
  // default constructor
  //
}
//_____________________________________________________________________
AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
  TObject(ped),
  fVersion(ped.fVersion),
  fNameCalibUsed(ped.fNameCalibUsed),
  fLinearFitterHistoArray(540),
  fLinearFitterPArray(540),
  fLinearFitterEArray(540),
  fRobustFit(kTRUE),
  fMinNpointsFit(10),
  fNbBindx(32),
  fNbBindy(70),
  fRangedx(0.8),
  fRangedy(1.4),
  fDebugStreamer(0x0),
  fDebugLevel(0),
  fSeeDetector(0)
{
    //
    // copy constructor
    //
  for (Int_t idet = 0; idet < 540; idet++){
   
    const TVectorD     *vectorE     = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
    const TVectorD     *vectorP     = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
    const TH2S         *hped        = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
    
    if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
    if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
    if ( hped != 0x0 ){
      TH2S *hNew = (TH2S *)hped->Clone();
      //hNew->SetDirectory(0);
      fLinearFitterHistoArray.AddAt(hNew,idet);
    }
  }
}
//_____________________________________________________________________
AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
  TObject(),
  fVersion(0),
  fNameCalibUsed(""),
  fLinearFitterHistoArray(540),
  fLinearFitterPArray(540),
  fLinearFitterEArray(540),
  fRobustFit(kTRUE),
  fMinNpointsFit(10),
  fNbBindx(32),
  fNbBindy(70),
  fRangedx(0.8),
  fRangedy(1.4),
  fDebugStreamer(0x0),
  fDebugLevel(0),
  fSeeDetector(0)
{
  //
  // constructor from a TObjArray
  //
  for (Int_t idet = 0; idet < 540; idet++){
    const TH2S         *hped        = (TH2S*)obja.UncheckedAt(idet);
    if ( hped != 0x0 ){
      TH2S *hNew = (TH2S *)hped->Clone();
      //hNew->SetDirectory(0);
      fLinearFitterHistoArray.AddAt(hNew,idet);
    }
  }
}
//_____________________________________________________________________
AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const  AliTRDCalibraVdriftLinearFit &source)
{
  //
  // assignment operator
  //
  if (&source == this) return *this;
  new (this) AliTRDCalibraVdriftLinearFit(source);

  return *this;
}
//_____________________________________________________________________
AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
{
  //
  // destructor
  //
  fLinearFitterHistoArray.SetOwner();
  fLinearFitterPArray.SetOwner();
  fLinearFitterEArray.SetOwner();

  fLinearFitterHistoArray.Delete();
  fLinearFitterPArray.Delete();
  fLinearFitterEArray.Delete();

  if ( fDebugStreamer ) delete fDebugStreamer;

}
//_____________________________________________________________________________
void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
{
  //
  // Copy function
  //

  AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;

  // Copy only the histos
  for (Int_t idet = 0; idet < 540; idet++){
    if(fLinearFitterHistoArray.UncheckedAt(idet)){
      TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
      //hped1->SetDirectory(0);
      hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
    }
  }
  
  TObject::Copy(c);

}
//_____________________________________________________________________________
Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list) 
{
  // Merge list of objects (needed by PROOF)

  if (!list)
    return 0;
  
  if (list->IsEmpty())
    return 1;
  
  TIterator* iter = list->MakeIterator();
  TObject* obj = 0;
  
  // collection of generated histograms
  Int_t count=0;
  while((obj = iter->Next()) != 0) 
    {
      AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
      if (entry == 0) continue; 
      
      // Copy only the histos
      for (Int_t idet = 0; idet < 540; idet++){
	if(entry->GetLinearFitterHisto(idet)){
	  TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
	  Double_t entriesa = hped1->GetEntries();
	  Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
	  if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
	}
      }
      
      count++;
    }
  
  return count;
}
//_____________________________________________________________________
void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
{
  //
  // Add histo
  //

  fVersion++;

  for (Int_t idet = 0; idet < 540; idet++){
    const TH2S         *hped        = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
    //printf("idet %d\n",idet);
    if ( hped != 0x0 ){
      //printf("add\n");
      TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
      Double_t entriesa = hped1->GetEntries();
      Double_t entriesb = hped->GetEntries();
      if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
    }
  }
}
//______________________________________________________________________________________
TH2S* AliTRDCalibraVdriftLinearFit::AddAll() 
{
    //
    // return pointer to TH2S of all added histos 
  //

  TH2S *histo = 0x0;
  for(Int_t k=0; k < 540; k++) {
    TH2S * u = GetLinearFitterHistoForce(k);
    if(k == 0) {
      histo = new TH2S(*u);
      histo->SetName("sum");
    }
    else histo->Add(u);
  }

  return histo;

}
//______________________________________________________________________________________
TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
{
    //
    // return pointer to TH2F histo 
    // if force is true create a new histo if it doesn't exist allready
    //
    if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
	return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);

    return GetLinearFitterHistoForce(detector);

}
//______________________________________________________________________________________
TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
{
  //
  // return pointer to TH2F histo 
  // if NULL create a new histo if it doesn't exist allready
  //
  if (fLinearFitterHistoArray.UncheckedAt(detector))
    return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
  
  // if we are forced and TLinearFitter doesn't yes exist create it
  
  // new TH2F
  TString name("LFDV");
  name += detector;
  name += "version";
  name +=  fVersion;
  
  TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
			,(Int_t)fNbBindx,-fRangedx,fRangedx,(Int_t)fNbBindy
			,-fRangedy,fRangedy);
  lfdv->SetXTitle("tan(phi_{track})");
  lfdv->SetYTitle("dy/dt");
  lfdv->SetZTitle("Number of clusters");
  lfdv->SetStats(0);
  lfdv->SetDirectory(0);
  
  fLinearFitterHistoArray.AddAt(lfdv,detector);
  return lfdv;
}
//______________________________________________________________________________________
Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
{
    //
    // return param for this detector
    //
  if ( fLinearFitterPArray.UncheckedAt(detector) ){
    const TVectorD     *vectorP     = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
    if(!param) param = new TVectorD(2);
    for(Int_t k = 0; k < 2; k++){
      (*param)[k] = (*vectorP)[k];
    }
    return kTRUE;
  }
  else return kFALSE;

}
//______________________________________________________________________________________
Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
{
    //
    // return error for this detector 
    //
  if ( fLinearFitterEArray.UncheckedAt(detector) ){
    const TVectorD     *vectorE     = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
    if(!error) error = new TVectorD(3);
    for(Int_t k = 0; k < 3; k++){
      (*error)[k] = (*vectorE)[k];
    }
    return kTRUE;
  }
  else return kFALSE;

}
//______________________________________________________________________________________
void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
{
    //
    // Fill the 2D histos for debugging
    //
  
  TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
  Double_t nbentries = h->GetEntries();
  if(nbentries < 5*32767) h->Fill(tnp,pars1);

}
//____________Functions fit Online CH2d________________________________________
void AliTRDCalibraVdriftLinearFit::FillPEArray()
{
  //
  // Fill fLinearFitterPArray and fLinearFitterEArray from inside
  //

  
  Int_t *arrayI = new Int_t[540];
  for(Int_t k = 0; k< 540; k++){
    arrayI[k] = 0; 
  }

  // Loop over histos 
  for(Int_t cb = 0; cb < 540; cb++){
    const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
    //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);    

    if ( linearfitterhisto != 0 ){
      
      // Fill a linearfitter
      const TAxis *xaxis = linearfitterhisto->GetXaxis();
      const TAxis *yaxis = linearfitterhisto->GetYaxis();
      TLinearFitter linearfitter = TLinearFitter(2,"pol1");
      //printf("test\n");
      Double_t integral = linearfitterhisto->Integral();
      //printf("Integral is %f\n",integral);
      Bool_t securitybreaking = kFALSE;
      if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
      for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
	for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
	  if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
	    Double_t x = xaxis->GetBinCenter(ibinx+1);
	    Double_t y = yaxis->GetBinCenter(ibiny+1);
	    
	    for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
	      if(!securitybreaking){
		linearfitter.AddPoint(&x,y);
		arrayI[cb]++;
	      }
	      else {
		if(arrayI[cb]< 1198){
		  linearfitter.AddPoint(&x,y);
		  arrayI[cb]++; 
		}
	      }
	    }
	    
	  }
	}
      }
      
      //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);

      // Eval the linear fitter
      if(arrayI[cb]>10){
	TVectorD  *par  = new TVectorD(2);
	TVectorD   pare = TVectorD(2);
	TVectorD  *parE = new TVectorD(3);
	//printf("Fit\n");
	if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (linearfitter.Eval()==0))) {
	  //if((linearfitter.Eval()==0)) {
	  //printf("Take the param\n");
	  linearfitter.GetParameters(*par);
	  //printf("Done\n");
	  //linearfitter.GetErrors(pare);
	  //Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
	  //(*parE)[0] = pare[0]*ppointError;
	  //(*parE)[1] = pare[1]*ppointError;
	  
	  (*parE)[0] = 0.0;
	  (*parE)[1] = 0.0;
	  (*parE)[2] = (Double_t) arrayI[cb];
	  fLinearFitterPArray.AddAt(par,cb);
	  fLinearFitterEArray.AddAt(parE,cb);
	  
	  //par->Print();
	  //parE->Print();
	}
	//printf("Finish\n");
      }
      
      //delete linearfitterhisto;
      
    }// if something

  }

  delete [] arrayI;
   
}

//____________Functions fit Online CH2d________________________________________
void AliTRDCalibraVdriftLinearFit::FillPEArray2()
{
  //
  // Fill fFitterPArray and fFitterEArray from inside
  //

  // Loop over histos 
  TF1 *f1 = new TF1("f1","[0]+[1]*x",-1,1);
      
  for(Int_t cb = 0; cb < 540; cb++){
    //const TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
    TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
    //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) fitterhisto);    
  
    if ( fitterhisto != 0 ){
      
      Int_t nEntries=0;
      TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
      if (!gg) continue;
      // Number of points of the TGraphErrors
      //printf("Number of points %d for detector %d\n",gg->GetN(),cb);
      if(gg->GetN() <  fMinNpointsFit) {
	if(gg) delete gg;
	continue;	
      }
      //printf("det %d, number of points %d, nEntries %d\n",cb,gg->GetN(),nEntries);
      gg->Fit(f1,"Q0");
      
      TVectorD  *par  = new TVectorD(2);
      TVectorD  *parE = new TVectorD(3);
      (*parE)[0] = 0.0;
      (*parE)[1] = 0.0;
      (*parE)[2] = (Double_t) nEntries;
      (*par)[0] = f1->GetParameter(0);
      (*par)[1] = f1->GetParameter(1);
      fLinearFitterPArray.AddAt(par,cb);
      fLinearFitterEArray.AddAt(parE,cb);
      //printf("Value %f for detector %d\n",(*par)[1],cb);
            
      if(fDebugLevel==0) {
	if(gg) delete gg;
      }
      else {
	if(cb==fSeeDetector) {
	  gStyle->SetPalette(1);
	  gStyle->SetOptStat(1111);
	  gStyle->SetPadBorderMode(0);
	  gStyle->SetCanvasColor(10);
	  gStyle->SetPadLeftMargin(0.13);
	  gStyle->SetPadRightMargin(0.01);
	  TCanvas *stat = new TCanvas("stat","",50,50,600,800);
	  stat->cd(1);
	  fitterhisto->Draw("colz");
	  gg->Draw("P");
	  f1->Draw("same");
	  break;
	}
      } 
    }
  }
  if(fDebugLevel==0) delete f1;
}

//_________Helper function__________________________________________________
TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
{
  //
  // Debug function
  //

  TF1 fg("fg", "gaus", -10., 30.);
  TGraphErrors *gp = new TGraphErrors();

  const TAxis *ax(h2->GetXaxis());
  const TAxis *ay(h2->GetYaxis());
  TH1D *h1(NULL);
  for(Int_t ipt(0), jpt(1), ig(0); ipt<ax->GetNbins(); ipt++, jpt++){
    h1 = h2->ProjectionY("py", jpt, jpt);
    fg.SetParameter(1, h1->GetMean());
    //Float_t x(ax->GetBinCenter(jpt)); 
    Int_t n=Int_t(h1->Integral(1, h1->GetNbinsX()));
    nEntries+=n;
    if(n < 15){
      //Warning("drawMS()", Form("reject x[%d]=%f on n=%d", jpt, x, n));
      continue;
    }
    h1->Fit(&fg, "WWQ0");
    if(fg.GetNDF()<2){
      //Warning("drawMS()", Form("reject x[%d]=%f on NDF=%d", jpt, x, fg.GetNDF()));
      continue;
    }
    if(((fg.GetParameter(1)+fg.GetParameter(2)/2)>ay->GetXmax()) || ((fg.GetParameter(1)-fg.GetParameter(2)/2)<ay->GetXmin()) || (TMath::Abs(fg.GetParameter(0))< 0.00001)) continue;
    gp->SetPoint(ig, ax->GetBinCenter(jpt), fg.GetParameter(1));
    gp->SetPointError(ig, 0, TMath::Sqrt(pow(fg.GetParError(1),2) + (1/pow(fg.GetParameter(0),2))));
    ig++;
  }
  delete h1;
  return gp;
}
Commit	Line	Data
3a0f6479	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
0bc7827a	18	////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// AliTRDCalibraVdriftLinearFit //
	21	// //
	22	// Does the Vdrift an ExB calibration by applying a linear fit //
	23	// //
	24	// Author: //
	25	// R. Bailhache (R.Bailhache@gsi.de) //
	26	// //
	27	////////////////////////////////////////////////////////////////////////////
	28
3a0f6479	29	//Root includes
	30	#include <TObjArray.h>
	31	#include <TH2F.h>
	32	#include <TString.h>
	33	#include <TVectorD.h>
	34	#include <TAxis.h>
	35	#include <TLinearFitter.h>
	36	#include <TMath.h>
661c7be6	37	#include <TStyle.h>
661c7be6	38	#include <TCanvas.h>
3a0f6479	39	#include <TTreeStream.h>
661c7be6	40	#include <TGraphErrors.h>
	41	#include <TDirectory.h>
	42	#include <TTreeStream.h>
	43	#include <TF1.h>
3a0f6479	44
	45	//header file
	46	#include "AliTRDCalibraVdriftLinearFit.h"
	47
	48	ClassImp(AliTRDCalibraVdriftLinearFit) /FOLD00/
	49
	50	//_____________________________________________________________________
	51	AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /FOLD00/
	52	TObject(),
	53	fVersion(0),
b2277aa2	54	fNameCalibUsed(""),
3a0f6479	55	fLinearFitterHistoArray(540),
3a0f6479	56	fLinearFitterPArray(540),
2a1a7b36	57	fLinearFitterEArray(540),
661c7be6	58	fRobustFit(kTRUE),
532be2d4	59	fMinNpointsFit(11),
4f3bd513	60	fNbBindx(32),
	61	fNbBindy(70),
	62	fRangedx(0.8),
	63	fRangedy(1.4),
661c7be6	64	fDebugStreamer(0x0),
	65	fDebugLevel(0),
	66	fSeeDetector(0)
3a0f6479	67	{
ba94744a	68	//
	69	// default constructor
	70	//
3a0f6479	71	}
3a0f6479	72	//_____________________________________________________________________
	73	AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /FOLD00/
	74	TObject(ped),
	75	fVersion(ped.fVersion),
b2277aa2	76	fNameCalibUsed(ped.fNameCalibUsed),
3a0f6479	77	fLinearFitterHistoArray(540),
3a0f6479	78	fLinearFitterPArray(540),
2a1a7b36	79	fLinearFitterEArray(540),
661c7be6	80	fRobustFit(kTRUE),
532be2d4	81	fMinNpointsFit(10),
4f3bd513	82	fNbBindx(32),
	83	fNbBindy(70),
	84	fRangedx(0.8),
	85	fRangedy(1.4),
661c7be6	86	fDebugStreamer(0x0),
	87	fDebugLevel(0),
	88	fSeeDetector(0)
3a0f6479	89	{
	90	//
	91	// copy constructor
	92	//
	93	for (Int_t idet = 0; idet < 540; idet++){
ba94744a	94
3a0f6479	95	const TVectorD vectorE = (TVectorD)ped.fLinearFitterEArray.UncheckedAt(idet);
3a0f6479	96	const TVectorD vectorP = (TVectorD)ped.fLinearFitterPArray.UncheckedAt(idet);
0ca1720e	97	const TH2S hped = (TH2S)ped.fLinearFitterHistoArray.UncheckedAt(idet);
3a0f6479	98
	99	if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
	100	if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
	101	if ( hped != 0x0 ){
0ca1720e	102	TH2S hNew = (TH2S )hped->Clone();
ba94744a	103	//hNew->SetDirectory(0);
3a0f6479	104	fLinearFitterHistoArray.AddAt(hNew,idet);
3a0f6479	105	}
3a0f6479	106	}
	107	}
	108	//_____________________________________________________________________
d0569428	109	AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /FOLD00/
	110	TObject(),
	111	fVersion(0),
b2277aa2	112	fNameCalibUsed(""),
d0569428	113	fLinearFitterHistoArray(540),
d0569428	114	fLinearFitterPArray(540),
2a1a7b36	115	fLinearFitterEArray(540),
661c7be6	116	fRobustFit(kTRUE),
532be2d4	117	fMinNpointsFit(10),
4f3bd513	118	fNbBindx(32),
	119	fNbBindy(70),
	120	fRangedx(0.8),
	121	fRangedy(1.4),
661c7be6	122	fDebugStreamer(0x0),
	123	fDebugLevel(0),
	124	fSeeDetector(0)
d0569428	125	{
	126	//
	127	// constructor from a TObjArray
	128	//
	129	for (Int_t idet = 0; idet < 540; idet++){
0ca1720e	130	const TH2S hped = (TH2S)obja.UncheckedAt(idet);
d0569428	131	if ( hped != 0x0 ){
0ca1720e	132	TH2S hNew = (TH2S )hped->Clone();
ba94744a	133	//hNew->SetDirectory(0);
d0569428	134	fLinearFitterHistoArray.AddAt(hNew,idet);
	135	}
	136	}
	137	}
	138	//_____________________________________________________________________
3a0f6479	139	AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const AliTRDCalibraVdriftLinearFit &source)
	140	{
	141	//
	142	// assignment operator
	143	//
	144	if (&source == this) return *this;
	145	new (this) AliTRDCalibraVdriftLinearFit(source);
	146
	147	return *this;
	148	}
	149	//_____________________________________________________________________
	150	AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /FOLD00/
	151	{
	152	//
	153	// destructor
	154	//
f5cda287	155	fLinearFitterHistoArray.SetOwner();
	156	fLinearFitterPArray.SetOwner();
	157	fLinearFitterEArray.SetOwner();
	158
1ca79a00	159	fLinearFitterHistoArray.Delete();
	160	fLinearFitterPArray.Delete();
	161	fLinearFitterEArray.Delete();
	162
661c7be6	163	if ( fDebugStreamer ) delete fDebugStreamer;
661c7be6	164
3a0f6479	165	}
64942b85	166	//_____________________________________________________________________________
	167	void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
	168	{
	169	//
	170	// Copy function
	171	//
	172
	173	AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
	174
	175	// Copy only the histos
	176	for (Int_t idet = 0; idet < 540; idet++){
	177	if(fLinearFitterHistoArray.UncheckedAt(idet)){
0ca1720e	178	TH2S hped1 = (TH2S )target.GetLinearFitterHisto(idet,kTRUE);
ba94744a	179	//hped1->SetDirectory(0);
0ca1720e	180	hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
64942b85	181	}
64942b85	182	}
ba94744a	183
64942b85	184	TObject::Copy(c);
	185
	186	}
	187	//_____________________________________________________________________________
6bbdc11a	188	Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list)
64942b85	189	{
	190	// Merge list of objects (needed by PROOF)
	191
	192	if (!list)
	193	return 0;
	194
	195	if (list->IsEmpty())
	196	return 1;
	197
	198	TIterator* iter = list->MakeIterator();
	199	TObject* obj = 0;
	200
	201	// collection of generated histograms
	202	Int_t count=0;
	203	while((obj = iter->Next()) != 0)
ba94744a	204	{
	205	AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
	206	if (entry == 0) continue;
	207
	208	// Copy only the histos
	209	for (Int_t idet = 0; idet < 540; idet++){
	210	if(entry->GetLinearFitterHisto(idet)){
0ca1720e	211	TH2S hped1 = (TH2S )GetLinearFitterHisto(idet,kTRUE);
ba94744a	212	Double_t entriesa = hped1->GetEntries();
0ca1720e	213	Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
ba94744a	214	if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
ba94744a	215	}
64942b85	216	}
64942b85	217
ba94744a	218	count++;
64942b85	219	}
64942b85	220
	221	return count;
	222	}
3a0f6479	223	//_____________________________________________________________________
fe4ee353	224	void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
3a0f6479	225	{
	226	//
	227	// Add histo
	228	//
	229
	230	fVersion++;
	231
	232	for (Int_t idet = 0; idet < 540; idet++){
fe4ee353	233	const TH2S hped = (TH2S)ped->GetLinearFitterHistoNoForce(idet);
3a0f6479	234	//printf("idet %d\n",idet);
	235	if ( hped != 0x0 ){
	236	//printf("add\n");
0ca1720e	237	TH2S hped1 = (TH2S )GetLinearFitterHisto(idet,kTRUE);
ba94744a	238	Double_t entriesa = hped1->GetEntries();
	239	Double_t entriesb = hped->GetEntries();
	240	if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
3a0f6479	241	}
	242	}
	243	}
	244	//______________________________________________________________________________________
4f3bd513	245	TH2S* AliTRDCalibraVdriftLinearFit::AddAll()
	246	{
	247	//
	248	// return pointer to TH2S of all added histos
	249	//
	250
	251	TH2S *histo = 0x0;
	252	for(Int_t k=0; k < 540; k++) {
	253	TH2S * u = GetLinearFitterHistoForce(k);
	254	if(k == 0) {
	255	histo = new TH2S(*u);
	256	histo->SetName("sum");
	257	}
	258	else histo->Add(u);
	259	}
	260
	261	return histo;
	262
	263	}
	264	//______________________________________________________________________________________
0ca1720e	265	TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
3a0f6479	266	{
	267	//
	268	// return pointer to TH2F histo
	269	// if force is true create a new histo if it doesn't exist allready
	270	//
	271	if ( !force \|\| fLinearFitterHistoArray.UncheckedAt(detector) )
0ca1720e	272	return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
3a0f6479	273
fe4ee353	274	return GetLinearFitterHistoForce(detector);
	275
	276	}
	277	//______________________________________________________________________________________
	278	TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
	279	{
	280	//
	281	// return pointer to TH2F histo
	282	// if NULL create a new histo if it doesn't exist allready
	283	//
	284	if (fLinearFitterHistoArray.UncheckedAt(detector))
	285	return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
	286
	287	// if we are forced and TLinearFitter doesn't yes exist create it
	288
	289	// new TH2F
	290	TString name("LFDV");
	291	name += detector;
	292	name += "version";
	293	name += fVersion;
	294
	295	TH2S lfdv = new TH2S((const Char_t )name,(const Char_t *) name
4f3bd513	296	,(Int_t)fNbBindx,-fRangedx,fRangedx,(Int_t)fNbBindy
4f3bd513	297	,-fRangedy,fRangedy);
fe4ee353	298	lfdv->SetXTitle("tan(phi_{track})");
	299	lfdv->SetYTitle("dy/dt");
	300	lfdv->SetZTitle("Number of clusters");
	301	lfdv->SetStats(0);
	302	lfdv->SetDirectory(0);
	303
	304	fLinearFitterHistoArray.AddAt(lfdv,detector);
	305	return lfdv;
3a0f6479	306	}
	307	//______________________________________________________________________________________
	308	Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
	309	{
	310	//
	311	// return param for this detector
	312	//
	313	if ( fLinearFitterPArray.UncheckedAt(detector) ){
	314	const TVectorD vectorP = (TVectorD)fLinearFitterPArray.UncheckedAt(detector);
	315	if(!param) param = new TVectorD(2);
	316	for(Int_t k = 0; k < 2; k++){
	317	(param)[k] = (vectorP)[k];
	318	}
	319	return kTRUE;
	320	}
	321	else return kFALSE;
	322
	323	}
	324	//______________________________________________________________________________________
	325	Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
	326	{
	327	//
	328	// return error for this detector
	329	//
	330	if ( fLinearFitterEArray.UncheckedAt(detector) ){
	331	const TVectorD vectorE = (TVectorD)fLinearFitterEArray.UncheckedAt(detector);
	332	if(!error) error = new TVectorD(3);
	333	for(Int_t k = 0; k < 3; k++){
	334	(error)[k] = (vectorE)[k];
	335	}
	336	return kTRUE;
	337	}
	338	else return kFALSE;
	339
	340	}
	341	//______________________________________________________________________________________
	342	void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
	343	{
	344	//
	345	// Fill the 2D histos for debugging
	346	//
ba94744a	347
0ca1720e	348	TH2S h = ((TH2S ) GetLinearFitterHisto(detector,kTRUE));
ba94744a	349	Double_t nbentries = h->GetEntries();
0ca1720e	350	if(nbentries < 5*32767) h->Fill(tnp,pars1);
3a0f6479	351
	352	}
	353	//____________Functions fit Online CH2d________________________________________
	354	void AliTRDCalibraVdriftLinearFit::FillPEArray()
	355	{
	356	//
	357	// Fill fLinearFitterPArray and fLinearFitterEArray from inside
	358	//
	359
	360
	361	Int_t *arrayI = new Int_t[540];
	362	for(Int_t k = 0; k< 540; k++){
	363	arrayI[k] = 0;
	364	}
	365
	366	// Loop over histos
	367	for(Int_t cb = 0; cb < 540; cb++){
0ca1720e	368	const TH2S linearfitterhisto = (TH2S)fLinearFitterHistoArray.UncheckedAt(cb);
3a0f6479	369	//printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);
3a0f6479	370
0ca1720e	371	if ( linearfitterhisto != 0 ){
3a0f6479	372
3a0f6479	373	// Fill a linearfitter
700a5a26	374	const TAxis *xaxis = linearfitterhisto->GetXaxis();
700a5a26	375	const TAxis *yaxis = linearfitterhisto->GetYaxis();
3a0f6479	376	TLinearFitter linearfitter = TLinearFitter(2,"pol1");
3a0f6479	377	//printf("test\n");
a5dcf618	378	Double_t integral = linearfitterhisto->Integral();
	379	//printf("Integral is %f\n",integral);
	380	Bool_t securitybreaking = kFALSE;
	381	if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
3a0f6479	382	for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
3a0f6479	383	for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
8a606eb5	384	if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
3a0f6479	385	Double_t x = xaxis->GetBinCenter(ibinx+1);
3a0f6479	386	Double_t y = yaxis->GetBinCenter(ibiny+1);
a5dcf618	387
3a0f6479	388	for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
a5dcf618	389	if(!securitybreaking){
	390	linearfitter.AddPoint(&x,y);
	391	arrayI[cb]++;
	392	}
	393	else {
	394	if(arrayI[cb]< 1198){
	395	linearfitter.AddPoint(&x,y);
	396	arrayI[cb]++;
	397	}
	398	}
3a0f6479	399	}
a5dcf618	400
3a0f6479	401	}
	402	}
	403	}
	404
	405	//printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
	406
	407	// Eval the linear fitter
	408	if(arrayI[cb]>10){
	409	TVectorD *par = new TVectorD(2);
	410	TVectorD pare = TVectorD(2);
	411	TVectorD *parE = new TVectorD(3);
a5dcf618	412	//printf("Fit\n");
2a1a7b36	413	if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) \|\| ((!fRobustFit) && (linearfitter.Eval()==0))) {
a5dcf618	414	//if((linearfitter.Eval()==0)) {
a5dcf618	415	//printf("Take the param\n");
c1105918	416	linearfitter.GetParameters(*par);
a5dcf618	417	//printf("Done\n");
c1105918	418	//linearfitter.GetErrors(pare);
	419	//Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
	420	//(parE)[0] = pare[0]ppointError;
	421	//(parE)[1] = pare[1]ppointError;
	422
	423	(*parE)[0] = 0.0;
	424	(*parE)[1] = 0.0;
	425	(*parE)[2] = (Double_t) arrayI[cb];
	426	fLinearFitterPArray.AddAt(par,cb);
	427	fLinearFitterEArray.AddAt(parE,cb);
	428
	429	//par->Print();
	430	//parE->Print();
	431	}
a5dcf618	432	//printf("Finish\n");
ba94744a	433	}
c1105918	434
07f9b7f2	435	//delete linearfitterhisto;
ba94744a	436
3a0f6479	437	}// if something
fdc15553	438
3a0f6479	439	}
a987273c	440
a987273c	441	delete [] arrayI;
3a0f6479	442
d1e1a064	443	}
c1105918	444
661c7be6	445	//____________Functions fit Online CH2d________________________________________
	446	void AliTRDCalibraVdriftLinearFit::FillPEArray2()
	447	{
	448	//
	449	// Fill fFitterPArray and fFitterEArray from inside
	450	//
	451
	452	// Loop over histos
	453	TF1 f1 = new TF1("f1","[0]+[1]x",-1,1);
	454
	455	for(Int_t cb = 0; cb < 540; cb++){
	456	//const TH2S fitterhisto = (TH2S)fLinearFitterHistoArray.UncheckedAt(cb);
	457	TH2S fitterhisto = (TH2S)fLinearFitterHistoArray.UncheckedAt(cb);
	458	//printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) fitterhisto);
	459
	460	if ( fitterhisto != 0 ){
	461
	462	Int_t nEntries=0;
	463	TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
ce42378b	464	if (!gg) continue;
661c7be6	465	// Number of points of the TGraphErrors
532be2d4	466	//printf("Number of points %d for detector %d\n",gg->GetN(),cb);
532be2d4	467	if(gg->GetN() < fMinNpointsFit) {
661c7be6	468	if(gg) delete gg;
	469	continue;
	470	}
	471	//printf("det %d, number of points %d, nEntries %d\n",cb,gg->GetN(),nEntries);
	472	gg->Fit(f1,"Q0");
	473
	474	TVectorD *par = new TVectorD(2);
	475	TVectorD *parE = new TVectorD(3);
	476	(*parE)[0] = 0.0;
	477	(*parE)[1] = 0.0;
	478	(*parE)[2] = (Double_t) nEntries;
	479	(*par)[0] = f1->GetParameter(0);
	480	(*par)[1] = f1->GetParameter(1);
	481	fLinearFitterPArray.AddAt(par,cb);
	482	fLinearFitterEArray.AddAt(parE,cb);
532be2d4	483	//printf("Value %f for detector %d\n",(*par)[1],cb);
661c7be6	484
	485	if(fDebugLevel==0) {
	486	if(gg) delete gg;
	487	}
	488	else {
	489	if(cb==fSeeDetector) {
	490	gStyle->SetPalette(1);
	491	gStyle->SetOptStat(1111);
	492	gStyle->SetPadBorderMode(0);
	493	gStyle->SetCanvasColor(10);
	494	gStyle->SetPadLeftMargin(0.13);
	495	gStyle->SetPadRightMargin(0.01);
	496	TCanvas *stat = new TCanvas("stat","",50,50,600,800);
	497	stat->cd(1);
	498	fitterhisto->Draw("colz");
	499	gg->Draw("P");
	500	f1->Draw("same");
	501	break;
	502	}
	503	}
	504	}
	505	}
	506	if(fDebugLevel==0) delete f1;
	507	}
	508
	509	//_________Helper function__________________________________________________
	510	TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
	511	{
92df0132	512	//
	513	// Debug function
	514	//
	515
661c7be6	516	TF1 fg("fg", "gaus", -10., 30.);
661c7be6	517	TGraphErrors *gp = new TGraphErrors();
c1105918	518
700a5a26	519	const TAxis *ax(h2->GetXaxis());
700a5a26	520	const TAxis *ay(h2->GetYaxis());
661c7be6	521	TH1D *h1(NULL);
	522	for(Int_t ipt(0), jpt(1), ig(0); ipt<ax->GetNbins(); ipt++, jpt++){
	523	h1 = h2->ProjectionY("py", jpt, jpt);
	524	fg.SetParameter(1, h1->GetMean());
	525	//Float_t x(ax->GetBinCenter(jpt));
	526	Int_t n=Int_t(h1->Integral(1, h1->GetNbinsX()));
	527	nEntries+=n;
	528	if(n < 15){
	529	//Warning("drawMS()", Form("reject x[%d]=%f on n=%d", jpt, x, n));
	530	continue;
	531	}
	532	h1->Fit(&fg, "WWQ0");
	533	if(fg.GetNDF()<2){
	534	//Warning("drawMS()", Form("reject x[%d]=%f on NDF=%d", jpt, x, fg.GetNDF()));
	535	continue;
	536	}
64c7fdf8	537	if(((fg.GetParameter(1)+fg.GetParameter(2)/2)>ay->GetXmax()) \|\| ((fg.GetParameter(1)-fg.GetParameter(2)/2)<ay->GetXmin()) \|\| (TMath::Abs(fg.GetParameter(0))< 0.00001)) continue;
661c7be6	538	gp->SetPoint(ig, ax->GetBinCenter(jpt), fg.GetParameter(1));
	539	gp->SetPointError(ig, 0, TMath::Sqrt(pow(fg.GetParError(1),2) + (1/pow(fg.GetParameter(0),2))));
	540	ig++;
	541	}
	542	delete h1;
	543	return gp;
	544	}