[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),
  fMinNpointsFit(11),
  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),
  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),
  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
      //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();

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