ROOT logo

// See http://en.wikipedia.org/wiki/Kalman_filter
//
/*
.x ~/UliStyle.C
gSystem->AddIncludePath("-I$ALICE_ROOT/STAT");
gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
gSystem->Load("libSTAT.so");

.L $ALICE_ROOT/TPC/CalibMacros/DriftKalman.C++ 
  //
  //
  //
  TFile fgoofie("goofieTree.root");
  TTree *chainTime = (TTree*)fgoofie.Get("drift");
  chainTime->SetAlias("mcut","abs(areaF/areaN-1.117)<0.1");
  chainTime->SetAlias("atpcor","(1000*press/(3.329*(tempF)))-1");
  chainTime->SetAlias("dvdrift","(((vdrift/2.636)-1))");
  chainTime->SetAlias("trCorr","0");
 
  Int_t npoints= chainTime->Draw("time:(dvdrift)-trCorr:(atpcor)","mcut");
  Double_t *timeV = new Double_t[npoints];
  Double_t *vdriftV = new Double_t[npoints];
  Double_t *cPTV = new Double_t[npoints];
  memcpy(timeV,chainTime->GetV1(),npoints*sizeof(Double_t));
  memcpy(vdriftV,chainTime->GetV2(),npoints*sizeof(Double_t));
  memcpy(cPTV,chainTime->GetV3(),npoints*sizeof(Double_t));

  Double_t rms=0.006, corr;
  covXk(0,0)=0.01; covXk(1,1)=0.1; covXk(2,2)=0.1; covXk(1,0)=0; covXk(0,1)=0;
  vecXk(0,0)=0; vecXk(1,0)=1; vecXk(2,0)=0;          //initial stat

  recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,10/(250*24*60*60.),kFALSE);
  recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,10/(250*24*60*60.),kTRUE); 
  TFile f("realvdrift.root");
  drift->Draw("vecXk.fElements[1]:time","isOK","")
 

  //
  // Real data test usage
  //

  AliXRDPROOFtoolkit tool;
  TChain * chainTime = tool.MakeChain("time.txt","timeInfo",0,10200);
  chainTime->Lookup();
  //
  
  Double_t rms, corr;
  chainTime->SetAlias("atpcorP","press0/(3.329*(temp0+273.15))-1");
  chainTime->SetAlias("mcutP","press0>0");
  chainTime->SetAlias("ncutP","press0<1");
  chainTime->SetAlias("atpcorG","temp1.fElements[14]/(0.003378*(temp0+273.15))-1");
  chainTime->SetAlias("mcutG","temp1.fElements[14]>0");

  chainTime->SetAlias("atpcor","(mcutP*atpcorP+(ncutP)*atpcorG)");
  chainTime->SetAlias("mcut","temp1.fElements[14]>0");
  chainTime->SetAlias("dvdrift","(fDz/500)");
  
  fitLinear(chainTime, rms,corr);
 

  Int_t npoints= chainTime->Draw("time:(dvdrift)-trCorr:(atpcor)","mcut&&abs(fDz-500*driftG)<1");
  Double_t *timeV = new Double_t[npoints];
  Double_t *vdriftV = new Double_t[npoints];
  Double_t *cPTV = new Double_t[npoints];
  memcpy(timeV,chainTime->GetV1(),npoints*sizeof(Double_t));
  memcpy(vdriftV,chainTime->GetV2(),npoints*sizeof(Double_t));
  memcpy(cPTV,chainTime->GetV3(),npoints*sizeof(Double_t));
  

  kalmanTime.Init(0,0,-1,0.01,0.01);

  recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,1/(250*24*60*60.),kFALSE);
  recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,1/(250*24*60*60.),kTRUE); 

  TFile f("realvdrift.root");
  drift->Draw("k.fState.fElements[0]:time","isOK","");
*/

#include "TMatrixD.h"
#include "TRandom.h"
#include "TTreeStream.h"
#include "TMath.h"
#include "TString.h"
#include "TChain.h"
#include "TStatToolkit.h"
#include "AliTPCkalmanTime.h"


