TPC/CalibMacros/CalibAlignKalman.C

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