[u/mrichter/AliRoot.git] / TPC / CalibMacros / CalibKalman.C

/// \file CalibKalman.C
///
/// ~~~
/// .x ~/rootlogon.C
/// gSystem->Load("libANALYSIS");
/// gSystem->Load("libTPCcalib");
/// 
/// .L $ALICE_ROOT/TPC/CalibMacros/CalibKalman.C+
/// ~~~


#include "TTreeStream.h"
#include "TFile.h"
#include "TH2D.h"
#include "THnSparse.h"
#include "AliExternalTrackParam.h"
#include "AliTPCcalibTime.h"
#include "AliTPCkalmanTime.h"


void FitKalman(Int_t maxCount=50){
  TFile fcalib("CalibObjects.root");
  TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  AliTPCcalibTime * calibTime = ( AliTPCcalibTime *)array->FindObject("calibTime");
  AliTPCkalmanTime kalman;
  //
  //

  TTreeSRedirector *pcstream = new TTreeSRedirector("debugFit.root");
  TH2D * hdz = calibTime->GetHistVdrift()->Projection(1,0);
  TH2D * hpress = calibTime->GetHistVdrift()->Projection(4,0);
  TH2D * htempA = calibTime->GetHistVdrift()->Projection(5,0);
  TH2D * htempC = calibTime->GetHistVdrift()->Projection(6,0);
  
  TH1 * hentries =  calibTime->GetHistVdrift()->Projection(0);
  Int_t nentries =hentries->GetNbinsX();
  TH1 *htemp=0;
  Int_t counter=0;
  for (Int_t ibin=0; ibin<nentries;ibin++){
    Double_t icon  = hentries->GetBinContent(ibin);
    if (icon<5) continue;
    if (counter>maxCount) break;
    Double_t time = hentries->GetBinCenter(ibin);
    //
    htemp   = hdz->ProjectionY("aaa",ibin,ibin+1);
    Double_t dzmean   = htemp->GetMean();
    Double_t dzrms    =  htemp->GetRMS();
    delete htemp;
    //
    htemp   = hpress->ProjectionY("aaa",ibin,ibin+1);
    Double_t dpress    = htemp->GetMean();
    Double_t rpress    =  htemp->GetRMS();
    delete htemp;
    htemp   = htempA->ProjectionY("aaa",ibin,ibin+1);
    Double_t mtemp    =  htemp->GetMean();
    Double_t rtemp    =  htemp->GetRMS();
    delete htemp;
    
    printf("ibin=%d\ttime=%f\tentries=%f\t%f\t%f\t%f\t%f\t%f\t%f\n",ibin,time,icon,dzmean,dzrms,dpress,rpress, mtemp,rtemp);
    if (counter==0) {
      kalman.Init(time,0,1,0.05,0.02);
    }
    Double_t sigma0t = 0.005/(24*60*60.);
    kalman.Propagate(time,sigma0t,pcstream);
    Float_t dvdriftrel = dzmean/500.;
    Float_t dverror    = (dzrms+0.1)/500.;
    Double_t  temperature = mtemp+273.15;
    Double_t  povertMeas = dpress/temperature;
    if (mtemp<15) continue;
    
    Float_t  nominalTemp = 19.03+273.15;
    Float_t  nominalPress= 973;
    Double_t povertNom = nominalPress/nominalTemp  ;
    Float_t dptratio   = (povertMeas-povertNom)/povertNom;
    kalman.Update(dvdriftrel,dverror,dptratio,pcstream);
    kalman.fState->Print();
    kalman.fCovariance->Print();
    (*pcstream)<<"fit"<<
      "time="<<time<<
      "dvdrift="<<dvdriftrel<<
      "dverror="<<dverror<<
      "dzrms="<<dzrms<<
      "dpt="<<dptratio<<
      "k.="<<&kalman<<
      "\n";
    counter++;
  }
  delete pcstream;
}
Commit	Line	Data
35d81915	1	/// \file CalibKalman.C
	2	///
	3	/// ~~~
	4	/// .x ~/rootlogon.C
	5	/// gSystem->Load("libANALYSIS");
	6	/// gSystem->Load("libTPCcalib");
	7	///
	8	/// .L $ALICE_ROOT/TPC/CalibMacros/CalibKalman.C+
	9	/// ~~~
	10
	11
