[u/mrichter/AliRoot.git] / TRD / 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),
  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),
  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),
  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::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
			,36,-0.9,0.9,48
			,-1.2,1.2);
  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
      TAxis *xaxis = linearfitterhisto->GetXaxis();
      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
      if(gg->GetN() < 20) {
	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);
            
      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)
{
  TF1 fg("fg", "gaus", -10., 30.);
  TGraphErrors *gp = new TGraphErrors();

  TAxis *ax(h2->GetXaxis());
  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),
	59	fDebugStreamer(0x0),
	60	fDebugLevel(0),
	61	fSeeDetector(0)
3a0f6479	62	{
ba94744a	63	//
	64	// default constructor
	65	//
3a0f6479	66	}
3a0f6479	67	//_____________________________________________________________________
	68	AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /FOLD00/
	69	TObject(ped),
	70	fVersion(ped.fVersion),
b2277aa2	71	fNameCalibUsed(ped.fNameCalibUsed),
3a0f6479	72	fLinearFitterHistoArray(540),
3a0f6479	73	fLinearFitterPArray(540),
2a1a7b36	74	fLinearFitterEArray(540),
661c7be6	75	fRobustFit(kTRUE),
	76	fDebugStreamer(0x0),
	77	fDebugLevel(0),
	78	fSeeDetector(0)
3a0f6479	79	{
	80	//
	81	// copy constructor
	82	//
	83	for (Int_t idet = 0; idet < 540; idet++){
ba94744a	84
3a0f6479	85	const TVectorD vectorE = (TVectorD)ped.fLinearFitterEArray.UncheckedAt(idet);
3a0f6479	86	const TVectorD vectorP = (TVectorD)ped.fLinearFitterPArray.UncheckedAt(idet);
0ca1720e	87	const TH2S hped = (TH2S)ped.fLinearFitterHistoArray.UncheckedAt(idet);
3a0f6479	88
	89	if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
	90	if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
	91	if ( hped != 0x0 ){
0ca1720e	92	TH2S hNew = (TH2S )hped->Clone();
ba94744a	93	//hNew->SetDirectory(0);
3a0f6479	94	fLinearFitterHistoArray.AddAt(hNew,idet);
3a0f6479	95	}
3a0f6479	96	}
	97	}
	98	//_____________________________________________________________________
d0569428	99	AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /FOLD00/
	100	TObject(),
	101	fVersion(0),
b2277aa2	102	fNameCalibUsed(""),
d0569428	103	fLinearFitterHistoArray(540),
d0569428	104	fLinearFitterPArray(540),
2a1a7b36	105	fLinearFitterEArray(540),
661c7be6	106	fRobustFit(kTRUE),
	107	fDebugStreamer(0x0),
	108	fDebugLevel(0),
	109	fSeeDetector(0)
d0569428	110	{
	111	//
	112	// constructor from a TObjArray
	113	//
	114	for (Int_t idet = 0; idet < 540; idet++){
0ca1720e	115	const TH2S hped = (TH2S)obja.UncheckedAt(idet);
d0569428	116	if ( hped != 0x0 ){
0ca1720e	117	TH2S hNew = (TH2S )hped->Clone();
ba94744a	118	//hNew->SetDirectory(0);
d0569428	119	fLinearFitterHistoArray.AddAt(hNew,idet);
	120	}
	121	}
	122	}
	123	//_____________________________________________________________________
3a0f6479	124	AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const AliTRDCalibraVdriftLinearFit &source)
	125	{
	126	//
	127	// assignment operator
	128	//
	129	if (&source == this) return *this;
	130	new (this) AliTRDCalibraVdriftLinearFit(source);
	131
	132	return *this;
	133	}
	134	//_____________________________________________________________________
	135	AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /FOLD00/
	136	{
	137	//
	138	// destructor
	139	//
f5cda287	140	fLinearFitterHistoArray.SetOwner();
	141	fLinearFitterPArray.SetOwner();
	142	fLinearFitterEArray.SetOwner();
	143
1ca79a00	144	fLinearFitterHistoArray.Delete();
	145	fLinearFitterPArray.Delete();
	146	fLinearFitterEArray.Delete();
	147
661c7be6	148	if ( fDebugStreamer ) delete fDebugStreamer;
661c7be6	149
3a0f6479	150	}
64942b85	151	//_____________________________________________________________________________
	152	void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
	153	{
	154	//
	155	// Copy function
	156	//
	157
	158	AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
	159
	160	// Copy only the histos
	161	for (Int_t idet = 0; idet < 540; idet++){
	162	if(fLinearFitterHistoArray.UncheckedAt(idet)){
0ca1720e	163	TH2S hped1 = (TH2S )target.GetLinearFitterHisto(idet,kTRUE);
ba94744a	164	//hped1->SetDirectory(0);
0ca1720e	165	hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
64942b85	166	}
64942b85	167	}
ba94744a	168
64942b85	169	TObject::Copy(c);
	170
	171	}
	172	//_____________________________________________________________________________
