[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraExbAltFit.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: AliTRDCalibraExbAltFit.cxx 46327 2011-01-10 13:29:56Z cblume $ */

////////////////////////////////////////////////////////////////////////////
//                                                                        //
// AliTRDCalibraExbAltFit                                                 //
//                                                                        //
// Does the ExB calibration by applying a quadratic 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 <TDirectory.h>
#include <TTreeStream.h>
#include <TGraphErrors.h>
#include <TF1.h>

//header file
#include "AliTRDCalibraExbAltFit.h"

ClassImp(AliTRDCalibraExbAltFit) /*FOLD00*/

//_____________________________________________________________________
AliTRDCalibraExbAltFit::AliTRDCalibraExbAltFit() : /*FOLD00*/
  TObject(),
  fVersion(0),
  fFitterHistoArray(540),
  fFitterPArray(540),
  fFitterEArray(540),
  fRobustFit(kFALSE),
  fDebugStreamer(0x0),
  fDebugLevel(0)
{
  //
  // default constructor
  //
}
//_____________________________________________________________________
AliTRDCalibraExbAltFit::AliTRDCalibraExbAltFit(const AliTRDCalibraExbAltFit &ped) : /*FOLD00*/
  TObject(ped),
  fVersion(ped.fVersion),
  fFitterHistoArray(540),
  fFitterPArray(540),
  fFitterEArray(540),
  fRobustFit(kFALSE),
  fDebugStreamer(0x0),
  fDebugLevel(0)
{
    //
    // copy constructor
    //
  for (Int_t idet = 0; idet < 540; idet++){
   
    const TVectorD     *vectorE     = (TVectorD*)ped.fFitterEArray.UncheckedAt(idet);
    const TVectorD     *vectorP     = (TVectorD*)ped.fFitterPArray.UncheckedAt(idet);
    const TH2S         *hped        = (TH2S*)ped.fFitterHistoArray.UncheckedAt(idet);
    
    if ( vectorE != 0x0 ) fFitterEArray.AddAt(new TVectorD(*vectorE), idet);
    if ( vectorP != 0x0 ) fFitterPArray.AddAt(new TVectorD(*vectorP), idet);
    if ( hped != 0x0 ){
      TH2S *hNew = (TH2S *)hped->Clone();
      //hNew->SetDirectory(0);
      fFitterHistoArray.AddAt(hNew,idet);
    }
  }
}
//_____________________________________________________________________
AliTRDCalibraExbAltFit::AliTRDCalibraExbAltFit(const TObjArray &obja) : /*FOLD00*/
  TObject(),
  fVersion(0),
  fFitterHistoArray(540),
  fFitterPArray(540),
  fFitterEArray(540),
  fRobustFit(kFALSE),
  fDebugStreamer(0x0),
  fDebugLevel(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);
      fFitterHistoArray.AddAt(hNew,idet);
    }
  }
}
//_____________________________________________________________________
AliTRDCalibraExbAltFit& AliTRDCalibraExbAltFit::operator = (const  AliTRDCalibraExbAltFit &source)
{
  //
  // assignment operator
  //
  if (&source == this) return *this;
  new (this) AliTRDCalibraExbAltFit(source);

  return *this;
}
//_____________________________________________________________________
AliTRDCalibraExbAltFit::~AliTRDCalibraExbAltFit() /*FOLD00*/
{
  //
  // destructor
  //
  fFitterHistoArray.SetOwner();
  fFitterPArray.SetOwner();
  fFitterEArray.SetOwner();

  fFitterHistoArray.Delete();
  fFitterPArray.Delete();
  fFitterEArray.Delete();

  if ( fDebugStreamer ) delete fDebugStreamer;

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

  AliTRDCalibraExbAltFit& target = (AliTRDCalibraExbAltFit &) c;

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

}
//_____________________________________________________________________________
Long64_t AliTRDCalibraExbAltFit::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) 
    {
      AliTRDCalibraExbAltFit* entry = dynamic_cast<AliTRDCalibraExbAltFit*>(obj);
      if (entry == 0) continue; 
      
      // Copy only the histos
      for (Int_t idet = 0; idet < 540; idet++){
	if(entry->GetFitterHisto(idet)){
	  TH2S *hped1 = (TH2S *)GetFitterHisto(idet,kTRUE);
	  Double_t entriesa = hped1->GetEntries();
	  Double_t entriesb = ((TH2S *)entry->GetFitterHisto(idet))->GetEntries();
	  if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetFitterHisto(idet));
	}
      }
      
      count++;
    }
  
  return count;
}
//_____________________________________________________________________
void AliTRDCalibraExbAltFit::Add(const AliTRDCalibraExbAltFit *ped)
{
  //
  // Add histo
  //

  fVersion++;

  for (Int_t idet = 0; idet < 540; idet++){
    const TH2S         *hped        = (TH2S*)ped->GetFitterHistoNoForce(idet);
    //printf("idet %d\n",idet);
    if ( hped != 0x0 ){
      //printf("add\n");
      TH2S *hped1 = (TH2S *)GetFitterHisto(idet,kTRUE);
      Double_t entriesa = hped1->GetEntries();
      Double_t entriesb = hped->GetEntries();
      if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
    }
  }
}
//______________________________________________________________________________________
TH2S* AliTRDCalibraExbAltFit::GetFitterHisto(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 || fFitterHistoArray.UncheckedAt(detector) )
	return (TH2S*)fFitterHistoArray.UncheckedAt(detector);

    return GetFitterHistoForce(detector);

}
//______________________________________________________________________________________
TH2S* AliTRDCalibraExbAltFit::GetFitterHistoForce(Int_t detector)
{
  //
  // return pointer to TH2F histo 
  // if NULL create a new histo if it doesn't exist allready
  //
  if (fFitterHistoArray.UncheckedAt(detector))
    return (TH2S*)fFitterHistoArray.UncheckedAt(detector);
  
  // if we are forced and TLinearFitter doesn't yes exist create it
  
  // new TH2F
  TString name("LFEXB");
  name += detector;
  name += "version";
  name +=  fVersion;
  
  TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name,
			30, -TMath::DegToRad()*45, TMath::DegToRad()*45, 
			30, 0.3, 1.4);
  lfdv->SetXTitle("tan(phi_{track})");
  lfdv->SetYTitle("rms");
  lfdv->SetZTitle("Number of tracklets");
  lfdv->SetStats(0);
  lfdv->SetDirectory(0);
  
  fFitterHistoArray.AddAt(lfdv,detector);
  return lfdv;
}
//______________________________________________________________________________________
Bool_t AliTRDCalibraExbAltFit::GetParam(Int_t detector, TVectorD *param)
{
    //
    // return param for this detector
    //
  if ( fFitterPArray.UncheckedAt(detector) ){
    const TVectorD     *vectorP     = (TVectorD*)fFitterPArray.UncheckedAt(detector);
    if(!param) param = new TVectorD(vectorP->GetNoElements());
    for(Int_t k = 0; k < vectorP->GetNoElements(); k++){
      (*param)[k] = (*vectorP)[k];
    }
    return kTRUE;
  }
  else return kFALSE;

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

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

}
//____________Functions fit Online CH2d________________________________________
void AliTRDCalibraExbAltFit::FillPEArray()
{
  //
  // Fill fFitterPArray and fFitterEArray 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 *fitterhisto = (TH2S*)fFitterHistoArray.UncheckedAt(cb);
    //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) fitterhisto);    

    if ( fitterhisto != 0 ){
      
      // Fill a fitter
      TAxis *xaxis = fitterhisto->GetXaxis();
      TAxis *yaxis = fitterhisto->GetYaxis();
      TLinearFitter fitter = TLinearFitter(3,"pol2");
      //printf("test\n");
      Double_t integral = fitterhisto->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 < fitterhisto->GetNbinsX(); ibinx++){
	for(Int_t ibiny = 0; ibiny < fitterhisto->GetNbinsY(); ibiny++){
	  if(fitterhisto->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)fitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
	      if(!securitybreaking){
		fitter.AddPoint(&x,y);
		arrayI[cb]++;
	      }
	      else {
		if(arrayI[cb]< 1198){
		  fitter.AddPoint(&x,y);
		  arrayI[cb]++; 
		}
	      }
	    }
	    
	  }
	}
      }
      
      //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);

      // Eval the  fitter
      if(arrayI[cb]>15){
	TVectorD  *par  = new TVectorD(3);
	TVectorD  pare = TVectorD(2);
	TVectorD  *parE = new TVectorD(3);
	//printf("Fit\n");
	//if((fitter.EvalRobust(0.8)==0)) {
	//if((fitter.EvalRobust()==0)) {
	if(((fRobustFit) && (fitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (fitter.Eval()==0))) {
	  //if((fitter.Eval()==0)) {
	  //printf("Take the param\n");
	  fitter.GetParameters(*par);
	  //printf("Done\n");
	  //fitter.GetErrors(*pare);
	  //Float_t  ppointError =  TMath::Sqrt(TMath::Abs(fitter.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];
	  fFitterPArray.AddAt(par,cb);
	  fFitterEArray.AddAt(parE,cb);
	  
	  //par->Print();
	  //parE->Print();
	}
	//printf("Finish\n");
      }
      
      //delete fitterhisto;
      
    }// if something

  }

  delete [] arrayI;
   
}

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

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

    if ( fitterhisto != 0 ){
      TF1 *f1 = new TF1("f1","[0]*(x-[1])**2+[2]",-1,1);
      
      Int_t nEntries=0;
      TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
      if (!gg) continue;
      // printf("N: %i\n",gg->GetN());
      // printf("entries: %i\n",nEntries);

      if(gg->GetN() < 20) {
	if(gg) delete gg;
	continue;	
      }
      gg->Fit(f1,"Q0");
      
      TVectorD  *par  = new TVectorD(3);
      TVectorD  *parE = new TVectorD(3);
      (*parE)[0] = 0.0;
      (*parE)[1] = 0.0;
      (*parE)[2] = (Double_t) nEntries;
      (*par)[0] = f1->GetParameter(0)*TMath::Power(f1->GetParameter(1),2)+f1->GetParameter(2);
      (*par)[1] = -2*f1->GetParameter(0)*f1->GetParameter(1);
      (*par)[2] = f1->GetParameter(0);
      fFitterPArray.AddAt(par,cb);
      fFitterEArray.AddAt(parE,cb);

      // if ( !fDebugStreamer ) {
      // 	//debug stream
      // 	TDirectory *backup = gDirectory;
      // 	fDebugStreamer = new TTreeSRedirector("TRDdebugCalibraFill.root");
      // 	if ( backup ) backup->cd();  //we don't want to be cd'd to the debug streamer
      // } 
    
      // if(cb==513){
      // 	fitterhisto->Draw();
      // 	gg->Draw("P");
      // 	f1->Draw("same");
      // }
      // double exb=f1->GetParameter(1);
      // (* fDebugStreamer) << "FindP1TrackPHtrackletV1"<<
      // 	//"snpright="<<snpright<<
      // 	"det="<<cb<<
      // 	"h2="<<(TObject*)fitterhisto<<
      // 	"gg="<<gg<<
      // 	"f1="<<f1<<
      // 	"exb="<<exb<<
      // 	"\n";

      //if(cb!=513){
      delete gg; 
      delete f1;
      //}
    }
  }
  //if(fDebugStreamer) delete fDebugStreamer;
}

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