TPC/CalibMacros/CalibKalman.C

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