6bbdc11a	173	Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list)
64942b85	174	{
	175	// Merge list of objects (needed by PROOF)
	176
	177	if (!list)
	178	return 0;
	179
	180	if (list->IsEmpty())
	181	return 1;
	182
	183	TIterator* iter = list->MakeIterator();
	184	TObject* obj = 0;
	185
	186	// collection of generated histograms
	187	Int_t count=0;
	188	while((obj = iter->Next()) != 0)
ba94744a	189	{
	190	AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
	191	if (entry == 0) continue;
	192
	193	// Copy only the histos
	194	for (Int_t idet = 0; idet < 540; idet++){
	195	if(entry->GetLinearFitterHisto(idet)){
0ca1720e	196	TH2S hped1 = (TH2S )GetLinearFitterHisto(idet,kTRUE);
ba94744a	197	Double_t entriesa = hped1->GetEntries();
0ca1720e	198	Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
ba94744a	199	if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
ba94744a	200	}
64942b85	201	}
64942b85	202
ba94744a	203	count++;
64942b85	204	}
64942b85	205
	206	return count;
	207	}
3a0f6479	208	//_____________________________________________________________________
fe4ee353	209	void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
3a0f6479	210	{
	211	//
	212	// Add histo
	213	//
	214
	215	fVersion++;
	216
	217	for (Int_t idet = 0; idet < 540; idet++){
fe4ee353	218	const TH2S hped = (TH2S)ped->GetLinearFitterHistoNoForce(idet);
3a0f6479	219	//printf("idet %d\n",idet);
	220	if ( hped != 0x0 ){
	221	//printf("add\n");
0ca1720e	222	TH2S hped1 = (TH2S )GetLinearFitterHisto(idet,kTRUE);
ba94744a	223	Double_t entriesa = hped1->GetEntries();
	224	Double_t entriesb = hped->GetEntries();
	225	if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
3a0f6479	226	}
	227	}
	228	}
	229	//______________________________________________________________________________________
