Update timestamps for new AMANDA simulation (17/02/2015)
[u/mrichter/AliRoot.git] / TPC / CalibMacros / CalibKalman.C
CommitLineData
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
22void 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/(24*60*60.);
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