TPC/CalibMacros/DriftKalman.C

   1
   2 /// \file DriftKalman.C
   3 /// See http://en.wikipedia.org/wiki/Kalman_filter
   4
   5 /*
   6 .x ~/UliStyle.C
   7 gSystem->AddIncludePath("-I$ALICE_ROOT/STAT");
   8 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
   9 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
  10 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  11 gSystem->Load("libSTAT");
  12
  13 .L $ALICE_ROOT/TPC/CalibMacros/DriftKalman.C++
  14   //
  15   //
  16   //
  17   TFile fgoofie("goofieTree.root");
  18   TTree *chainTime = (TTree*)fgoofie.Get("drift");
  19   chainTime->SetAlias("mcut","abs(areaF/areaN-1.117)<0.1");
  20   chainTime->SetAlias("atpcor","(1000*press/(3.329*(tempF)))-1");
  21   chainTime->SetAlias("dvdrift","(((vdrift/2.636)-1))");
  22   chainTime->SetAlias("trCorr","0");
  23
  24   Int_t npoints= chainTime->Draw("time:(dvdrift)-trCorr:(atpcor)","mcut");
  25   Double_t *timeV = new Double_t[npoints];
  26   Double_t *vdriftV = new Double_t[npoints];
  27   Double_t *cPTV = new Double_t[npoints];
  28   memcpy(timeV,chainTime->GetV1(),npoints*sizeof(Double_t));
  29   memcpy(vdriftV,chainTime->GetV2(),npoints*sizeof(Double_t));
  30   memcpy(cPTV,chainTime->GetV3(),npoints*sizeof(Double_t));
  31
  32   Double_t rms=0.006, corr;
  33   covXk(0,0)=0.01; covXk(1,1)=0.1; covXk(2,2)=0.1; covXk(1,0)=0; covXk(0,1)=0;
  34   vecXk(0,0)=0; vecXk(1,0)=1; vecXk(2,0)=0;          //initial stat
  35
  36   recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,10/(250*24*60*60.),kFALSE);
  37   recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,10/(250*24*60*60.),kTRUE);
  38   TFile f("realvdrift.root");
  39   drift->Draw("vecXk.fElements[1]:time","isOK","")
  40
  41
  42   //
  43   // Real data test usage
  44   //
  45
  46   AliXRDPROOFtoolkit tool;
  47   TChain * chainTime = tool.MakeChain("time.txt","timeInfo",0,10200);
  48   chainTime->Lookup();
  49   //
  50
  51   Double_t rms, corr;
  52   chainTime->SetAlias("atpcorP","press0/(3.329*(temp0+273.15))-1");
  53   chainTime->SetAlias("mcutP","press0>0");
  54   chainTime->SetAlias("ncutP","press0<1");
  55   chainTime->SetAlias("atpcorG","temp1.fElements[14]/(0.003378*(temp0+273.15))-1");
  56   chainTime->SetAlias("mcutG","temp1.fElements[14]>0");
  57
  58   chainTime->SetAlias("atpcor","(mcutP*atpcorP+(ncutP)*atpcorG)");
  59   chainTime->SetAlias("mcut","temp1.fElements[14]>0");
  60   chainTime->SetAlias("dvdrift","(fDz/500)");
  61
  62   fitLinear(chainTime, rms,corr);
  63
  64
  65   Int_t npoints= chainTime->Draw("time:(dvdrift)-trCorr:(atpcor)","mcut&&abs(fDz-500*driftG)<1");
  66   Double_t *timeV = new Double_t[npoints];
  67   Double_t *vdriftV = new Double_t[npoints];
  68   Double_t *cPTV = new Double_t[npoints];
  69   memcpy(timeV,chainTime->GetV1(),npoints*sizeof(Double_t));
  70   memcpy(vdriftV,chainTime->GetV2(),npoints*sizeof(Double_t));
  71   memcpy(cPTV,chainTime->GetV3(),npoints*sizeof(Double_t));
  72
  73
  74   kalmanTime.Init(0,0,-1,0.01,0.01);
  75
  76   recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,1/(250*24*60*60.),kFALSE);
  77   recVDriftSort(npoints,timeV,vdriftV,cPTV,rms,1/(250*24*60*60.),kTRUE);
  78
  79   TFile f("realvdrift.root");
  80   drift->Draw("k.fState.fElements[0]:time","isOK","");
  81 */
  82
  83 #include "TMatrixD.h"
  84 #include "TRandom.h"
  85 #include "TTreeStream.h"
  86 #include "TMath.h"
  87 #include "TString.h"
  88 #include "TChain.h"
  89 #include "TStatToolkit.h"
  90 #include "AliTPCkalmanTime.h"
  91
  92
  93 Double_t cshift  = 0;
  94 Double_t cPT     = 0;
  95 Double_t floatS  = 0.00005;   //
  96 Double_t floatTP = 0.0005;    //
  97 Double_t noiseV  = 0.3/250.;  // cosmic drift precission
  98 Double_t vdrift    = 0;
  99
 100 AliTPCkalmanTime kalmanTime;
 101
 102
 103 void recVDrift(Int_t npoints, Double_t *timeV, Double_t *vdriftV, Double_t * cPTV, Double_t noise, Double_t sshift){
 104   //
 105   // chainTime->Draw("time:(dvdrift):(atpcor)","mcut&&abs(((dvdrift))-(atpcor))<0.004");
 106   //
 107   TTreeSRedirector *pcstream = new TTreeSRedirector("realvdrift.root");
 108   //
 109   Int_t nused=0;
 110   for (Int_t ipoint=0; ipoint<npoints;ipoint++){
 111     //
 112     if (ipoint==0) kalmanTime.fTime=timeV[ipoint];
 113     TMatrixD &vecXk = *(kalmanTime.fState);
 114     TMatrixD &covXk = *(kalmanTime.fCovariance);
 115     Double_t fdrift = vecXk(0,0)+ cPTV[ipoint]*vecXk(1,0);
 116     Double_t ddrift = fdrift-vdriftV[ipoint];
 117     Double_t sdrift = TMath::Sqrt(covXk(0,0)+noise*noise);
 118     Bool_t isOK = nused<10||TMath::Abs(ddrift)<5*sdrift;
 119     kalmanTime.Propagate(timeV[ipoint],sshift);
 120     kalmanTime.Update(vdriftV[ipoint],noise,cPTV[ipoint],pcstream);
 121     nused++;
 122     (*pcstream)<<"drift"<<
 123       "ipoint="<<ipoint<<
 124       "time="<<timeV[ipoint]<<
 125       "isOK="<<isOK<<
 126       "noise="<<noise<<
 127       "cPT="<<cPTV[ipoint]<<
 128       "vdrift="<<vdriftV[ipoint]<<
 129       "fdrift="<<fdrift<<
 130       "k.="<<&kalmanTime<<
 131       "\n";
 132   }
 133   delete pcstream;
 134 }
 135
 136
 137 void recVDriftSort(Int_t npoints, Double_t *timeV, Double_t *vdriftV, Double_t * cPTV, Double_t noise, Double_t sshift, Bool_t sort){
 138   //
 139   // sort entries before fit
 140   //
 141   Int_t *index = new Int_t[npoints];
 142   Double_t *tmpsort = new Double_t[npoints];
 143   TMath::Sort(npoints, timeV, index,sort);
 144   //
 145   for (Int_t i=0; i<npoints;i++){
 146     tmpsort[i]= timeV[index[i]];
 147   }
 148   for (Int_t i=0; i<npoints;i++){
 149     timeV[i] =tmpsort[i];
 150   }
 151   //
 152   for (Int_t i=0; i<npoints;i++){
 153     tmpsort[i]= vdriftV[index[i]];
 154   }
 155   for (Int_t i=0; i<npoints;i++){
 156     vdriftV[i] =tmpsort[i];
 157   }
 158   //
 159   for (Int_t i=0; i<npoints;i++){
 160     tmpsort[i]= cPTV[index[i]];
 161   }
 162   for (Int_t i=0; i<npoints;i++){
 163     cPTV[i] =tmpsort[i];
 164   }
 165   delete [] index;
 166   delete [] tmpsort;
 167   recVDrift(npoints,timeV,vdriftV,cPTV,noise,sshift);
 168 }
 169
 170
 171 void  fitLinear(TChain * chainTime, Double_t& rms, Double_t & corr){
 172   //
 173   // Fit globaly - make aliases;
 174   //
 175   TStatToolkit toolkit;
 176   Double_t chi2=0;
 177   Int_t    npoints=0;
 178   TVectorD fitParam;
 179   TMatrixD covMatrix;
 180   TString  fitStr="";
 181   fitStr+="(atpcor)++";
 182   fitStr+="(trigger==1)++";
 183   fitStr+="(trigger==2)++";
 184   fitStr+="(trigger==4)++";
 185   fitStr+="(trigger==8)++";
 186
 187   TString *strvdpt1 = 0;
 188   strvdpt1 = toolkit.FitPlane(chainTime,"(dvdrift)",fitStr.Data(),"mcut" , chi2,npoints,fitParam,covMatrix);
 189   chainTime->SetAlias("driftG",strvdpt1->Data());
 190
 191   strvdpt1 =toolkit.FitPlane(chainTime,"(dvdrift)",fitStr.Data(),"mcut&&abs(fDz-500*driftG)<2" , chi2,npoints,fitParam,covMatrix);
 192   chainTime->SetAlias("driftG",strvdpt1->Data());
 193
 194
 195   rms = TMath::Sqrt(chi2/npoints);
 196   corr = fitParam[1];  // cPT parameter
 197   //
 198   //
 199   fitStr="";
 200   fitStr+="(trigger==1)++";
 201   fitStr+="(trigger==2)++";
 202   fitStr+="(trigger==4)++";
 203   fitStr+="(trigger==8)++";
 204   TString fitExpr = "(dvdrift)-";
 205   fitExpr+=fitParam[1];
 206   fitExpr+="* (atpcor)";
 207   TString *strCorr = toolkit.FitPlane(chainTime,fitExpr.Data(),fitStr.Data(),"abs(fDz-500*driftG)<1" , chi2,npoints,fitParam,covMatrix);
 208   chainTime->SetAlias("trCorr",strCorr->Data());
 209
 210
 211   fitStr="";
 212   fitStr+="(atpcor)++";
 213   fitStr+="(atpcor)^2++";
 214   fitStr+="(trigger==1)++";
 215   fitStr+="(trigger==2)++";
 216   fitStr+="(trigger==4)++";
 217   fitStr+="(trigger==8)++";
 218   TString *strvdpt2 = 0;
 219   strvdpt2 =toolkit.FitPlane(chainTime,"(dvdrift)",fitStr.Data(),"mcut&&abs(fDz-500*driftG)<2" , chi2,npoints,fitParam,covMatrix);
 220   chainTime->SetAlias("driftG2",strvdpt2->Data());
 221
 222
 223 }