Double_t cshift  = 0;
Double_t cPT     = 0;
Double_t floatS  = 0.00005;   // 
Double_t floatTP = 0.0005;    //
Double_t noiseV  = 0.3/250.;  // cosmic drift precission
Double_t vdrift    = 0;

AliTPCkalmanTime kalmanTime;


void recVDrift(Int_t npoints, Double_t *timeV, Double_t *vdriftV, Double_t * cPTV, Double_t noise, Double_t sshift){
  //
  // chainTime->Draw("time:(dvdrift):(atpcor)","mcut&&abs(((dvdrift))-(atpcor))<0.004");
  //
  TTreeSRedirector *pcstream = new TTreeSRedirector("realvdrift.root");
  //
  Int_t nused=0;
  for (Int_t ipoint=0; ipoint<npoints;ipoint++){    
    //
    if (ipoint==0) kalmanTime.fTime=timeV[ipoint];
    TMatrixD &vecXk = *(kalmanTime.fState);
    TMatrixD &covXk = *(kalmanTime.fCovariance);
    Double_t fdrift = vecXk(0,0)+ cPTV[ipoint]*vecXk(1,0);
    Double_t ddrift = fdrift-vdriftV[ipoint];
    Double_t sdrift = TMath::Sqrt(covXk(0,0)+noise*noise);
    Bool_t isOK = nused<10||TMath::Abs(ddrift)<5*sdrift;
    kalmanTime.Propagate(timeV[ipoint],sshift);
    kalmanTime.Update(vdriftV[ipoint],noise,cPTV[ipoint],pcstream);
    nused++;
    (*pcstream)<<"drift"<<
      "ipoint="<<ipoint<<      
      "time="<<timeV[ipoint]<<
      "isOK="<<isOK<<
      "noise="<<noise<<
      "cPT="<<cPTV[ipoint]<<
      "vdrift="<<vdriftV[ipoint]<<
      "fdrift="<<fdrift<<
      "k.="<<&kalmanTime<<
      "\n";
  }
  delete pcstream;
}


void recVDriftSort(Int_t npoints, Double_t *timeV, Double_t *vdriftV, Double_t * cPTV, Double_t noise, Double_t sshift, Bool_t sort){
  //
  // sort entries before fit
  //
  Int_t *index = new Int_t[npoints];
  Double_t *tmpsort = new Double_t[npoints];
  TMath::Sort(npoints, timeV, index,sort);
  //
  for (Int_t i=0; i<npoints;i++){
    tmpsort[i]= timeV[index[i]];
  }
  for (Int_t i=0; i<npoints;i++){
    timeV[i] =tmpsort[i];
  }
  //
  for (Int_t i=0; i<npoints;i++){
    tmpsort[i]= vdriftV[index[i]];
  }
  for (Int_t i=0; i<npoints;i++){
    vdriftV[i] =tmpsort[i];
  }
  //
  for (Int_t i=0; i<npoints;i++){
    tmpsort[i]= cPTV[index[i]];
  }
  for (Int_t i=0; i<npoints;i++){
    cPTV[i] =tmpsort[i];
  }
  delete [] index;
  delete [] tmpsort;
  recVDrift(npoints,timeV,vdriftV,cPTV,noise,sshift);
}