0ca1720e	230	TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
3a0f6479	231	{
	232	//
	233	// return pointer to TH2F histo
	234	// if force is true create a new histo if it doesn't exist allready
	235	//
	236	if ( !force \|\| fLinearFitterHistoArray.UncheckedAt(detector) )
0ca1720e	237	return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
3a0f6479	238
fe4ee353	239	return GetLinearFitterHistoForce(detector);
	240
	241	}
	242	//______________________________________________________________________________________
	243	TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
	244	{
	245	//
	246	// return pointer to TH2F histo
	247	// if NULL create a new histo if it doesn't exist allready
	248	//
	249	if (fLinearFitterHistoArray.UncheckedAt(detector))
	250	return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
	251
	252	// if we are forced and TLinearFitter doesn't yes exist create it
	253
	254	// new TH2F
	255	TString name("LFDV");
	256	name += detector;
	257	name += "version";
	258	name += fVersion;
	259
	260	TH2S lfdv = new TH2S((const Char_t )name,(const Char_t *) name
	261	,36,-0.9,0.9,48
	262	,-1.2,1.2);
	263	lfdv->SetXTitle("tan(phi_{track})");
	264	lfdv->SetYTitle("dy/dt");
	265	lfdv->SetZTitle("Number of clusters");
	266	lfdv->SetStats(0);
	267	lfdv->SetDirectory(0);
	268
	269	fLinearFitterHistoArray.AddAt(lfdv,detector);
	270	return lfdv;
3a0f6479	271	}
	272	//______________________________________________________________________________________
	273	Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
	274	{
	275	//
	276	// return param for this detector
	277	//
	278	if ( fLinearFitterPArray.UncheckedAt(detector) ){
	279	const TVectorD vectorP = (TVectorD)fLinearFitterPArray.UncheckedAt(detector);
	280	if(!param) param = new TVectorD(2);
	281	for(Int_t k = 0; k < 2; k++){
	282	(param)[k] = (vectorP)[k];
	283	}
	284	return kTRUE;
	285	}
	286	else return kFALSE;
	287
	288	}
	289	//______________________________________________________________________________________
	290	Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
	291	{
	292	//
	293	// return error for this detector
	294	//
	295	if ( fLinearFitterEArray.UncheckedAt(detector) ){
	296	const TVectorD vectorE = (TVectorD)fLinearFitterEArray.UncheckedAt(detector);
	297	if(!error) error = new TVectorD(3);
	298	for(Int_t k = 0; k < 3; k++){
	299	(error)[k] = (vectorE)[k];
	300	}
	301	return kTRUE;
	302	}
	303	else return kFALSE;
	304
	305	}
	306	//______________________________________________________________________________________
	307	void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
	308	{
	309	//
	310	// Fill the 2D histos for debugging
	311	//
ba94744a	312
0ca1720e	313	TH2S h = ((TH2S ) GetLinearFitterHisto(detector,kTRUE));
ba94744a	314	Double_t nbentries = h->GetEntries();
0ca1720e	315	if(nbentries < 5*32767) h->Fill(tnp,pars1);
3a0f6479	316
	317	}
	318	//____________Functions fit Online CH2d________________________________________
	319	void AliTRDCalibraVdriftLinearFit::FillPEArray()
	320	{
	321	//
	322	// Fill fLinearFitterPArray and fLinearFitterEArray from inside
	323	//
	324
	325
	326	Int_t *arrayI = new Int_t[540];
	327	for(Int_t k = 0; k< 540; k++){
	328	arrayI[k] = 0;
	329	}
	330
	331	// Loop over histos
	332	for(Int_t cb = 0; cb < 540; cb++){
0ca1720e	333	const TH2S linearfitterhisto = (TH2S)fLinearFitterHistoArray.UncheckedAt(cb);
3a0f6479	334	//printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);
3a0f6479	335
0ca1720e	336	if ( linearfitterhisto != 0 ){
3a0f6479	337
	338	// Fill a linearfitter
	339	TAxis *xaxis = linearfitterhisto->GetXaxis();
	340	TAxis *yaxis = linearfitterhisto->GetYaxis();
	341	TLinearFitter linearfitter = TLinearFitter(2,"pol1");
	342	//printf("test\n");
a5dcf618	343	Double_t integral = linearfitterhisto->Integral();
	344	//printf("Integral is %f\n",integral);
	345	Bool_t securitybreaking = kFALSE;
	346	if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
3a0f6479	347	for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
3a0f6479	348	for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
8a606eb5	349	if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
3a0f6479	350	Double_t x = xaxis->GetBinCenter(ibinx+1);
3a0f6479	351	Double_t y = yaxis->GetBinCenter(ibiny+1);
a5dcf618	352
3a0f6479	353	for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
a5dcf618	354	if(!securitybreaking){
	355	linearfitter.AddPoint(&x,y);
	356	arrayI[cb]++;
	357	}
	358	else {
	359	if(arrayI[cb]< 1198){
	360	linearfitter.AddPoint(&x,y);
	361	arrayI[cb]++;
	362	}
	363	}
3a0f6479	364	}
a5dcf618	365
3a0f6479	366	}
	367	}
	368	}
	369
	370	//printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
	371
	372	// Eval the linear fitter
	373	if(arrayI[cb]>10){
	374	TVectorD *par = new TVectorD(2);
	375	TVectorD pare = TVectorD(2);
	376	TVectorD *parE = new TVectorD(3);
a5dcf618	377	//printf("Fit\n");
2a1a7b36	378	if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) \|\| ((!fRobustFit) && (linearfitter.Eval()==0))) {
a5dcf618	379	//if((linearfitter.Eval()==0)) {
a5dcf618	380	//printf("Take the param\n");
c1105918	381	linearfitter.GetParameters(*par);
a5dcf618	382	//printf("Done\n");
c1105918	383	//linearfitter.GetErrors(pare);
	384	//Float_t ppointError = TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
	385	//(parE)[0] = pare[0]ppointError;
	386	//(parE)[1] = pare[1]ppointError;
	387
	388	(*parE)[0] = 0.0;
	389	(*parE)[1] = 0.0;
	390	(*parE)[2] = (Double_t) arrayI[cb];
	391	fLinearFitterPArray.AddAt(par,cb);
	392	fLinearFitterEArray.AddAt(parE,cb);
	393
	394	//par->Print();
	395	//parE->Print();
	396	}
