TPC/CalibMacros/CalibAlignKalman.C

   1 /// \file CalibAlignKalman.C
   2 ///
   3 /// ~~~{.cpp}
   4 /// gSystem->AddIncludePath("-I$ALICE_ROOT/TPC");
   5 /// gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
   6 /// gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
   7 ///
   8 /// gROOT->LoadMacro("$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C+");
   9 ///
  10 /// AliTPCTransformation::BuildBasicFormulas();
  11 ///
  12 /// AliXRDPROOFtoolkit tool;
  13 /// chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0,50000);
  14 /// chainPoints->Lookup();
  15 ///
  16 /// chainMS = tool.MakeChainRandom("kalmanFit.list","kf",0,50000);
  17 /// chainMS->Lookup();
  18 ///
  19 /// chainFP = tool.MakeChainRandom("kalmanFit.list","filter",0,50000);
  20 /// chainFP->Lookup();
  21 ///
  22 /// CalibAlignKalmanFit(40000,1);
  23 /// kalmanFit0->DumpCorelation(0.8);
  24 /// TFile f("kalmanfitTPC.root");
  25 /// ~~~
  26
  27
  28 #ifdef __CINT__
  29 #else
  30 #include <fstream>
  31 #include "TSystem.h"
  32 #include "TROOT.h"
  33 #include "TRandom.h"
  34 #include "TMath.h"
  35 #include "TBits.h"
  36 #include "TFormula.h"
  37 #include "TF1.h"
  38 #include "TLinearFitter.h"
  39 #include "TFile.h"
  40 #include "TChain.h"
  41 #include "TCut.h"
  42 #include "TEntryList.h"
  43 #include "TH1F.h"
  44 #include "THnSparse.h"
  45
  46
  47 #include "AliSysInfo.h"
  48 #include "AliExternalTrackParam.h"
  49 #include "TTreeStream.h"
  50 #include "AliTrackPointArray.h"
  51 #include "AliLog.h"
  52 #include "AliTPCTransformation.h"
  53 #include "AliTPCkalmanFit.h"
  54 #include "AliMathBase.h"
  55 #include "AliXRDPROOFtoolkit.h"
  56 #endif
  57
  58
  59 //
  60 //
  61 // track point cuts
  62 //
  63 const Float_t krmsYcut      = 0.2;   // cluster RMS cut in Y - residuals form local  fit
  64 const Float_t krmsZcut      = 0.2;   // cluster RMS cut in Z - residuals from local  fit
  65 const Float_t krmsYcutGlobal= 0.2;   // cluster RMS cut in Y - residuals form global fit
  66 const Float_t krmsZcutGlobal= 2.0;   // cluster RMS cut in Z - residuals from global fit
  67 const Float_t kSigmaCut     = 5.;    // clusters to be removed
  68 const Int_t   knclCut       =  80;   // minimal number of clusters
  69 const Double_t kArmCut      = 50.;   // minimal level arm
  70 const Double_t kNsigma      = 5.;    // cut on track level
  71 //
  72 // mult. scattering cuts
  73 /*
  74 TCut cutClY("rms09.fElements[0]<0.15");
  75 TCut cuttY("rms09.fElements[1]/rms09.fElements[0]<0.9");
  76 TCut cutkY("rms09.fElements[2]<0.015");
  77 TCut cutClZ("rms09.fElements[3]<0.2");
  78 TCut cuttZ("rms09.fElements[4]/rms09.fElements[3]<0.9");
  79 TCut cutkZ("rms09.fElements[5]<0.015");
  80 TCut cutMS=cutClY+cuttY+cutkY+cutClZ+cuttZ+cutkZ
  81 */
  82 TMatrixD cutMatrix(4*7,2);
  83 const Double_t rmsCut09[6]={0.15,0.9,0.015, 0.2, 0.9, 0.015};
  84
  85
  86
  87 //
  88 //
  89 Int_t          toSkip       = 2;     //
  90 Int_t          toSkipOffset = 0;
  91 Int_t          toSkipTrack       = 2;
  92 Int_t          toSkipTrackOffset = 0;
  93 Int_t          isFilterTest = 0;
  94 //
  95 // Track cuts
  96 //
  97 TCut * cSide[4]={0,0,0,0};
  98 TCut *cP3[4]={0,0,0,0};
  99 TCut *cSP3[4]={0,0,0,0};
 100 TCut *cP4[4]={0,0,0,0};
 101 TCut *cM4[4]={0,0,0,0};
 102 TCut *cA[4]={0,0,0,0};
 103
 104 TCut cutAll="";
 105
 106
 107 //
 108 //
 109 //
 110 TChain *chainPoints=0;
 111 TEntryList *elist=0;
 112 AliTPCkalmanFit * kalmanFitNew=0;
 113 AliTPCkalmanFit * kalmanFitOrig=0;
 114 AliTPCkalmanFit * kalmanFitApply=0;
 115
 116 AliTPCkalmanFit *  CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump);
 117 void FilterTracks();
 118 AliTPCkalmanFit * SetupFit();
 119 //
 120 void              AddFitFieldCage(AliTPCkalmanFit *kalmanFit);
 121 void              AddPhiScaling(AliTPCkalmanFit *kalmanFit);
 122 void              AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin );
 123 void              AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
 124 void              AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit);
 125 void              AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit);
 126 void              AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
 127 void              AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin);
 128 void              AddDrift(AliTPCkalmanFit *kalmanFit);
 129
 130
 131 AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream);
 132
 133 TVectorD * EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1 , TTreeSRedirector *pcstream);
 134 AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset);
 135
 136 AliTPCkalmanFit *   FitPointsLinear(Int_t maxTracks, Int_t trackDump);
 137
 138 void CalibAlignKalman(Int_t npoints, Int_t maxFiles, Int_t startFile, Int_t trackDump, Int_t nSkipTrack, Int_t nSkipTrackOffset, Int_t nSkip, Int_t nSkipOffset, Int_t bfilterTest){
 139   ///
 140
 141   AliTPCTransformation::BuildBasicFormulas();
 142   toSkip=nSkip;
 143   toSkipOffset= nSkipOffset;
 144   toSkipTrack = nSkipTrack;
 145   toSkipTrackOffset = nSkipTrackOffset;
 146   isFilterTest = bfilterTest;
 147   //
 148   // read the transformation to be applied
 149   TFile ftrafo("kalmanFitApply.root");
 150   kalmanFitApply = (AliTPCkalmanFit *)ftrafo.Get("kalmanFitNew");
 151   if (kalmanFitApply) {
 152     printf("Loaded transforamtion\n");
 153     kalmanFitApply->DumpCalib("IROCOROC");
 154     kalmanFitApply->InitTransformation();
 155   }else{
 156     printf("Not trnasformation specified\n");
 157   }
 158   //
 159   //
 160   //
 161   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
 162   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
 163   //
 164   AliXRDPROOFtoolkit tool;
 165   chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0, maxFiles, startFile);
 166   CalibAlignKalmanFit(npoints, trackDump);
 167 }
 168
 169
 170
 171
 172 AliTPCkalmanFit *  CalibAlignKalmanFit(Int_t maxTracks, Int_t trackDump){
 173   /// Fitting procedure
 174
 175   AliTPCTransformation::BuildBasicFormulas();
 176   FilterTracks();
 177   kalmanFitNew     = SetupFit();
 178   kalmanFitOrig  = SetupFit();
 179   kalmanFitNew->Init();
 180   kalmanFitOrig->Init();
 181   return FitPointsLinear(maxTracks,trackDump);
 182 }
 183
 184
 185
 186 AliTPCkalmanFit * SetupFit(){
 187   ///
 188
 189   AliTPCkalmanFit *kalmanFit =  new AliTPCkalmanFit;
 190   AddFitFieldCage(kalmanFit);
 191   AddPhiScaling(kalmanFit);
 192   AddDrift(kalmanFit);
 193   AddZShift(kalmanFit,3,3);
 194   AddZTilting(kalmanFit,3,3);
 195   //  AddLocalXYMisalignment(kalmanFit);
 196   //  AddLocalXYMisalignmentSector(kalmanFit);
 197   AddAlignSectorFourier(kalmanFit,4,4);
 198   AddAlignOROCIROCFourier(kalmanFit,5,5);
 199   kalmanFit->Init();
 200   return kalmanFit;
 201 }
 202
 203
 204 void FilterTracks(){
 205   ///
 206
 207   cSide[0] = new TCut("cutAA","p0In.fP[1]>0&&p1In.fP[1]>0");
 208   cSide[1] = new TCut("cutCC","p0In.fP[1]<0&&p1In.fP[1]<0");
 209   cSide[2] = new TCut("cutAC","p0In.fP[1]>0&&p1In.fP[1]<0");
 210   cSide[3] = new TCut("cutCA","p0In.fP[1]<0&&p1In.fP[1]>0");
 211   //
 212   TH1F * phisP3 = new TH1F("hhisP3","hhisP3",100,-0.01,0.01);
 213   TH1F * phisSP3 = new TH1F("hhisSP3","hhisSP3",100,-0.001,0.001);
 214   TH1F * phisP4 = new TH1F("hhisP4","hhisP4",100,-0.1,0.1);
 215   TH1F * phisM4 = new TH1F("hhisM4","hhisM4",100,-0.01,0.01);
 216
 217   //
 218
 219   TF1 *fg = new TF1("fg","gaus");
 220   for (Int_t iter=0; iter<2; iter++){
 221     for (Int_t ic=0;ic<4;ic++){
 222       if (!cA[ic]){
 223         cA[ic]=new TCut;
 224         *cA[ic]= *cSide[ic];
 225       }
 226       //
 227       // cutP3
 228       //
 229       chainPoints->Draw("p0In.fP[3]+p1In.fP[3]>>hhisP3",*cA[ic],"goff");
 230       phisP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisP3->GetMaximumBin())-0.003,phisP3->GetBinCenter(phisP3->GetMaximumBin())+0.003);
 231       cP3[ic]=new TCut(Form("abs(p0In.fP[3]+p1In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
 232       cutMatrix(7*ic+0,0) = fg->GetParameter(1);
 233       cutMatrix(7*ic+0,1) = fg->GetParameter(2);
 234       //
 235       // cutSP3
 236       //
 237       chainPoints->Draw("p0Out.fP[3]-p0In.fP[3]>>hhisSP3",*cA[ic],"goff");
 238       phisSP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisSP3->GetMaximumBin())-0.0015,phisP3->GetBinCenter(phisSP3->GetMaximumBin())+0.0015);
 239       cSP3[ic]=new TCut(Form("abs(p0Out.fP[3]-p0In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
 240       cutMatrix(7*ic+1,0) = fg->GetParameter(1);
 241       cutMatrix(7*ic+1,1) = fg->GetParameter(2);
 242       //
 243       chainPoints->Draw("p1Out.fP[3]-p1In.fP[3]>>hhisSP3",*cA[ic],"goff");
 244       phisSP3->Fit(fg,"QNR","QNR",phisP3->GetBinCenter(phisSP3->GetMaximumBin())-0.0015,phisP3->GetBinCenter(phisSP3->GetMaximumBin())+0.0015);
 245       *cSP3[ic]+=Form("abs(p1Out.fP[3]-p1In.fP[3]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
 246       cutMatrix(7*ic+2,0) = fg->GetParameter(1);
 247       cutMatrix(7*ic+2,1) = fg->GetParameter(2);
 248       //
 249       // cutP4
 250       //
 251       chainPoints->Draw("p0Out.fP[4]>>hhisP4",*cA[ic],"goff");
 252       phisP4->Fit(fg,"QNR","QNR",phisP4->GetBinCenter(phisP4->GetMaximumBin())-0.03,phisP4->GetBinCenter(phisP4->GetMaximumBin())+0.03);
 253       cP4[ic]=new TCut(Form("abs(p0Out.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
 254       cutMatrix(7*ic+3,0) = fg->GetParameter(1);
 255       cutMatrix(7*ic+3,1) = fg->GetParameter(2);
 256       chainPoints->Draw("p1Out.fP[4]>>hhisP4",*cA[ic],"goff");
 257       phisP4->Fit(fg,"QNR","QNR",phisP4->GetBinCenter(phisP4->GetMaximumBin())-0.03,phisP4->GetBinCenter(phisP4->GetMaximumBin())+0.03);
 258       *cP4[ic]+=Form("abs(p1Out.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
 259       cutMatrix(7*ic+4,0) = fg->GetParameter(1);
 260       cutMatrix(7*ic+4,1) = fg->GetParameter(2);
 261
 262     //
 263       // cutM4
 264       //
 265       chainPoints->Draw("p0Out.fP[4]-p0In.fP[4]>>hhisM4",*cA[ic],"goff");
 266       phisM4->Fit(fg,"QNR","QNR",phisM4->GetBinCenter(phisM4->GetMaximumBin())-0.03,phisM4->GetBinCenter(phisM4->GetMaximumBin())+0.03);
 267       cM4[ic]=new TCut(Form("abs(p0Out.fP[4]-p0In.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma));
 268       cutMatrix(7*ic+5,0) = fg->GetParameter(1);
 269       cutMatrix(7*ic+5,1) = fg->GetParameter(2);
 270
 271       chainPoints->Draw("p1Out.fP[4]-p1In.fP[4]>>hhisM4",*cA[ic],"goff");
 272       phisM4->Fit(fg,"QNR","QNR",phisM4->GetBinCenter(phisM4->GetMaximumBin())-0.03,phisM4->GetBinCenter(phisM4->GetMaximumBin())+0.03);
 273       *cM4[ic]+=Form("abs(p1Out.fP[4]-p1In.fP[4]-%f)<%f",fg->GetParameter(1), fg->GetParameter(2)*kNsigma);
 274       cutMatrix(7*ic+6,0) = fg->GetParameter(1);
 275       cutMatrix(7*ic+6,1) = fg->GetParameter(2);
 276       //
 277       //
 278       //
 279       cA[ic]=new TCut;
 280       *cA[ic]= *cSide[ic]+*cP3[ic]+*cSP3[ic]+*cP4[ic]+*cM4[ic];
 281     }
 282   }
 283   cutMatrix.Print();
 284
 285   cutAll = (*cA[0])||(*cA[1])||(*cA[2])||(*cA[3])+"abs(mag)<0.01&&ncont>0&&p.fNPoints>120";
 286
 287   delete  phisP3; // = new TH1F("hhisP3","hhisP3",100,-0.01,0.01);
 288   delete  phisSP3; // = new TH1F("hhisSP3","hhisSP3",100,-0.001,0.001);
 289   delete  phisP4;// = new TH1F("hhisP4","hhisP4",100,-0.1,0.1);
 290   delete  phisM4;// = new TH1F("hhisM4","hhisM4",100,-0.01,0.01);
 291
 292
 293
 294   chainPoints->Draw(">>listEL",cutAll,"entryList");
 295   elist = (TEntryList*)gDirectory->Get("listEL");
 296   chainPoints->SetEntryList(elist);
 297   elist->SetDirectory(0);
 298 }
 299
 300
 301
 302 AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks, Int_t trackDump){
 303   // create debug streamers
 304
 305   TTreeSRedirector *pcstream      = new TTreeSRedirector("kalmanfitTPC.root");
 306   TTreeSRedirector *pcstreamOrig = new TTreeSRedirector("kalmanfitTPCOrig.root");
 307   pcstream->GetFile()->cd();
 308   cutMatrix.Write("cutMarix");
 309   elist->Write("eventList");
 310   //
 311   //
 312   AliTrackPointArray *pointsNS=0;
 313   Float_t mag=0;
 314   Int_t   time=0;
 315   AliExternalTrackParam *param0=0;
 316   AliExternalTrackParam *param1=0;
 317   chainPoints->SetBranchAddress("p.",&pointsNS);
 318   chainPoints->SetBranchAddress("p0In.",&param0);
 319   chainPoints->SetBranchAddress("p1In.",&param1);
 320   chainPoints->SetBranchAddress("mag",&mag);
 321   chainPoints->SetBranchAddress("time",&time);
 322   Int_t accepted=0;
 323   printf("\n*\n*\n*Selected entries = %d\n*\n*\n*",Int_t(elist->GetN()));
 324
 325   //
 326   for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
 327     if (itrack%toSkipTrack!=toSkipTrackOffset) continue;
 328     Int_t entry=chainPoints->GetEntryNumber(itrack);
 329     chainPoints->GetEntry(entry);
 330     if (accepted>maxTracks) break;
 331     //
 332     AliTrackPointArray *points = AliTPCkalmanFit::SortPoints(*pointsNS);
 333     if (kalmanFitApply)  kalmanFitApply->ApplyCalibration(points,-1.);
 334     //
 335     // estimate and filter scattering
 336     //
 337     TVectorD *vecRMS09 = EstimateScatteringKalmanLinear(*points,*param0,*param1,pcstream);
 338     if (!vecRMS09) continue;
 339     Bool_t isOK=kTRUE;
 340     if ((*vecRMS09)[0] >rmsCut09[0]) isOK=kFALSE;
 341     if ((*vecRMS09)[1]/(*vecRMS09)[0] >rmsCut09[1]) isOK=kFALSE;;
 342     if ((*vecRMS09)[2] >rmsCut09[2]) isOK=kFALSE;
 343     if ((*vecRMS09)[3] >rmsCut09[3]) isOK=kFALSE;
 344     if ((*vecRMS09)[4]/(*vecRMS09)[0] >rmsCut09[4]) isOK=kFALSE;
 345     if ((*vecRMS09)[5] >rmsCut09[5]) isOK=kFALSE;
 346     if (!isOK || isFilterTest) {
 347       delete points;
 348       continue;
 349     }
 350     kalmanFitNew->PropagateTime(time);
 351     //
 352     //
 353     for (Int_t idir=-1; idir<=1; idir++){
 354       AliTrackPointArray *pointsF = FilterPoints(*points,idir, pcstream);
 355       if (!pointsF) continue;
 356       AliTrackPointArray *spointsF = 0;
 357       // we skip points for alignemnt  but not for QA
 358       if (idir==0)  spointsF = SkipPoints(*pointsF, toSkip*2, toSkipOffset);
 359       if (idir!=0)  spointsF = SkipPoints(*pointsF, toSkip,   toSkipOffset);
 360       if (idir==0)  accepted++;
 361       //
 362       if (accepted%50==0) {
 363         kalmanFitNew->FitTrackLinear(*pointsF, pcstream);
 364       }else{
 365         if (idir==0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
 366         if (idir!=0) kalmanFitNew->FitTrackLinear(*spointsF, 0);
 367       }
 368       if (idir==0) kalmanFitNew->DumpTrackLinear(*pointsF,pcstream);
 369       if (idir==0) kalmanFitOrig->DumpTrackLinear(*pointsF,pcstreamOrig);
 370       if (accepted%trackDump==0) {
 371         printf("%d\n", accepted);
 372       }
 373       AliSysInfo::AddStamp("trackFit", accepted,itrack);
 374       delete pointsF;
 375       delete spointsF;
 376     }
 377     delete points;
 378   }
 379   pcstream->GetFile()->cd();
 380   kalmanFitNew->Write("kalmanFit");
 381   pcstreamOrig->GetFile()->cd();
 382   kalmanFitOrig->Write("kalmanFitOrig");
 383   pcstreamOrig->GetFile()->cd();
 384   if (kalmanFitApply) kalmanFitApply->Write("kalmanFitApply");
 385
 386   delete pcstream;
 387   delete pcstreamOrig;
 388   return kalmanFitNew;
 389 }
 390
 391 void  QAPointsLinear(Int_t maxTracks, Int_t trackDump){
 392   /// check  the consistency of kalman fit
 393   /// Apply transformation
 394
 395   // create debug streeamers
 396   TTreeSRedirector *pcstreamNonCalib      = new TTreeSRedirector("kalmanfitTPCQANonCalib.root");
 397   TTreeSRedirector *pcstreamCalib         = new TTreeSRedirector("kalmanfitTPCQACalib.root");
 398   TTreeSRedirector *pcstream=0;
 399   AliTPCkalmanFit  *kalmanFitters[6]={0,0,0,0,0,0};
 400   for (Int_t i=0;i<6;i++){
 401     kalmanFitters[i]=SetupFit();
 402   }
 403   //
 404   AliTrackPointArray *points=0;
 405   AliExternalTrackParam *param0=0;
 406   AliExternalTrackParam *param1=0;
 407   Float_t mag=0;
 408   Int_t   time=0;
 409   chainPoints->SetBranchAddress("p.",&points);
 410   chainPoints->SetBranchAddress("mag",&mag);
 411   chainPoints->SetBranchAddress("time",&time);
 412   chainPoints->SetBranchAddress("p0In.",&param0);
 413   chainPoints->SetBranchAddress("p1In.",&param1);
 414
 415   Int_t accepted=0;
 416   //
 417
 418
 419   for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
 420     if (itrack%toSkipTrack!=toSkipTrackOffset) continue;
 421     Int_t entry=chainPoints->GetEntryNumber(itrack);
 422     chainPoints->GetEntry(entry);
 423     if (accepted>maxTracks) break;
 424     //
 425
 426     AliTrackPointArray pointsCalib(*points);
 427     for (Int_t iscalib=0; iscalib<1;iscalib++){
 428       if (iscalib>0)  kalmanFitNew->ApplyCalibration(&pointsCalib,-1.);
 429       if (iscalib==0) pcstream=pcstreamNonCalib;
 430       if (iscalib>0)  pcstream=pcstreamCalib;
 431       for (Int_t idir=-1; idir<=1; idir++){
 432         AliTrackPointArray *pointsF = FilterPoints(pointsCalib,idir, pcstream);
 433         if (!pointsF) continue;
 434         //
 435         if (idir==0) accepted++;
 436         kalmanFitters[iscalib*3+idir+1]->DumpTrackLinear(*pointsF,0);
 437         EstimateScatteringKalmanLinear(*pointsF,*param0,*param1,pcstream);
 438         delete pointsF;
 439       }
 440     }
 441     if (accepted%trackDump==0) {
 442       printf("%d\n", accepted);
 443     }
 444   }
 445   pcstreamCalib->GetFile()->cd();
 446   kalmanFitters[0]->Write("fitUpNonCalib");
 447   kalmanFitters[1]->Write("fitUpDownNonCalib");
 448   kalmanFitters[2]->Write("fitDownNonCalib");
 449   kalmanFitters[3]->Write("fitUpCalib");
 450   kalmanFitters[4]->Write("fitUpDownCalib");
 451   kalmanFitters[5]->Write("fitDownCalib");
 452   delete pcstreamCalib;
 453   delete pcstreamNonCalib;
 454 }
 455
 456
 457 void  TestScattering(Int_t maxTracks, Int_t trackDump){
 458   /// test Multiple scattering algorithm
 459   /// Apply transformation
 460   ///
 461   /// create debug streeamers
 462
 463   TTreeSRedirector *pcstream      = new TTreeSRedirector("kalmanfitTPCMS.root");
 464   //
 465   //
 466   AliTrackPointArray *points=0;
 467   AliExternalTrackParam *param0=0;
 468   AliExternalTrackParam *param1=0;
 469   Float_t mag=0;
 470   Int_t   time=0;
 471   chainPoints->SetBranchAddress("p.",&points);
 472   chainPoints->SetBranchAddress("mag",&mag);
 473   chainPoints->SetBranchAddress("time",&time);
 474   chainPoints->SetBranchAddress("p0In.",&param0);
 475   chainPoints->SetBranchAddress("p1In.",&param1);
 476   Int_t accepted=0;
 477   //
 478   for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
 479     Int_t entry=chainPoints->GetEntryNumber(itrack);
 480     chainPoints->GetEntry(entry);
 481     if (accepted>maxTracks) break;
 482     //
 483     AliTrackPointArray *pointsSorted = AliTPCkalmanFit::SortPoints(*points);
 484     EstimateScatteringKalmanLinear(*pointsSorted,*param0,*param1,pcstream);
 485     accepted++;
 486     if (accepted%trackDump==0) {
 487       printf("%d\n", accepted);
 488     }
 489     delete pointsSorted;
 490   }
 491   delete pcstream;
 492 }
 493
 494
 495
 496
 497
 498 AliTrackPointArray *SkipPoints(AliTrackPointArray &points, Int_t nskip, Int_t nskipOffset){
 499   /// create new array with skipped points
 500
 501   Int_t npoints = points.GetNPoints();
 502   Int_t npointsF = (npoints-nskipOffset-1)/nskip;
 503   AliTrackPoint point;
 504   AliTrackPointArray *pointsF= new AliTrackPointArray(npointsF);
 505   Int_t used=0;
 506   for (Int_t ipoint=nskipOffset; ipoint<npoints; ipoint+=nskip){
 507     //
 508     if (!points.GetPoint(point,ipoint)) continue;
 509     pointsF->AddPoint(used,&point);
 510     used++;
 511     if (used==npointsF) break;
 512   }
 513   return pointsF;
 514 }
 515
 516
 517
 518 AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream){
 519   ///  Filter points - input points for KalmanFilter
 520
 521   TLinearFitter lfitY(2,"pol1");
 522   TLinearFitter lfitZ(2,"pol1");
 523   TVectorD vecZ(2);
 524   TVectorD vecY(2);
 525   //
 526   lfitY.StoreData(kTRUE);
 527   lfitZ.StoreData(kTRUE);
 528   Int_t npoints = points.GetNPoints();
 529   if (npoints<2) return 0;
 530   Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
 531   Double_t ca = TMath::Cos(currentAlpha);
 532   Double_t sa = TMath::Sin(currentAlpha);
 533   //
 534   // 1.b Fit the track in the rotated frame - MakeSeed
 535   //
 536   Double_t maxX =-10000, minX=10000;
 537   for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 538     Double_t rx =   ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
 539     Double_t ry =  -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
 540     Double_t rz =  points.GetZ()[ipoint];
 541     if (dir== 1 && rx<0) continue;
 542     if (dir==-1 && rx>0) continue;
 543     if (maxX<rx) maxX=rx;
 544     if (minX>rx) minX=rx;
 545     lfitY.AddPoint(&rx,ry,1);
 546     lfitZ.AddPoint(&rx,rz,1);
 547   }
 548   if (TMath::Abs(maxX-minX)<kArmCut) return 0;
 549   if (lfitY.GetNpoints()<knclCut) return 0;
 550   //
 551   lfitY.Eval();
 552   lfitZ.Eval();
 553   lfitY.GetParameters(vecY);
 554   lfitZ.GetParameters(vecZ);
 555   //
 556   Double_t chi2Y = lfitY.GetChisquare()/lfitY.GetNpoints();
 557   Double_t chi2Z = lfitZ.GetChisquare()/lfitZ.GetNpoints();
 558   if (TMath::Sqrt(chi2Y)>krmsYcutGlobal) return 0;
 559   if (TMath::Sqrt(chi2Y)>krmsYcutGlobal) return 0;
 560   //
 561   //
 562   Int_t accepted=0;
 563   Int_t toBeUsed    =0;
 564   AliTrackPoint point;
 565   AliTrackPointArray *pointsF=0;
 566   for (Int_t iter=0; iter<2;iter++){
 567     for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 568       //
 569       if (!points.GetPoint(point,ipoint)) continue;
 570       Double_t rx =   ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
 571       Double_t ry =  -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
 572       Double_t rz =  points.GetZ()[ipoint];
 573       if (dir== 1 && rx<0) continue;
 574       if (dir==-1 && rx>0) continue;
 575       Double_t erry = TMath::Sqrt(chi2Y);
 576       Double_t errz = TMath::Sqrt(chi2Z);
 577       Double_t fy = vecY[0]+vecY[1]*rx;
 578       Double_t fz = vecZ[0]+vecZ[1]*rx;
 579       if (TMath::Abs(fy-ry)>erry*kSigmaCut) continue;
 580       if (TMath::Abs(fz-rz)>errz*kSigmaCut) continue;
 581       accepted++;
 582       if (pointsF) pointsF->AddPoint(toBeUsed,&point);
 583       toBeUsed++;
 584     }
 585     if (pcstream){
 586       (*pcstream)<<"filter"<<
 587         "iter="<<iter<<
 588         "accepted="<<accepted<<
 589         "minX="<<minX<<
 590         "maxX="<<maxX<<
 591         "vY.="<<&vecY<<
 592         "vZ.="<<&vecZ<<
 593         "chi2Y="<<chi2Y<<
 594         "chi2Z="<<chi2Z<<
 595         "\n";
 596     }
 597     if (accepted<knclCut) break;
 598     if (iter==0) pointsF = new AliTrackPointArray(toBeUsed);
 599     accepted=0;
 600     toBeUsed=0;
 601   }
 602   return pointsF;
 603 }
 604
 605
 606 AliTrackPointArray * SortPoints(AliTrackPointArray &points){
 607   /// Creates the array  - points sorted according radius - neccessay for kalman fit
 608   ///
 609   /// 0. choose the frame - rotation angle
 610
 611   Int_t npoints = points.GetNPoints();
 612   if (npoints<1) return 0;
 613   Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
 614   Double_t ca = TMath::Cos(currentAlpha);
 615   Double_t sa = TMath::Sin(currentAlpha);
 616   //
 617   // 1. sort the points
 618   //
 619   Double_t *rxvector = new Double_t[npoints];
 620   Int_t    *indexes  = new Int_t[npoints];
 621   for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 622     rxvector[ipoint]=ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
 623   }
 624   TMath::Sort(npoints, rxvector,indexes,kFALSE);
 625   AliTrackPoint point;
 626   AliTrackPointArray *pointsSorted= new AliTrackPointArray(npoints);
 627   for (Int_t ipoint=0; ipoint<npoints; ipoint++){
 628     if (!points.GetPoint(point,indexes[ipoint])) continue;
 629     pointsSorted->AddPoint(ipoint,&point);
 630   }
 631   delete [] rxvector;
 632   delete [] indexes;
 633   return pointsSorted;
 634 }
 635
 636 TVectorD *  EstimateScatteringKalmanLinear(AliTrackPointArray &points, AliExternalTrackParam &p0, AliExternalTrackParam &p1 , TTreeSRedirector *pcstream){
 637   /// Algorithm - 0. Fit the track forward and backward
 638   ///           - 1. Store the current parameters in each point
 639
 640   const Int_t kMinPoints= 70;
 641   const Double_t kResY  = 0.1;
 642   const Double_t kResZ  = 0.1;
 643   const Double_t kMisY  = 0.02;
 644   const Double_t kMisZ  = 0.02;
 645   const Double_t kLArm  = 120.;
 646   const Double_t kACsideFac = 10.;
 647   const Double_t kMS0   = kResY/20.;
 648
 649   Int_t npoints = points.GetNPoints();
 650   if (npoints<kMinPoints) return 0;
 651
 652   TVectorD *vecPos[11]={0,0,0,0,0,0,0,0,0,0,0};
 653   for (Int_t i=0;i<11;i++){
 654     vecPos[i]=new TVectorD(npoints);
 655   }
 656
 657   //
 658   Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
 659   Double_t ca  = TMath::Cos(currentAlpha);
 660   Double_t sa  = TMath::Sin(currentAlpha);
 661   //
 662   //
 663   TMatrixD trParamY(2,1),trCovarY(2,2);
 664   TMatrixD trParamZ(2,1),trCovarZ(2,2);
 665   Int_t kNmeas  = 1;
 666   Int_t nelem   = 2;
 667   //
 668   TMatrixD matHk(kNmeas,nelem);     // vector to mesurement
 669   TMatrixD matHkT(nelem,kNmeas);    // helper matrix Hk transpose
 670   TMatrixD matFk(nelem,nelem);
 671   TMatrixD matFkT(nelem,nelem);
 672   TMatrixD vecYk(kNmeas,1);         // Innovation or measurement residual
 673   TMatrixD vecZk(kNmeas,1);         // Innovation or measurement residual
 674   TMatrixD measR(kNmeas,kNmeas);
 675   TMatrixD matSk(kNmeas,kNmeas);    // Innovation (or residual) covariance
 676   TMatrixD matKk(nelem,kNmeas);     // Optimal Kalman gain
 677   TMatrixD covXk2(nelem,nelem);     // helper matrix
 678   TMatrixD covXk3(nelem,nelem);     // helper matrix
 679   TMatrixD mat1(nelem,nelem);
 680   mat1(0,0)=1.; mat1(0,1)=0.;
 681   mat1(1,0)=0.; mat1(1,1)=1.;
 682   matHk(0,0)=1;
 683   matHk(0,1)=0;
 684
 685   Double_t lastX    = 0;
 686   Int_t lastVolId=-1;
 687   for (Int_t idir=0; idir<2;idir++){
 688     // fit direction
 689     //
 690     for (Int_t ip=0; ip<npoints; ip++){
 691       Int_t ipoint= ip;
 692       if (idir>0) ipoint = npoints-ip-1;
 693       Double_t rx =   ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
 694       Double_t ry =  -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
 695       Double_t rz =  points.GetZ()[ipoint];
 696       Int_t  volId= points.GetVolumeID()[ipoint];
 697       //
 698       if (ip==0){
 699         // set initital parameters and covariance - use first and middle point
 700         Double_t rxm =   ca*points.GetX()[npoints/2]+sa*points.GetY()[npoints/2];
 701         Double_t rym =  -sa*points.GetX()[npoints/2]+ca*points.GetY()[npoints/2];
 702         Double_t rzm =  points.GetZ()[npoints/2];
 703         trParamY(0,0) = ry;
 704         trParamY(1,0) = (rym-ry)/(rxm-rx);
 705         trParamZ(0,0) = rz;
 706         trParamZ(1,0) = (rzm-rz)/(rxm-rx);
 707         //
 708         trCovarY(0,0) = kResY*kResY;
 709         trCovarY(1,1) = (kResY*kResY)/((rxm-rx)*(rxm-rx));
 710         trCovarZ(0,0) = kResZ*kResZ;
 711         trCovarZ(1,1) = (kResZ*kResZ)/((rxm-rx)*(rxm-rx));
 712         lastX     = rx;
 713         lastVolId = volId;
 714       }
 715       //
 716       // Propagate
 717       //
 718       if ((volId%36)<18 && (lastVolId%36)>=18){
 719         // A - C side cross
 720         trCovarY(0,0)+=kMisY*kMisY*kACsideFac;
 721         trCovarZ(0,0)+=kMisZ*kMisZ*kACsideFac;
 722         trCovarY(1,1)+=kACsideFac*(kMisY*kMisY)/(kLArm*kLArm);;
 723         trCovarZ(1,1)+=kACsideFac*(kMisZ*kMisZ)/(kLArm*kLArm);
 724         lastVolId=volId;
 725       }
 726       if (volId!=lastVolId){
 727         // volumeID change
 728         trCovarY(0,0)+=kMisY*kMisY;
 729         trCovarZ(0,0)+=kMisZ*kMisZ;
 730         trCovarY(1,1)+=(kMisY*kMisY)/(kLArm*kLArm);;
 731         trCovarZ(1,1)+=(kMisZ*kMisZ)/(kLArm*kLArm);
 732         lastVolId=volId;
 733       }
 734       //
 735       // Propagate
 736       //
 737       Double_t deltaX=rx-lastX;
 738       trParamY(0,0)+=deltaX*trParamY(1,0);
 739       trParamZ(0,0)+=deltaX*trParamZ(1,0);
 740       matFk(0,0)=1; matFk(0,1)=deltaX;
 741       matFk(1,0)=0; matFk(1,1)=1;
 742       matFkT=matFk.T(); matFk.T();
 743       covXk2=matFk*trCovarY*matFkT;
 744       trCovarY=covXk2;
 745       covXk2=matFk*trCovarZ*matFkT;
 746       trCovarZ=covXk2;
 747
 748       // multiple scattering
 749       trCovarY(1,1)+=TMath::Abs(deltaX)*kMS0*kMS0;
 750       trCovarZ(1,1)+=TMath::Abs(deltaX)*kMS0*kMS0;
 751       lastX=rx;
 752       //
 753       // Update
 754       //
 755       for (Int_t coord=0; coord<2;coord++){
 756         TMatrixD* pvecXk = (coord==0)? &trParamY: &trParamZ;
 757         TMatrixD* pcovXk = (coord==0)? &trCovarY: &trCovarZ;
 758         //
 759         TMatrixD& vecXk = *pvecXk;
 760         TMatrixD& covXk = *pcovXk;
 761         measR(0,0) = (coord==0) ? kResY:kResZ;
 762         vecZk(0,0) = (coord==0) ? ry:rz;
 763         //
 764         vecYk = vecZk-matHk*vecXk;               // Innovation or measurement residual
 765         matHkT=matHk.T(); matHk.T();
 766         matSk = (matHk*(covXk*matHkT))+measR;    // Innovation (or residual) covariance
 767         matSk.Invert();
 768         matKk = (covXk*matHkT)*matSk;            //  Optimal Kalman gain
 769         //
 770         covXk2= (mat1-(matKk*matHk));
 771         covXk3 =  covXk2*covXk;
 772         vecXk += matKk*vecYk;                    //  updated vector
 773         covXk = covXk3;
 774       }
 775       //
 776       // store parameters
 777       //
 778       (*vecPos[0])[ipoint]=rx;
 779       (*vecPos[1])[ipoint]=ry;
 780       (*vecPos[2])[ipoint]=rz;
 781
 782       (*vecPos[4*idir+0+3])[ipoint]=trParamY(0,0);
 783       (*vecPos[4*idir+1+3])[ipoint]=trParamY(1,0);
 784       //
 785       (*vecPos[4*idir+2+3])[ipoint]=trParamZ(0,0);
 786       (*vecPos[4*idir+3+3])[ipoint]=trParamZ(1,0);
 787     }
 788   }
 789   //
 790   //
 791   //
 792   TVectorD vec(npoints);
 793   TVectorD rms(6);
 794   TVectorD rms09(6); // robust RMS - fraction 0.9
 795   TVectorD mean09(6); // robust RMS - fraction 0.9
 796   Double_t meanR,rmsR;
 797   Int_t npoints09 = Int_t(npoints*0.9);
 798   //
 799   vec=(*(vecPos[3])-*(vecPos[1]));
 800   rms[0]=TMath::RMS(npoints, vec.GetMatrixArray());   // rms cluster y
 801   AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
 802   rms09[0]=rmsR;
 803   mean09[0]=meanR;
 804   vec=(*(vecPos[7])-*(vecPos[3]));
 805   rms[1]=TMath::RMS(npoints, vec.GetMatrixArray());   // rms track y
 806   AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
 807   rms09[1]=rmsR;
 808   mean09[1]=meanR;
 809   vec=(*(vecPos[8])-*(vecPos[4]));
 810   rms[2]=TMath::RMS(npoints, vec.GetMatrixArray());   // rms track ky
 811   AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
 812   rms09[2]=rmsR;
 813   mean09[2]=meanR;
 814   //
 815   vec=(*(vecPos[5])-*(vecPos[2]));
 816   rms[3]=TMath::RMS(npoints, vec.GetMatrixArray());   // rms cluster z
 817   AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
 818   rms09[3]=rmsR;
 819   mean09[3]=meanR;
 820   vec=(*(vecPos[9])-*(vecPos[5]));
 821   rms[4]=TMath::RMS(npoints, vec.GetMatrixArray());   // rms track z
 822   AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
 823   rms09[4]=rmsR;
 824   mean09[4]=meanR;
 825   vec=(*(vecPos[10])-*(vecPos[6]));
 826   rms[5]=TMath::RMS(npoints, vec.GetMatrixArray());   // rms track kz
 827   AliMathBase::EvaluateUni(npoints, vec.GetMatrixArray(),meanR,rmsR,npoints09);
 828   rms09[5]=rmsR;
 829   mean09[5]=meanR;
 830
 831
 832   (*pcstream)<<"kf"<<
 833     "p.="<<&points<<
 834     "p0.="<<&p0<<
 835     "p1.="<<&p1<<
 836     "rms.="<<&rms<<
 837     "rms09.="<<&rms09<<
 838     "mean09.="<<&mean09<<
 839     "py.="<<&trParamY<<
 840     "pz.="<<&trParamZ<<
 841     "cy.="<<&trCovarY<<
 842     "cz.="<<&trCovarY<<
 843     "\n";
 844   for (Int_t i=0;i<11;i++){
 845     delete vecPos[i];
 846   }
 847   return new TVectorD(rms09);
 848 }
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858 void  AddFitFieldCage(AliTPCkalmanFit *kalmanFit){
 859   /// Add radial scaling due field cage
 860
 861   TVectorD fpar(10);
 862   AliTPCTransformation * transformation=0;
 863   char tname[100];
 864   //
 865   // linear R scaling and shift
 866   //
 867   for (Int_t iside=0; iside<=1; iside++)
 868     for (Int_t ipolR=0; ipolR<2; ipolR++){
 869       for (Int_t ipolZ=0; ipolZ<3; ipolZ++){
 870         fpar[0]=ipolR;
 871         fpar[1]=ipolZ;
 872         if (ipolR+ipolZ==0) continue;
 873         sprintf(tname,"tTPCscalingRPolR%dDr%dSide%d",ipolR,ipolZ,iside);
 874         transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,  1);
 875         transformation->SetParams(0,0.2,0,&fpar);
 876         kalmanFit->AddCalibration(transformation);
 877         //
 878       }
 879     }
 880   //
 881   //
 882   //Inner field cage
 883   for (Int_t iside=0; iside<=1; iside++)
 884     for (Int_t ipol=0; ipol<3; ipol++){
 885       fpar[0]=ipol;
 886       sprintf(tname,"tTPCscalingRIFC%dSide%d",ipol,iside);
 887       transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRIFC",0,0,   1);
 888       transformation->SetParams(0,0.2,0,&fpar);
 889       kalmanFit->AddCalibration(transformation);
 890     }
 891   //
 892   //
 893   //Outer field cage
 894   for (Int_t iside=0; iside<=1; iside++)
 895     for (Int_t ipol=0; ipol<3; ipol++){
 896       fpar[0]=ipol;
 897       //Outer field cage
 898       sprintf(tname,"tTPCscalingROFC%dSide%d",ipol,iside);
 899       transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingROFC",0,0,  1);
 900       transformation->SetParams(0,0.2,0,&fpar);
 901       kalmanFit->AddCalibration(transformation);
 902     }
 903 }
 904
 905
 906 void AddPhiScaling(AliTPCkalmanFit *kalmanFit){
 907   /// Add linear local phi scaling -
 908   /// separate IROC/OROC  - A side/C side
 909   ///    "tscalingLocalPhiIROC"
 910   ///    "tscalingLocalPhiOROC"
 911
 912   TBits maskInner(72);
 913   TBits maskOuter(72);
 914   for (Int_t i=0;i<72;i++){
 915     if (i<36){
 916       maskInner[i]=kTRUE;
 917     }
 918     if (i>=36){
 919       maskOuter[i]=kTRUE;
 920     }
 921   }
 922   TVectorD fpar(10);
 923   AliTPCTransformation * transformation=0;
 924   fpar[0]=1;
 925   transformation = new AliTPCTransformation("tscalingLocalPhiIROC", new TBits(maskInner), 0,"TPCscalingPhiLocal",0,  1);
 926   transformation->SetParams(0,0.02,0,&fpar);
 927   kalmanFit->AddCalibration(transformation);
 928   transformation = new AliTPCTransformation("tscalingLocalPhiOROC", new TBits(maskOuter), 0,"TPCscalingPhiLocal",0,  1);
 929   transformation->SetParams(0,0.02,0,&fpar);
 930   kalmanFit->AddCalibration(transformation);
 931   //
 932 }
 933
 934 void AddDrift(AliTPCkalmanFit *kalmanFit){
 935   /// Add drift velocity transformation
 936
 937   TVectorD fpar(10);
 938   AliTPCTransformation * transformation=0;
 939   fpar[0]=1;
 940   transformation = new AliTPCTransformation("tTPCscalingZDrift", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDr",  0);
 941   transformation->SetParams(0,4.,1.,&fpar);
 942   kalmanFit->AddCalibration(transformation);
 943   //
 944   transformation = new AliTPCTransformation("tTPCscalingZDriftGy", AliTPCTransformation::BitsAll(), 0,0,"TPCscalingZDrGy",  0);
 945   transformation->SetParams(0,0.2,0.0,&fpar);
 946   kalmanFit->AddCalibration(transformation);
 947 }
 948
 949
 950
 951
 952 void  AddZShift(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
 953   ///
 954
 955   TVectorD fpar(10);
 956   fpar[0]=0; fpar[1]=0; fpar[2]=0;
 957   char tname[1000];
 958   AliTPCTransformation * transformation=0;
 959   //
 960   //
 961   //
 962   for (Int_t i=0; i<=ncos;i++){
 963     fpar[0]=i;  // cos frequency
 964     fpar[1]=0;  // no sinus
 965     fpar[2]=1;  // relative misalignment
 966     //
 967     sprintf(tname,"tTPCDeltaZIROCOROCSideA_Cos%d",i);
 968     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" ,  0);
 969     transformation->SetParams(0,0.03,0,&fpar);
 970     kalmanFit->AddCalibration(transformation);
 971     //
 972     sprintf(tname,"tTPCDeltaZIROCOROCSideC_Cos%d",i);
 973     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" ,  0);
 974     transformation->SetParams(0,0.03,0,&fpar);
 975     kalmanFit->AddCalibration(transformation);
 976     //
 977     //
 978     //
 979     fpar[2]=0;  // absolute  misalignment
 980     sprintf(tname,"tTPCDeltaZSectorSideA_Cos%d",i);
 981     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" ,  0);
 982     transformation->SetParams(0,0.1,0,&fpar);
 983     kalmanFit->AddCalibration(transformation);
 984     //
 985     sprintf(tname,"tTPCDeltaZSectorSideC_Cos%d",i);
 986     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" ,  0);
 987     transformation->SetParams(0,0.1,0,&fpar);
 988     kalmanFit->AddCalibration(transformation);
 989   }
 990
 991   for (Int_t i=1; i<=nsin;i++){
 992     fpar[0]=0;  // cos frequency
 993     fpar[1]=i;  // sinus frequncy
 994     fpar[2]=1;  // relative misalignment
 995     //
 996     sprintf(tname,"tTPCDeltaZIROCOROCSideA_Sin%d",i);
 997     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" ,  0);
 998     transformation->SetParams(0,0.03,0,&fpar);
 999     kalmanFit->AddCalibration(transformation);
1000     //
1001     sprintf(tname,"tTPCDeltaZIROCOROCSideC_Sin%d",i);
1002     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" ,  0);
1003     transformation->SetParams(0,0.03,0,&fpar);
1004     kalmanFit->AddCalibration(transformation);
1005     //
1006     //
1007     //
1008     fpar[2]=0;  // absolute  misalignment
1009     sprintf(tname,"tTPCDeltaZSectorSideA_Sin%d",i);
1010     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCDeltaZ" ,  0);
1011     transformation->SetParams(0,0.1,0,&fpar);
1012     kalmanFit->AddCalibration(transformation);
1013     //
1014     sprintf(tname,"tTPCDeltaZSectorSideC_Sin%d",i);
1015     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCDeltaZ" ,  0);
1016     transformation->SetParams(0,0.1,0,&fpar);
1017     kalmanFit->AddCalibration(transformation);
1018   }
1019
1020 }
1021
1022
1023
1024
1025 void AddZTilting(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1026   /// z tilting absolute (sector) and relative (IROC-OROC)
1027
1028   TVectorD fpar(10);
1029   fpar[0]=0; fpar[1]=0; fpar[2]=0;
1030   char tname[1000];
1031   AliTPCTransformation * transformation=0;
1032   //
1033   //
1034   //
1035   for (Int_t i=0; i<=ncos;i++){
1036     fpar[0]=i;  // cos frequency
1037     fpar[1]=0;  // sinus frequency
1038     fpar[2]=1;  // relative misalignment
1039     //
1040     sprintf(tname,"tTPCTiltingZIROCOROCSideA_Cos%d",i);
1041     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" ,  0);
1042     transformation->SetParams(0,0.03,0,&fpar);
1043     kalmanFit->AddCalibration(transformation);
1044     //
1045     sprintf(tname,"tTPCTiltingZIROCOROCSideC_Cos%d",i);
1046     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" ,  0);
1047     transformation->SetParams(0,0.03,0,&fpar);
1048     kalmanFit->AddCalibration(transformation);
1049     //
1050     //
1051     //
1052     fpar[2]=0;  // absolute  misalignment
1053     sprintf(tname,"tTPCTiltingZSectorSideA_Cos%d",i);
1054     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" ,  0);
1055     transformation->SetParams(0,0.1,0,&fpar);
1056     kalmanFit->AddCalibration(transformation);
1057     //
1058     sprintf(tname,"tTPCTiltingZSectorSideC_Cos%d",i);
1059     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" ,  0);
1060     transformation->SetParams(0,0.1,0,&fpar);
1061     kalmanFit->AddCalibration(transformation);
1062   }
1063
1064   for (Int_t i=1; i<=nsin;i++){
1065     fpar[0]=0;  // cos frequency
1066     fpar[1]=i;  // sinus frequncy
1067     fpar[2]=1;  // relative misalignment
1068     //
1069     sprintf(tname,"tTPCTiltingZIROCOROCSideA_Sin%d",i);
1070     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" ,  0);
1071     transformation->SetParams(0,0.03,0,&fpar);
1072     kalmanFit->AddCalibration(transformation);
1073     //
1074     sprintf(tname,"tTPCTiltingZIROCOROCSideC_Sin%d",i);
1075     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" ,  0);
1076     transformation->SetParams(0,0.03,0,&fpar);
1077     kalmanFit->AddCalibration(transformation);
1078     //
1079     //
1080     //
1081     fpar[2]=0;  // absolute  misalignment
1082     sprintf(tname,"tTPCTiltingZSectorSideA_Sin%d",i);
1083     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), 0,0, "TPCTiltingZ" ,  0);
1084     transformation->SetParams(0,0.1,0,&fpar);
1085     kalmanFit->AddCalibration(transformation);
1086     //
1087     sprintf(tname,"tTPCTiltingZSectorSideC_Sin%d",i);
1088     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), 0,0, "TPCTiltingZ" ,  0);
1089     transformation->SetParams(0,0.1,0,&fpar);
1090     kalmanFit->AddCalibration(transformation);
1091   }
1092 }
1093
1094
1095
1096 void  AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit){
1097   ///
1098
1099   TVectorD fpar(10);
1100   AliTPCTransformation * transformation=0;
1101   TBits maskInnerA(72);
1102   TBits maskInnerC(72);
1103   for (Int_t i=0;i<72;i++){
1104     if (i<36){
1105       if (i%36<18)  maskInnerA[i]=kTRUE;
1106       if (i%36>=18) maskInnerC[i]=kTRUE;
1107     }
1108   }
1109   //
1110   //
1111   transformation = new AliTPCTransformation("tTPCDeltaLxIROCA", new TBits(maskInnerA), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
1112   transformation->SetParams(0,0.2,0,&fpar);
1113   kalmanFit->AddCalibration(transformation);
1114   transformation = new AliTPCTransformation("tTPCDeltaLxIROCC", new TBits(maskInnerC), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
1115   transformation->SetParams(0,0.2,0,&fpar);
1116   kalmanFit->AddCalibration(transformation);
1117   //
1118   transformation = new AliTPCTransformation("tTPCDeltaLyIROCA", new TBits(maskInnerA), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1119   transformation->SetParams(0,0.2,0,&fpar);
1120   kalmanFit->AddCalibration(transformation);
1121   transformation = new AliTPCTransformation("tTPCDeltaLyIROCC", new TBits(maskInnerC), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1122   transformation->SetParams(0,0.2,0,&fpar);
1123   kalmanFit->AddCalibration(transformation);
1124 }
1125
1126 void  AddLocalXYMisalignmentSector(AliTPCkalmanFit *kalmanFit){
1127   ///
1128
1129   TVectorD fpar(10);
1130   AliTPCTransformation * transformation=0;
1131   Int_t fixSector =4;
1132   //
1133   for (Int_t isec=0; isec<36;isec++){
1134     TBits mask(72);
1135     mask[isec]=kTRUE;
1136     char tname[1000];
1137     //
1138     sprintf(tname,"tTPClocalLxIROC%d",isec);
1139     transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1140     transformation->SetParams(0,0.2,0,&fpar);
1141     kalmanFit->AddCalibration(transformation);
1142     //
1143     sprintf(tname,"tTPClocalLyIROC%d",isec);
1144     transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1145     transformation->SetParams(0,0.2,0,&fpar);
1146     kalmanFit->AddCalibration(transformation);
1147
1148     sprintf(tname,"tTPClocalRzIROC%d",isec);
1149     transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldRzdGY",0,  0);
1150     transformation->SetParams(0,0.2,0,&fpar);
1151     kalmanFit->AddCalibration(transformation);
1152
1153   }
1154   //
1155   for (Int_t isec=0; isec<36;isec++){
1156     if (isec%18==fixSector) continue;
1157     TBits mask(72);
1158     mask[isec]   =kTRUE;
1159     mask[isec+36]=kTRUE;
1160     char tname[1000];
1161     //
1162     sprintf(tname,"tTPClocalLxSector%d",isec);
1163     transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1164     transformation->SetParams(0,0.2,0,&fpar);
1165     kalmanFit->AddCalibration(transformation);
1166     //
1167     sprintf(tname,"tTPClocalLySector%d",isec);
1168     transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1169     transformation->SetParams(0,0.2,0,&fpar);
1170     kalmanFit->AddCalibration(transformation);
1171
1172     sprintf(tname,"tTPClocalRzSector%d",isec);
1173     transformation = new AliTPCTransformation(tname, new TBits(mask), "TPClocaldLydGX","TPClocaldRzdGY",0,  0);
1174     transformation->SetParams(0,0.2,0,&fpar);
1175     kalmanFit->AddCalibration(transformation);
1176   }
1177
1178
1179 }
1180
1181
1182 void  AddAlignOROCIROCFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1183   ///
1184
1185   TVectorD fpar(10);
1186   AliTPCTransformation * transformation=0;
1187
1188   for (Int_t i=0; i<=ncos;i++){
1189     char tname[1000];
1190     fpar[0]=i;  // cos frequency
1191     fpar[1]=0;  // no sinus
1192     fpar[2]=1;  // relative misalignment
1193     //
1194     // Local X shift
1195     //
1196     sprintf(tname,"tTPClocalLxIROCOROCSideA_Cos%d",i);
1197     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
1198     transformation->SetParams(0,0.03,0,&fpar);
1199     kalmanFit->AddCalibration(transformation);
1200     sprintf(tname,"tTPClocalLxIROCOROCSideC_Cos%d",i);
1201     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
1202     transformation->SetParams(0,0.03,0,&fpar);
1203     kalmanFit->AddCalibration(transformation);
1204     //
1205     //
1206     // Local Y shift
1207     //
1208     sprintf(tname,"tTPClocalLyIROCOROCSideA_Cos%d",i);
1209     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1210     transformation->SetParams(0,0.03,0,&fpar);
1211     kalmanFit->AddCalibration(transformation);
1212     sprintf(tname,"tTPClocalLyIROCOROCSideC_Cos%d",i);
1213     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1214     transformation->SetParams(0,0.03,0,&fpar);
1215     kalmanFit->AddCalibration(transformation);
1216     //
1217     //
1218     //
1219     // Z rotation
1220     //
1221     sprintf(tname,"tTPClocalRzIROCOROCSideA_Cos%d",i);
1222     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0,  0);
1223     transformation->SetParams(0,0.3,0,&fpar);
1224     kalmanFit->AddCalibration(transformation);
1225     sprintf(tname,"tTPClocalRzIROCOROCSideC_Cos%d",i);
1226     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0,  0);
1227     transformation->SetParams(0,0.3,0,&fpar);
1228     kalmanFit->AddCalibration(transformation);
1229   }
1230   //
1231   for (Int_t i=1; i<=nsin;i++){
1232     char tname[1000];
1233     fpar[0]=0;  // cos frequency
1234     fpar[1]=i;  // sinus frequency
1235     fpar[2]=1;  // relative misalignment
1236     //
1237     // Local X shift
1238     //
1239     sprintf(tname,"tTPClocalLxIROCOROCSideA_Sin%d",i);
1240     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
1241     transformation->SetParams(0,0.03,0,&fpar);
1242     kalmanFit->AddCalibration(transformation);
1243     sprintf(tname,"tTPClocalLxIROCOROCSideC_Sin%d",i);
1244     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
1245     transformation->SetParams(0,0.03,0,&fpar);
1246     kalmanFit->AddCalibration(transformation);
1247     //
1248     //
1249     // Local Y shift
1250     //
1251     sprintf(tname,"tTPClocalLyIROCOROCSideA_Sin%d",i);
1252     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1253     transformation->SetParams(0,0.03,0,&fpar);
1254     kalmanFit->AddCalibration(transformation);
1255     sprintf(tname,"tTPClocalLyIROCOROCSideC_Sin%d",i);
1256     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1257     transformation->SetParams(0,0.03,0,&fpar);
1258     kalmanFit->AddCalibration(transformation);
1259     //
1260     //
1261     //
1262     // Z rotation
1263     //
1264     sprintf(tname,"tTPClocalRzIROCOROCSideA_Sin%d",i);
1265     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0,  0);
1266     transformation->SetParams(0,0.3,0,&fpar);
1267     kalmanFit->AddCalibration(transformation);
1268     sprintf(tname,"tTPClocalRzIROCOROCSideC_Sin%d",i);
1269     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0,  0);
1270     transformation->SetParams(0,0.3,0,&fpar);
1271     kalmanFit->AddCalibration(transformation);
1272   }
1273 }
1274
1275 void  AddAlignSectorFourier(AliTPCkalmanFit *kalmanFit, Int_t ncos, Int_t nsin){
1276   ///
1277
1278   TVectorD fpar(10);
1279   AliTPCTransformation * transformation=0;
1280
1281   for (Int_t i=0; i<=ncos;i++){
1282     char tname[1000];
1283     fpar[0]=i;  // cos frequency
1284     fpar[1]=0;  // no sinus
1285     fpar[2]=0;  // absolute misalignment
1286     if (i>0){
1287       //
1288       // A side is reference
1289       // local x
1290       sprintf(tname,"tTPClocalLxSectorSideA_Cos%d",i);
1291       transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
1292       transformation->SetParams(0,0.03,0,&fpar);
1293       kalmanFit->AddCalibration(transformation);
1294       if (i>1){
1295         //
1296         // Local Y shift
1297         //
1298         sprintf(tname,"tTPClocalLySectorSideA_Cos%d",i);
1299         transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1300         transformation->SetParams(0,0.03,0,&fpar);
1301         kalmanFit->AddCalibration(transformation);
1302       }
1303       //
1304       //
1305     }
1306     //
1307     // C side to vary
1308     // local x
1309     sprintf(tname,"tTPClocalLxSectorSideC_Cos%d",i);
1310     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
1311     transformation->SetParams(0,0.03,0,&fpar);
1312     kalmanFit->AddCalibration(transformation);
1313     //
1314     // Local Y shift
1315     //
1316     sprintf(tname,"tTPClocalLySectorSideC_Cos%d",i);
1317     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1318     transformation->SetParams(0,0.03,0,&fpar);
1319     kalmanFit->AddCalibration(transformation);
1320     //
1321     //
1322     // Z rotation - independent
1323     //
1324     sprintf(tname,"tTPClocalRzSectorSideA_Cos%d",i);
1325     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0,  0);
1326     transformation->SetParams(0,0.3,0,&fpar);
1327     kalmanFit->AddCalibration(transformation);
1328     sprintf(tname,"tTPClocalRzSectorSideC_Cos%d",i);
1329     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0,  0);
1330     transformation->SetParams(0,0.3,0,&fpar);
1331     kalmanFit->AddCalibration(transformation);
1332   }
1333
1334
1335
1336
1337   //
1338   //
1339   //
1340   for (Int_t i=1; i<=nsin;i++){
1341     char tname[1000];
1342     fpar[0]=0;  // non cos frequency
1343     fpar[1]=i;  // sinus frequency
1344     fpar[2]=0;  // absolute misalignment
1345     //
1346     // Local X shift
1347     //
1348     sprintf(tname,"tTPClocalLxSectorSideA_Sin%d",i);
1349     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
1350     transformation->SetParams(0,0.03,0,&fpar);
1351     kalmanFit->AddCalibration(transformation);
1352     sprintf(tname,"tTPClocalLxSectorSideC_Sin%d",i);
1353     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
1354     transformation->SetParams(0,0.03,0,&fpar);
1355     kalmanFit->AddCalibration(transformation);
1356     //
1357     //
1358     // Local Y shift
1359     //
1360     sprintf(tname,"tTPClocalLySectorSideA_Sin%d",i);
1361     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1362     transformation->SetParams(0,0.03,0,&fpar);
1363     kalmanFit->AddCalibration(transformation);
1364     sprintf(tname,"tTPClocalLySectorSideC_Sin%d",i);
1365     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
1366     transformation->SetParams(0,0.03,0,&fpar);
1367     kalmanFit->AddCalibration(transformation);
1368     //
1369     //
1370     //
1371     // Z rotation
1372     //
1373     sprintf(tname,"tTPClocalRzSectorSideA_Sin%d",i);
1374     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(0), "TPClocaldRzdGX","TPClocaldRzdGY",0,  0);
1375     transformation->SetParams(0,0.3,0,&fpar);
1376     kalmanFit->AddCalibration(transformation);
1377     sprintf(tname,"tTPClocalRzSectorSideC_Sin%d",i);
1378     transformation = new AliTPCTransformation(tname, AliTPCTransformation::BitsSide(1), "TPClocaldRzdGX","TPClocaldRzdGY",0,  0);
1379     transformation->SetParams(0,0.3,0,&fpar);
1380     kalmanFit->AddCalibration(transformation);
1381   }
1382 }
1383
1384 void SelectPixel(){
1385   for (Int_t i=0;i<8;i++){
1386     kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(-10,10);
1387     kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(-10,10);
1388     kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-15,15);
1389     kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-15,15);
1390   }
1391 }
1392
1393 void SelectNonPixelA(){
1394   for (Int_t i=0;i<8;i++){
1395     kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1396     kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1397     kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(0,250);
1398     kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(0,250);
1399   }
1400 }
1401
1402 void SelectNonPixelC(){
1403   for (Int_t i=0;i<8;i++){
1404     kalmanFitNew->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1405     kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(2)->SetRangeUser(10,250);
1406     kalmanFitNew->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-250,0);
1407     kalmanFitOrig->fLinearTrackDelta[i]->GetAxis(3)->SetRangeUser(-250,0);
1408   }
1409 }
1410
1411
1412 void DumpQA1D(  TObjArray &arrayOut){
1413   ///
1414
1415   TF1 fg("fg","gaus");
1416   TMatrixD sideARMS(8,2);
1417   TMatrixD sideCRMS(8,2);
1418   TMatrixD sideACRMS(8,2);
1419   TH1 *his=0;
1420   //
1421   // A side
1422   //
1423   SelectNonPixelA();
1424   for (Int_t i=0; i<8;i++){
1425     his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1426     his->Fit(&fg,"","", -0.15,0.15);
1427     sideARMS(i,0) = fg.GetParameter(2);
1428     his->SetDirectory(0);
1429     his->SetName(Form("Original SideA_%s",his->GetName()));
1430     his->SetTitle(Form("Original SideA_%s",his->GetTitle()));
1431     arrayOut.AddLast(his);
1432     his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1433     his->Fit(&fg,"","", -0.15,0.15);
1434     sideARMS(i,1) = fg.GetParameter(2);
1435     his->SetDirectory(0);
1436     his->SetName(Form("Aligned SideA_%s",his->GetName()));
1437     his->SetTitle(Form("Aligned SideA_%s",his->GetTitle()));
1438     arrayOut.AddLast(his);
1439   }
1440   //
1441   // C side
1442   //
1443   SelectNonPixelC();
1444   for (Int_t i=0; i<8;i++){
1445     his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1446     his->Fit(&fg,"","", -0.15,0.15);
1447     sideCRMS(i,0) = fg.GetParameter(2);
1448     his->SetDirectory(0);
1449     his->SetName(Form("Original SideC_%s",his->GetName()));
1450     his->SetTitle(Form("Original SideC_%s",his->GetTitle()));
1451     arrayOut.AddLast(his);
1452     his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1453     his->Fit(&fg,"","", -0.15,0.15);
1454     sideCRMS(i,1) = fg.GetParameter(2);
1455     his->SetDirectory(0);
1456     his->SetName(Form("Aligned SideC_%s",his->GetName()));
1457     his->SetTitle(Form("Aligned SideC_%s",his->GetTitle()));
1458     arrayOut.AddLast(his);
1459   }
1460   //
1461   // AC side
1462   //
1463   SelectPixel();
1464   for (Int_t i=0; i<8;i++){
1465     his = kalmanFitOrig->fLinearTrackDelta[i]->Projection(0);
1466     his->Fit(&fg,"","", -0.15,0.15);
1467     sideACRMS(i,0) = fg.GetParameter(2);
1468     his->SetDirectory(0);
1469     his->SetName(Form("Original SideAC_%s",his->GetName()));
1470     his->SetTitle(Form("Original SideAC_%s",his->GetTitle()));
1471     arrayOut.AddLast(his);
1472     his = kalmanFitNew->fLinearTrackDelta[i]->Projection(0);
1473     his->Fit(&fg,"","", -0.15,0.15);
1474     sideACRMS(i,1) = fg.GetParameter(2);
1475     his->SetDirectory(0);
1476     his->SetName(Form("Aligned SideC_%s",his->GetName()));
1477     his->SetTitle(Form("Aligned SideC_%s",his->GetTitle()));
1478     arrayOut.AddLast(his);
1479   }
1480   printf("DumQA\n");
1481   sideARMS.Print();
1482   sideCRMS.Print();
1483   sideACRMS.Print();
1484 }
1485
1486 void MakeFolders(TObjArray * arrayOut){
1487   ///
1488
1489   TFolder *folderBase = new TFolder("TPC align","TPC align");
1490   //
1491   //
1492   TString maskDelta[8];
1493   TString maskAlign[2]={"Orig","Alig");
1494   for (Int_t i=0; i<8;i++){
1495     maskDelta[i] = kalmanFitNew->fLinearTrackDelta[i]->GetName();
1496   }
1497
1498   Int_t entries = arrayOut->GetEntries();
1499   //
1500 }
1501
1502
1503
1504
1505 void MergeKalman(const char * list = "kalmanFit.list"){
1506   ///
1507
1508   ifstream in;
1509   in.open(list);
1510   TString currentFile;
1511   kalmanFitNew= 0;
1512   Int_t counter=0;
1513   while(in.good()) {
1514     //
1515     // calibrated
1516     //
1517     in >> currentFile;
1518     printf("%d\t%d\t%s\n", counter,currentFile.Length(),currentFile.Data());
1519     if (currentFile.Length()==0) continue;
1520     TFile * ffit = TFile::Open(currentFile.Data());
1521     TEntryList * tlist = (TEntryList*) ffit->Get("eventList");
1522     TMatrixD * cMatrix = (TMatrixD*) ffit->Get("cutMarix");
1523     if (tlist&&cMatrix){
1524       printf("Track entries=%d\n",tlist->GetN());
1525       if (cMatrix->Sum()<=0.00000001) {
1526         printf("Problem with track selection\n");
1527         continue;
1528       }
1529     }else{
1530       printf("Problem with track selection\n");
1531       continue;
1532     }
1533     //
1534
1535     AliTPCkalmanFit * fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFit");
1536     if (!fit) continue;
1537     if (!kalmanFitNew) {kalmanFitNew= fit; continue;};
1538     kalmanFitNew->Add(fit);
1539     printf("Selected entries=%f\n",fit->fLinearTrackDelta[0]->GetEntries());
1540     //delete tlist;
1541     //delete cMatrix;
1542     delete fit;
1543     delete ffit;
1544     //
1545     // original
1546     //
1547     currentFile.ReplaceAll("TPC","TPCOrig");
1548     printf("%d\t%d\t%s\n", counter,currentFile.Length(),currentFile.Data());
1549     if (currentFile.Length()==0) continue;
1550     ffit = TFile::Open(currentFile.Data());
1551     fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFitOrig");
1552     if (!fit) continue;
1553     if (!kalmanFitOrig) {kalmanFitOrig= fit; continue;};
1554     kalmanFitOrig->Add(fit);
1555     delete fit;
1556     delete ffit;
1557     counter++;
1558   }
1559   //
1560   // save merged results
1561   //
1562   TFile f("mergeKalmanFit.root","recreate");
1563   kalmanFitNew->Write("kalmanFitNew");
1564   kalmanFitOrig->Write("kalmanFitOrig");
1565   f.Close();
1566 }
1567
1568
1569
1570 /*
1571   Example - submit alignment as batch jobs
1572   ifile=0;
1573   ntracksSkip=200
1574   ntracksSkipOffset=0
1575   nclustersSkip=2
1576   nclustersSkipOffset=0
1577   ntracks=1000000
1578   ndump=5
1579   isTest=0
1580   bDir=`pwd`
1581   ver=aaa;
1582
1583   while [ $ntracksSkipOffset -lt $ntracksSkip ] ; do
1584     nclustersSkipOffset=0;
1585     while [ $nclustersSkipOffset -lt $nclustersSkip ] ; do
1586       echo Tr $ntracksSkipOffset  Cl $nclustersSkipOffset;
1587       ver=kalmanDirTrack$ntracksSkipOffset$nclustersSkipOffset
1588       echo New Directory $ver
1589       mkdir $ver;
1590       cd $ver;
1591       cp $bDir/align.txt . ;
1592       ln -sf $bDir/kalmanFitApply.root . ;
1593       echo  aliroot  -q -b  "$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C($ntracks,10000,0,$ndump,$ntracksSkip,$ntracksSkipOffset, $nclustersSkip,$nclustersSkipOffset,$isTest)" ;
1594       bsub -q alice-t3_8h -o `pwd`/output.log command aliroot  -q -b  "$ALICE_ROOT/TPC/CalibMacros/CalibAlignKalman.C($ntracks,10000,0,$ndump,$ntracksSkip,$ntracksSkipOffset, $nclustersSkip,$nclustersSkipOffset,$isTest)"  ;
1595       nclustersSkipOffset=$(( $nclustersSkipOffset + 1 ))
1596       cd $bDir;
1597       echo $bDir;
1598     done;
1599     cd $bDir;
1600     echo $bDir;
1601     ntracksSkipOffset=$(( $ntracksSkipOffset + 1 ))
1602     echo Tr $ntracksSkipOffset;
1603 done
1604
1605
1606
1607 */
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618