void  fitLinear(TChain * chainTime, Double_t& rms, Double_t & corr){
  //
  // Fit globaly - make aliases;
  //
  TStatToolkit toolkit;
  Double_t chi2=0;
  Int_t    npoints=0;
  TVectorD fitParam;
  TMatrixD covMatrix;
  TString  fitStr="";  
  fitStr+="(atpcor)++";
  fitStr+="(trigger==1)++";
  fitStr+="(trigger==2)++";
  fitStr+="(trigger==4)++";
  fitStr+="(trigger==8)++";

  TString *strvdpt1 = 0;
  strvdpt1 = toolkit.FitPlane(chainTime,"(dvdrift)",fitStr.Data(),"mcut" , chi2,npoints,fitParam,covMatrix); 
  chainTime->SetAlias("driftG",strvdpt1->Data());

  strvdpt1 =toolkit.FitPlane(chainTime,"(dvdrift)",fitStr.Data(),"mcut&&abs(fDz-500*driftG)<2" , chi2,npoints,fitParam,covMatrix);    
  chainTime->SetAlias("driftG",strvdpt1->Data());
  
 
  rms = TMath::Sqrt(chi2/npoints);
  corr = fitParam[1];  // cPT parameter
  //
  //
  fitStr="";
  fitStr+="(trigger==1)++";
  fitStr+="(trigger==2)++";
  fitStr+="(trigger==4)++";
  fitStr+="(trigger==8)++";
  TString fitExpr = "(dvdrift)-";
  fitExpr+=fitParam[1];
  fitExpr+="* (atpcor)";
  TString *strCorr = toolkit.FitPlane(chainTime,fitExpr.Data(),fitStr.Data(),"abs(fDz-500*driftG)<1" , chi2,npoints,fitParam,covMatrix);  
  chainTime->SetAlias("trCorr",strCorr->Data());


  fitStr="";  
  fitStr+="(atpcor)++";
  fitStr+="(atpcor)^2++";
  fitStr+="(trigger==1)++";
  fitStr+="(trigger==2)++";
  fitStr+="(trigger==4)++";
  fitStr+="(trigger==8)++";
  TString *strvdpt2 = 0;
  strvdpt2 =toolkit.FitPlane(chainTime,"(dvdrift)",fitStr.Data(),"mcut&&abs(fDz-500*driftG)<2" , chi2,npoints,fitParam,covMatrix);    
  chainTime->SetAlias("driftG2",strvdpt2->Data());


}
 DriftKalman.C:1
 DriftKalman.C:2
 DriftKalman.C:3
 DriftKalman.C:4
 DriftKalman.C:5
 DriftKalman.C:6
 DriftKalman.C:7
 DriftKalman.C:8
 DriftKalman.C:9
 DriftKalman.C:10
 DriftKalman.C:11
 DriftKalman.C:12
 DriftKalman.C:13
 DriftKalman.C:14
 DriftKalman.C:15
 DriftKalman.C:16
 DriftKalman.C:17
 DriftKalman.C:18
 DriftKalman.C:19
 DriftKalman.C:20
 DriftKalman.C:21
 DriftKalman.C:22
 DriftKalman.C:23
 DriftKalman.C:24
 DriftKalman.C:25
 DriftKalman.C:26
 DriftKalman.C:27
 DriftKalman.C:28
 DriftKalman.C:29
 DriftKalman.C:30
 DriftKalman.C:31
 DriftKalman.C:32
 DriftKalman.C:33
 DriftKalman.C:34
 DriftKalman.C:35
 DriftKalman.C:36
 DriftKalman.C:37
 DriftKalman.C:38
 DriftKalman.C:39
 DriftKalman.C:40
 DriftKalman.C:41
 DriftKalman.C:42
 DriftKalman.C:43
 DriftKalman.C:44
 DriftKalman.C:45
 DriftKalman.C:46
 DriftKalman.C:47
 DriftKalman.C:48
 DriftKalman.C:49
 DriftKalman.C:50
 DriftKalman.C:51
 DriftKalman.C:52
 DriftKalman.C:53
 DriftKalman.C:54
 DriftKalman.C:55
 DriftKalman.C:56
 DriftKalman.C:57
 DriftKalman.C:58
 DriftKalman.C:59
 DriftKalman.C:60
 DriftKalman.C:61
 DriftKalman.C:62
 DriftKalman.C:63
 DriftKalman.C:64
 DriftKalman.C:65
 DriftKalman.C:66
 DriftKalman.C:67
 DriftKalman.C:68
 DriftKalman.C:69
 DriftKalman.C:70
 DriftKalman.C:71
 DriftKalman.C:72
 DriftKalman.C:73
 DriftKalman.C:74
 DriftKalman.C:75
 DriftKalman.C:76
 DriftKalman.C:77
 DriftKalman.C:78
 DriftKalman.C:79
 DriftKalman.C:80
 DriftKalman.C:81
 DriftKalman.C:82
 DriftKalman.C:83
 DriftKalman.C:84
 DriftKalman.C:85
 DriftKalman.C:86
 DriftKalman.C:87
 DriftKalman.C:88
 DriftKalman.C:89
 DriftKalman.C:90
 DriftKalman.C:91
 DriftKalman.C:92
 DriftKalman.C:93
 DriftKalman.C:94
 DriftKalman.C:95
 DriftKalman.C:96
 DriftKalman.C:97
 DriftKalman.C:98
 DriftKalman.C:99
 DriftKalman.C:100
 DriftKalman.C:101
 DriftKalman.C:102
 DriftKalman.C:103
 DriftKalman.C:104
 DriftKalman.C:105
 DriftKalman.C:106
 DriftKalman.C:107
 DriftKalman.C:108
 DriftKalman.C:109
 DriftKalman.C:110
 DriftKalman.C:111
 DriftKalman.C:112
 DriftKalman.C:113
 DriftKalman.C:114
 DriftKalman.C:115
 DriftKalman.C:116
 DriftKalman.C:117
 DriftKalman.C:118
 DriftKalman.C:119
 DriftKalman.C:120
 DriftKalman.C:121
 DriftKalman.C:122
 DriftKalman.C:123
 DriftKalman.C:124
 DriftKalman.C:125
 DriftKalman.C:126
 DriftKalman.C:127
 DriftKalman.C:128
 DriftKalman.C:129
 DriftKalman.C:130
 DriftKalman.C:131
 DriftKalman.C:132
 DriftKalman.C:133
 DriftKalman.C:134
 DriftKalman.C:135
 DriftKalman.C:136
 DriftKalman.C:137
 DriftKalman.C:138
 DriftKalman.C:139
 DriftKalman.C:140
 DriftKalman.C:141
 DriftKalman.C:142
 DriftKalman.C:143
 DriftKalman.C:144
 DriftKalman.C:145
 DriftKalman.C:146
 DriftKalman.C:147
 DriftKalman.C:148
 DriftKalman.C:149
 DriftKalman.C:150
 DriftKalman.C:151
 DriftKalman.C:152
 DriftKalman.C:153
 DriftKalman.C:154
 DriftKalman.C:155
 DriftKalman.C:156
 DriftKalman.C:157
 DriftKalman.C:158
 DriftKalman.C:159
 DriftKalman.C:160
 DriftKalman.C:161
 DriftKalman.C:162
 DriftKalman.C:163
 DriftKalman.C:164
 DriftKalman.C:165
 DriftKalman.C:166
 DriftKalman.C:167
 DriftKalman.C:168
 DriftKalman.C:169
 DriftKalman.C:170
 DriftKalman.C:171
 DriftKalman.C:172
 DriftKalman.C:173
 DriftKalman.C:174
 DriftKalman.C:175
 DriftKalman.C:176
 DriftKalman.C:177
 DriftKalman.C:178
 DriftKalman.C:179
 DriftKalman.C:180
 DriftKalman.C:181
 DriftKalman.C:182
 DriftKalman.C:183
 DriftKalman.C:184
 DriftKalman.C:185
 DriftKalman.C:186
 DriftKalman.C:187
 DriftKalman.C:188
 DriftKalman.C:189
 DriftKalman.C:190
 DriftKalman.C:191
 DriftKalman.C:192
 DriftKalman.C:193
 DriftKalman.C:194
 DriftKalman.C:195
 DriftKalman.C:196
 DriftKalman.C:197
 DriftKalman.C:198
 DriftKalman.C:199
 DriftKalman.C:200
 DriftKalman.C:201
 DriftKalman.C:202
 DriftKalman.C:203
 DriftKalman.C:204
 DriftKalman.C:205
 DriftKalman.C:206
 DriftKalman.C:207
 DriftKalman.C:208
 DriftKalman.C:209
 DriftKalman.C:210
 DriftKalman.C:211
 DriftKalman.C:212
 DriftKalman.C:213
 DriftKalman.C:214
 DriftKalman.C:215
 DriftKalman.C:216
 DriftKalman.C:217
 DriftKalman.C:218
 DriftKalman.C:219
 DriftKalman.C:220
 DriftKalman.C:221
 DriftKalman.C:222
 DriftKalman.C:223