a5dcf618	397	//printf("Finish\n");
ba94744a	398	}
c1105918	399
07f9b7f2	400	//delete linearfitterhisto;
ba94744a	401
3a0f6479	402	}// if something
fdc15553	403
3a0f6479	404	}
a987273c	405
a987273c	406	delete [] arrayI;
3a0f6479	407
d1e1a064	408	}
c1105918	409
661c7be6	410	//____________Functions fit Online CH2d________________________________________
	411	void AliTRDCalibraVdriftLinearFit::FillPEArray2()
	412	{
	413	//
	414	// Fill fFitterPArray and fFitterEArray from inside
	415	//
	416
	417	// Loop over histos
	418	TF1 f1 = new TF1("f1","[0]+[1]x",-1,1);
	419
	420	for(Int_t cb = 0; cb < 540; cb++){
	421	//const TH2S fitterhisto = (TH2S)fLinearFitterHistoArray.UncheckedAt(cb);
	422	TH2S fitterhisto = (TH2S)fLinearFitterHistoArray.UncheckedAt(cb);
	423	//printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) fitterhisto);
	424
	425	if ( fitterhisto != 0 ){
	426
	427	Int_t nEntries=0;
	428	TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
ce42378b	429	if (!gg) continue;
661c7be6	430	// Number of points of the TGraphErrors
	431	if(gg->GetN() < 20) {
	432	if(gg) delete gg;
	433	continue;
	434	}
	435	//printf("det %d, number of points %d, nEntries %d\n",cb,gg->GetN(),nEntries);
	436	gg->Fit(f1,"Q0");
	437
	438	TVectorD *par = new TVectorD(2);
	439	TVectorD *parE = new TVectorD(3);
	440	(*parE)[0] = 0.0;
	441	(*parE)[1] = 0.0;
	442	(*parE)[2] = (Double_t) nEntries;
	443	(*par)[0] = f1->GetParameter(0);
	444	(*par)[1] = f1->GetParameter(1);
	445	fLinearFitterPArray.AddAt(par,cb);
	446	fLinearFitterEArray.AddAt(parE,cb);
	447
	448	if(fDebugLevel==0) {
	449	if(gg) delete gg;
	450	}
	451	else {
	452	if(cb==fSeeDetector) {
	453	gStyle->SetPalette(1);
	454	gStyle->SetOptStat(1111);
	455	gStyle->SetPadBorderMode(0);
	456	gStyle->SetCanvasColor(10);
	457	gStyle->SetPadLeftMargin(0.13);
	458	gStyle->SetPadRightMargin(0.01);
	459	TCanvas *stat = new TCanvas("stat","",50,50,600,800);
	460	stat->cd(1);
	461	fitterhisto->Draw("colz");
	462	gg->Draw("P");
	463	f1->Draw("same");
	464	break;
	465	}
	466	}
	467	}
	468	}
	469	if(fDebugLevel==0) delete f1;
	470	}
	471
	472	//_________Helper function__________________________________________________
	473	TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
	474	{
	475	TF1 fg("fg", "gaus", -10., 30.);
	476	TGraphErrors *gp = new TGraphErrors();
c1105918	477
661c7be6	478	TAxis *ax(h2->GetXaxis());
	479	TAxis *ay(h2->GetYaxis());
	480	TH1D *h1(NULL);
	481	for(Int_t ipt(0), jpt(1), ig(0); ipt<ax->GetNbins(); ipt++, jpt++){
	482	h1 = h2->ProjectionY("py", jpt, jpt);
	483	fg.SetParameter(1, h1->GetMean());
	484	//Float_t x(ax->GetBinCenter(jpt));
	485	Int_t n=Int_t(h1->Integral(1, h1->GetNbinsX()));
	486	nEntries+=n;
	487	if(n < 15){
	488	//Warning("drawMS()", Form("reject x[%d]=%f on n=%d", jpt, x, n));
	489	continue;
	490	}
	491	h1->Fit(&fg, "WWQ0");
	492	if(fg.GetNDF()<2){
	493	//Warning("drawMS()", Form("reject x[%d]=%f on NDF=%d", jpt, x, fg.GetNDF()));
	494	continue;
	495	}
64c7fdf8	496	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	497	gp->SetPoint(ig, ax->GetBinCenter(jpt), fg.GetParameter(1));
	498	gp->SetPointError(ig, 0, TMath::Sqrt(pow(fg.GetParError(1),2) + (1/pow(fg.GetParameter(0),2))));
	499	ig++;
	500	}
	501	delete h1;
	502	return gp;
	503	}