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 // AliTPCTransformation::BuildBasicFormulas();
   9
  10 // AliXRDPROOFtoolkit tool;
  11 // chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0,50000);
  12 // chainPoints->Lookup();
  13
  14 // CalibAlignKalman(40000);
  15 // kalmanFit0->DumpCorelation(0.8);
  16 // TFile f("kalmanfitTPC.root");
  17
  18
  19 */
  20 #include <fstream>
  21
  22 #include "TSystem.h"
  23 #include "TROOT.h"
  24 #include "TRandom.h"
  25 #include "TMath.h"
  26 #include "TBits.h"
  27 #include "TFormula.h"
  28 #include "TF1.h"
  29 #include "TLinearFitter.h"
  30 #include "TFile.h"
  31 #include "TChain.h"
  32 #include "TCut.h"
  33 #include "TEntryList.h"
  34 #include "TH1F.h"
  35
  36
  37 #include "TTreeStream.h"
  38 #include "AliTrackPointArray.h"
  39 #include "AliLog.h"
  40 #include "AliTPCTransformation.h"
  41 #include "AliTPCkalmanFit.h"
  42 #include "AliXRDPROOFtoolkit.h"
  43
  44
  45 //
  46 TChain *chainPoints=0;
  47 TEntryList *elist=0;
  48 AliTPCkalmanFit * kalmanFit0=0;
  49 AliTPCkalmanFit * kalmanFitIdeal=0;
  50
  51 AliTPCkalmanFit *  CalibAlignKalmanFit(Int_t maxTracks);
  52 void FilterTracks();
  53 AliTPCkalmanFit * SetupFit();
  54 void              AddFitFieldCage(AliTPCkalmanFit *kalmanFit);
  55 void              AddPhiScaling(AliTPCkalmanFit *kalmanFit);
  56 void              AddZShift(AliTPCkalmanFit *kalmanFit);
  57 void              AddZRotation(AliTPCkalmanFit *kalmanFit);
  58 void              AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit);
  59 AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector *pcstream);
  60 AliTPCkalmanFit *   FitPointsLinear(Int_t maxTracks);
  61
  62 void CalibAlignKalman(Int_t npoints, Int_t maxFiles, Int_t startFile){
  63   //
  64   //
  65   //
  66   gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  67   gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+");
  68
  69   AliTPCTransformation::BuildBasicFormulas();
  70   AliXRDPROOFtoolkit tool;
  71   chainPoints = tool.MakeChainRandom("align.txt","trackPoints",0, maxFiles, startFile);
  72   chainPoints->Lookup();
  73   CalibAlignKalmanFit(npoints);
  74 }
  75
  76
  77
  78 //
  79 //
  80 // track point cuts
  81 //
  82 const Float_t krmsYcut    = 0.3;
  83 const Float_t krmsZcut    = 0.3;
  84 const Float_t kSigmaCut   = 10.;
  85 const Int_t   knclCut     =  80;
  86 const Double_t kArmCut    = 50.;
  87 //
  88 // Track cuts
  89 //
  90 TCut cutsP3="abs(p0Out.fP[3]-p0In.fP[3])<0.0002&&abs(p1Out.fP[3]-p1In.fP[3])<0.0002";
  91 TCut cutsP4="abs(p0Out.fP[4]-p0In.fP[4])<0.005&&abs(p1Out.fP[4]-p1In.fP[4])<0.005";
  92 TCut cutP3 ="abs(p1In.fP[3]+p0In.fP[3]-0.0015)<0.0025";
  93 TCut cutP4 ="abs(p1In.fP[4])<0.03&&abs(p0In.fP[4])<0.03";
  94 TCut cutSide="p1Out.fP[1]*p0Out.fP[1]>0&&p1In.fP[1]*p0In.fP[1]>0";
  95 TCut cutAll=cutSide+cutsP4+cutsP3+cutP3+cutP4+"abs(mag)<0.01&&ncont>0&&p.fNPoints>120";
  96
  97
  98
  99
 100 AliTPCkalmanFit *  CalibAlignKalmanFit(Int_t maxTracks){
 101   //
 102   //
 103   AliTPCTransformation::BuildBasicFormulas();
 104   FilterTracks();
 105   kalmanFit0     = SetupFit();
 106   kalmanFitIdeal = SetupFit();
 107   kalmanFit0->Init();
 108   kalmanFitIdeal->Init();
 109   return FitPointsLinear(maxTracks);
 110 }
 111
 112
 113
 114 AliTPCkalmanFit * SetupFit(){
 115   //
 116   AliTPCkalmanFit *kalmanFit =  new AliTPCkalmanFit;
 117   AddFitFieldCage(kalmanFit);
 118   AddPhiScaling(kalmanFit);
 119   AddZShift(kalmanFit);
 120   AddZRotation(kalmanFit);
 121   AddLocalXYMisalignment(kalmanFit);
 122   return kalmanFit;
 123 }
 124
 125
 126 void FilterTracks(){
 127   //
 128   //
 129   //
 130   chainPoints->Draw(">>listEL",cutAll,"entryList");
 131   elist = (TEntryList*)gDirectory->Get("listEL");
 132   chainPoints->SetEntryList(elist);
 133   elist->SetDirectory(0);
 134 }
 135
 136
 137
 138 AliTPCkalmanFit * FitPointsLinear(Int_t maxTracks){
 139   //
 140   //
 141   //
 142   // create debug streeamers
 143   TTreeSRedirector *pcstream      = new TTreeSRedirector("kalmanfitTPC.root");
 144   TTreeSRedirector *pcstreamIdeal = new TTreeSRedirector("kalmanfitTPCOrig.root");
 145   //
 146   //
 147   AliTrackPointArray *points=0;
 148   Float_t mag=0;
 149   Int_t   time=0;
 150   chainPoints->SetBranchAddress("p.",&points);
 151   chainPoints->SetBranchAddress("mag",&mag);
 152   chainPoints->SetBranchAddress("time",&time);
 153   Int_t accepted=0;
 154   //
 155   for (Int_t itrack=0;itrack<elist->GetN(); itrack++){
 156     Int_t entry=chainPoints->GetEntryNumber(itrack);
 157     chainPoints->GetEntry(entry);
 158     if (TMath::Abs(mag)>0.01) {printf("mag- Cut not accpeted\n"); continue;}
 159     if (points->GetNPoints()<knclCut) {printf("ncl - Cut not accpeted\n"); continue;}
 160     if (accepted>maxTracks) break;
 161     for (Int_t idir=-1; idir<=1; idir++){
 162       AliTrackPointArray *pointsF = FilterPoints(*points,idir, pcstream);
 163       if (!pointsF) continue;
 164       if (idir==0)  accepted++;
 165       //
 166       if (accepted%50==0) {
 167         kalmanFit0->FitTrackLinear(*pointsF, 10, pcstream);
 168       }else{
 169         kalmanFit0->FitTrackLinear(*pointsF, 1, 0);
 170       }
 171       if (idir==0) kalmanFit0->DumpTrackLinear(*pointsF,pcstream);
 172       if (idir==0) kalmanFitIdeal->DumpTrackLinear(*pointsF,pcstreamIdeal);
 173       if (accepted%25==0) printf("%d\n", accepted);
 174       delete pointsF;
 175     }
 176   }
 177   pcstream->GetFile()->cd();
 178   kalmanFit0->Write("kalmanFit");
 179   pcstreamIdeal->GetFile()->cd();
 180   kalmanFitIdeal->Write("kalmanFitIdeal");
 181   delete pcstream;
 182   delete pcstreamIdeal;
 183   return kalmanFit0;
 184 }
 185
 186
 187
 188
 189 AliTrackPointArray *FilterPoints(AliTrackPointArray &points, Int_t dir, TTreeSRedirector */*pcstream*/){
 190   //
 191   //
 192   //
 193   TLinearFitter lfitY(2,"pol1");
 194   TLinearFitter lfitZ(2,"pol1");
 195   TVectorD vecZ(2);
 196   TVectorD vecY(2);
 197   //
 198   lfitY.StoreData(kTRUE);
 199   lfitZ.StoreData(kTRUE);
 200   Int_t npoints = points.GetNPoints();
 201   if (npoints<2) return 0;
 202   Double_t currentAlpha = TMath::ATan2(points.GetY()[npoints-1]-points.GetY()[0], points.GetX()[npoints-1]-points.GetX()[0]);
 203   Double_t ca = TMath::Cos(currentAlpha);
 204   Double_t sa = TMath::Sin(currentAlpha);
 205   //
 206   // 1.b Fit the track in the rotated frame - MakeSeed
 207   //
 208   Double_t maxX =-10000, minX=10000;
 209   for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 210     Double_t rx =   ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
 211     Double_t ry =  -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
 212     Double_t rz =  points.GetZ()[ipoint];
 213     if (dir== 1 && rx<0) continue;
 214     if (dir==-1 && rx>0) continue;
 215     if (maxX<rx) maxX=rx;
 216     if (minX>rx) minX=rx;
 217     lfitY.AddPoint(&rx,ry,1);
 218     lfitZ.AddPoint(&rx,rz,1);
 219   }
 220   if (TMath::Abs(maxX-minX)<kArmCut) return 0;
 221   if (lfitY.GetNpoints()<knclCut) return 0;
 222   //
 223   lfitY.Eval();
 224   lfitZ.Eval();
 225   lfitY.GetParameters(vecY);
 226   lfitZ.GetParameters(vecZ);
 227   //
 228   Double_t chi2Y = lfitY.GetChisquare()/lfitY.GetNpoints();
 229   Double_t chi2Z = lfitZ.GetChisquare()/lfitZ.GetNpoints();
 230   if (TMath::Sqrt(chi2Y)>krmsYcut) return 0;
 231   if (TMath::Sqrt(chi2Z)>krmsZcut) return 0;
 232   //
 233   //
 234   Int_t accepted=0;
 235   AliTrackPoint point;
 236   AliTrackPointArray *pointsF=0;
 237   for (Int_t iter=0; iter<2;iter++){
 238     for (Int_t ipoint=0; ipoint<npoints-1; ipoint++){
 239       //
 240       if (!points.GetPoint(point,ipoint)) continue;
 241       Double_t rx =   ca*points.GetX()[ipoint]+sa*points.GetY()[ipoint];
 242       Double_t ry =  -sa*points.GetX()[ipoint]+ca*points.GetY()[ipoint];
 243       Double_t rz =  points.GetZ()[ipoint];
 244       if (dir== 1 && rx<0) continue;
 245       if (dir==-1 && rx>0) continue;
 246       Double_t erry = TMath::Sqrt(chi2Y);
 247       Double_t errz = TMath::Sqrt(chi2Z);
 248       Double_t fy = vecY[0]+vecY[1]*rx;
 249       Double_t fz = vecZ[0]+vecZ[1]*rx;
 250       if (TMath::Abs(fy-ry)>erry*kSigmaCut) continue;
 251       if (TMath::Abs(fz-rz)>errz*kSigmaCut) continue;
 252       if (pointsF) pointsF->AddPoint(accepted,&point);
 253       accepted++;
 254     }
 255     if (accepted<knclCut) break;
 256     if (iter==0) pointsF = new AliTrackPointArray(accepted);
 257     accepted=0;
 258   }
 259   return pointsF;
 260 }
 261
 262
 263
 264
 265 void  AddFitFieldCage(AliTPCkalmanFit *kalmanFit){
 266   //
 267   // Add radial scaling due field cage
 268   //
 269   TVectorD fpar(10);
 270   AliTPCTransformation * transformation=0;
 271   char tname[100];
 272   //
 273   // linear R scaling and shift
 274   //
 275   for (Int_t iside=0; iside<=1; iside++)
 276     for (Int_t ipolR=0; ipolR<2; ipolR++){
 277       for (Int_t ipolZ=0; ipolZ<3; ipolZ++){
 278         fpar[0]=ipolR;
 279         fpar[1]=ipolZ;
 280         if (ipolR+ipolZ==0) continue;
 281         sprintf(tname,"tTPCscalingRPolR%dDr%dSide%d",ipolR,ipolZ,iside);
 282         transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRPol",0,0,  1);
 283         transformation->SetParams(0,0.2,0,&fpar);
 284         kalmanFit->AddCalibration(transformation);
 285         //
 286       }
 287     }
 288   //
 289   //
 290   //Inner field cage
 291   for (Int_t iside=0; iside<=1; iside++)
 292     for (Int_t ipol=0; ipol<3; ipol++){
 293       fpar[0]=ipol;
 294       sprintf(tname,"tTPCscalingRIFC%dSide%d",ipol,iside);
 295       transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingRIFC",0,0,   1);
 296       transformation->SetParams(0,0.2,0,&fpar);
 297       kalmanFit->AddCalibration(transformation);
 298     }
 299   //
 300   //
 301   //Outer field cage
 302   for (Int_t iside=0; iside<=1; iside++)
 303     for (Int_t ipol=0; ipol<3; ipol++){
 304       fpar[0]=ipol;
 305       //Outer field cage
 306       sprintf(tname,"tTPCscalingROFC%dSide%d",ipol,iside);
 307       transformation = new AliTPCTransformation(tname,AliTPCTransformation::BitsSide(iside),"TPCscalingROFC",0,0,  1);
 308       transformation->SetParams(0,0.2,0,&fpar);
 309       kalmanFit->AddCalibration(transformation);
 310     }
 311 }
 312
 313
 314 void AddPhiScaling(AliTPCkalmanFit *kalmanFit){
 315   //
 316   // Add local phi scaling
 317   //
 318   TVectorD fpar(10);
 319   AliTPCTransformation * transformation=0;
 320   fpar[0]=1;
 321   transformation = new AliTPCTransformation("tscalingLocalPhi", AliTPCTransformation::BitsAll(), 0,"TPCscalingPhiLocal",0,  1);
 322   transformation->SetParams(0,0.2,0,&fpar);
 323   kalmanFit->AddCalibration(transformation);
 324   //
 325 }
 326
 327 void  AddZShift(AliTPCkalmanFit *kalmanFit){
 328   TVectorD fpar(10);
 329   AliTPCTransformation * transformation=0;
 330   TBits maskInnerA(72);
 331   TBits maskInnerC(72);
 332   TBits maskOuter(72);
 333   for (Int_t i=0;i<72;i++){
 334     if (i<36){
 335       if (i%36<18)  maskInnerA[i]=kTRUE;
 336       if (i%36>=18) maskInnerC[i]=kTRUE;
 337     }
 338     if (i>=36){
 339       maskOuter[i]=kTRUE;
 340     }
 341   }
 342   //
 343   transformation = new AliTPCTransformation("tTPCDeltaZIROCA", new TBits(maskInnerA), 0,0, "TPCDeltaZ",  0);
 344   transformation->SetParams(0,0.2,0,&fpar);
 345   kalmanFit->AddCalibration(transformation);
 346   transformation = new AliTPCTransformation("tTPCDeltaZIROCC", new TBits(maskInnerC), 0,0, "TPCDeltaZ",  0);
 347   transformation->SetParams(0,0.2,0,&fpar);
 348   kalmanFit->AddCalibration(transformation);
 349   //
 350   transformation = new AliTPCTransformation("tTPCDeltaZOROCML", new TBits(maskOuter), 0,0, "TPCDeltaZMediumLong",  0);
 351   transformation->SetParams(0,0.2,0,&fpar);
 352   kalmanFit->AddCalibration(transformation);
 353   //
 354 }
 355
 356 void AddZRotation(AliTPCkalmanFit *kalmanFit){
 357   //
 358   //
 359   //
 360   TVectorD fpar(10);
 361   //  AliTPCTransformation::BitsSide(iside)
 362   fpar[0]=0; fpar[1]=0;
 363   AliTPCTransformation * transformation=0;
 364   transformation = new AliTPCTransformation("tTPCTiltingZAside00",AliTPCTransformation::BitsSide(0) , 0,0, "TPCTiltingZ",  0);
 365   transformation->SetParams(0,0.4,0,&fpar);
 366   kalmanFit->AddCalibration(transformation);
 367   transformation = new AliTPCTransformation("tTPCTiltingZCside00",AliTPCTransformation::BitsSide(1) , 0,0, "TPCTiltingZ",  0);
 368   transformation->SetParams(0,0.4,0,&fpar);
 369   kalmanFit->AddCalibration(transformation);
 370   //
 371   fpar[0]=1; fpar[1]=0;
 372   transformation = new AliTPCTransformation("tTPCTiltingZAside10",AliTPCTransformation::BitsSide(0) , 0,0, "TPCTiltingZ",  0);
 373   transformation->SetParams(0,0.4,0,&fpar);
 374   kalmanFit->AddCalibration(transformation);
 375   transformation = new AliTPCTransformation("tTPCTiltingZCside10",AliTPCTransformation::BitsSide(1) , 0,0, "TPCTiltingZ",  0);
 376   transformation->SetParams(0,0.4,0,&fpar);
 377   kalmanFit->AddCalibration(transformation);
 378   //
 379   //
 380   fpar[0]=0; fpar[1]=1;
 381   transformation = new AliTPCTransformation("tTPCTiltingZAside01",AliTPCTransformation::BitsSide(0) , 0,0, "TPCTiltingZ",  0);
 382   transformation->SetParams(0,0.4,0,&fpar);
 383   kalmanFit->AddCalibration(transformation);
 384   transformation = new AliTPCTransformation("tTPCTiltingZCside01",AliTPCTransformation::BitsSide(1) , 0,0, "TPCTiltingZ",  0);
 385   transformation->SetParams(0,0.4,0,&fpar);
 386   kalmanFit->AddCalibration(transformation);
 387 }
 388
 389
 390
 391 void  AddLocalXYMisalignment(AliTPCkalmanFit *kalmanFit){
 392   //
 393   //
 394   //
 395   TVectorD fpar(10);
 396   AliTPCTransformation * transformation=0;
 397   TBits maskInnerA(72);
 398   TBits maskInnerC(72);
 399   for (Int_t i=0;i<72;i++){
 400     if (i<36){
 401       if (i%36<18)  maskInnerA[i]=kTRUE;
 402       if (i%36>=18) maskInnerC[i]=kTRUE;
 403     }
 404   }
 405   //
 406   //
 407   transformation = new AliTPCTransformation("tTPCDeltaLxIROCA", new TBits(maskInnerA), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
 408   transformation->SetParams(0,0.2,0,&fpar);
 409   kalmanFit->AddCalibration(transformation);
 410   transformation = new AliTPCTransformation("tTPCDeltaLxIROCC", new TBits(maskInnerC), "TPClocaldLxdGX","TPClocaldLxdGY",0,  0);
 411   transformation->SetParams(0,0.2,0,&fpar);
 412   kalmanFit->AddCalibration(transformation);
 413   //
 414   transformation = new AliTPCTransformation("tTPCDeltaLyIROCA", new TBits(maskInnerA), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
 415   transformation->SetParams(0,0.2,0,&fpar);
 416   kalmanFit->AddCalibration(transformation);
 417   transformation = new AliTPCTransformation("tTPCDeltaLyIROCC", new TBits(maskInnerC), "TPClocaldLydGX","TPClocaldLydGY",0,  0);
 418   transformation->SetParams(0,0.2,0,&fpar);
 419   kalmanFit->AddCalibration(transformation);
 420 }
 421
 422 void MergeKalman(const char * list = "kalmanFit.list"){
 423   //
 424   //
 425   //
 426   ifstream in;
 427   in.open(list);
 428   TString currentFile;
 429   kalmanFit0= 0;
 430   Int_t counter=0;
 431   while(in.good()) {
 432     in >> currentFile;
 433     printf("%d\t%d\t%s\n", counter,currentFile.Length(),currentFile.Data());
 434     if (currentFile.Length()==0) continue;
 435     TFile * ffit = TFile::Open(currentFile.Data());
 436     AliTPCkalmanFit * fit = ( AliTPCkalmanFit *)ffit->Get("kalmanFit");
 437     if (!fit) continue;
 438     if (!kalmanFit0) {kalmanFit0= fit; continue;};
 439     kalmanFit0->Add(fit);
 440     delete fit;
 441     delete ffit;
 442     counter++;
 443   }
 444 }
 445
 446
 447 /*
 448   myvar=0;
 449   ntracks=10000000
 450   bDir=`pwd`
 451   while [ $myvar -ne 190 ] ; do mkdir kalmanDir$myvar; cd kalmanDir$myvar; cp $bDir/align.txt .;  bsub -q proof command aliroot  -q -b  "CalibAlignKalman.C($ntracks,5,$myvar)" ; myvar=$(( $myvar + 5 )) ; echo $myvar ; cd $bDir; echo $bDir; done
 452
 453 */
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465