258112da	12
	13	#include "TTreeStream.h"
	14	#include "TFile.h"
	15	#include "TH2D.h"
	16	#include "THnSparse.h"
	17	#include "AliExternalTrackParam.h"
	18	#include "AliTPCcalibTime.h"
	19	#include "AliTPCkalmanTime.h"
	20
	21
	22	void FitKalman(Int_t maxCount=50){
	23	TFile fcalib("CalibObjects.root");
	24	TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
	25	AliTPCcalibTime * calibTime = ( AliTPCcalibTime *)array->FindObject("calibTime");
	26	AliTPCkalmanTime kalman;
	27	//
	28	//
	29
	30	TTreeSRedirector *pcstream = new TTreeSRedirector("debugFit.root");
	31	TH2D * hdz = calibTime->GetHistVdrift()->Projection(1,0);
	32	TH2D * hpress = calibTime->GetHistVdrift()->Projection(4,0);
	33	TH2D * htempA = calibTime->GetHistVdrift()->Projection(5,0);
	34	TH2D * htempC = calibTime->GetHistVdrift()->Projection(6,0);
	35
	36	TH1 * hentries = calibTime->GetHistVdrift()->Projection(0);
	37	Int_t nentries =hentries->GetNbinsX();
	38	TH1 *htemp=0;
	39	Int_t counter=0;
	40	for (Int_t ibin=0; ibin<nentries;ibin++){
	41	Double_t icon = hentries->GetBinContent(ibin);
	42	if (icon<5) continue;
	43	if (counter>maxCount) break;
	44	Double_t time = hentries->GetBinCenter(ibin);
	45	//
	46	htemp = hdz->ProjectionY("aaa",ibin,ibin+1);
	47	Double_t dzmean = htemp->GetMean();
	48	Double_t dzrms = htemp->GetRMS();
	49	delete htemp;
	50	//
	51	htemp = hpress->ProjectionY("aaa",ibin,ibin+1);
	52	Double_t dpress = htemp->GetMean();
	53	Double_t rpress = htemp->GetRMS();
	54	delete htemp;
	55	htemp = htempA->ProjectionY("aaa",ibin,ibin+1);
	56	Double_t mtemp = htemp->GetMean();
	57	Double_t rtemp = htemp->GetRMS();
	58	delete htemp;
	59
	60	printf("ibin=%d\ttime=%f\tentries=%f\t%f\t%f\t%f\t%f\t%f\t%f\n",ibin,time,icon,dzmean,dzrms,dpress,rpress, mtemp,rtemp);
	61	if (counter==0) {
	62	kalman.Init(time,0,1,0.05,0.02);
	63	}
	64	Double_t sigma0t = 0.005/(246060.);
	65	kalman.Propagate(time,sigma0t,pcstream);
	66	Float_t dvdriftrel = dzmean/500.;
	67	Float_t dverror = (dzrms+0.1)/500.;
	68	Double_t temperature = mtemp+273.15;
	69	Double_t povertMeas = dpress/temperature;
	70	if (mtemp<15) continue;
	71
	72	Float_t nominalTemp = 19.03+273.15;
	73	Float_t nominalPress= 973;
	74	Double_t povertNom = nominalPress/nominalTemp ;
	75	Float_t dptratio = (povertMeas-povertNom)/povertNom;
76	kalman.Update(dvdriftrel,dverror,dptratio,pcstream);
77	kalman.fState->Print();
78	kalman.fCovariance->Print();
79	(*pcstream)<<"fit"<<
80	"time="<<time<<
81	"dvdrift="<<dvdriftrel<<
82	"dverror="<<dverror<<
83	"dzrms="<<dzrms<<
84	"dpt="<<dptratio<<
85	"k.="<<&kalman<<
86	"\n";
87	counter++;
88	}
89	delete pcstream;
90	}
91