TPC/AliTPCcalibCosmic.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /*
  17     Comments to be written here:
  18     1. What do we calibrate.
  19     2. How to interpret results
  20     3. Simple example
  21     4. Analysis using debug streamers.
  22
  23
  24
  25     3.Simple example
  26     // To make cosmic scan the user interaction neccessary
  27     //
  28     .x ~/UliStyle.C
  29     gSystem->Load("libANALYSIS");
  30     gSystem->Load("libTPCcalib");
  31     TFile fcalib("CalibObjects.root");
  32     TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
  33     AliTPCcalibCosmic * cosmic = ( AliTPCcalibCosmic *)array->FindObject("cosmicTPC");
  34
  35
  36     //
  37     //
  38     gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  39     gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  40     AliXRDPROOFtoolkit tool;
  41     TChain * chainCosmic = tool.MakeChainRandom("cosmic.txt","Track0",0,10000);
  42     TChain * chainBudget = tool.MakeChainRandom("cosmic.txt","material",0,1000);
  43     TCut cutptV="abs(1/pt0V-1/pt1V)<0.1";
  44     TCut cutptI="abs(1/pt0In-1/pt1In)<0.5";
  45     TCut cutncl="nclmin>120";
  46     TCut cutDz="abs(p0.fP[1])<50";
  47     TCut cutDr="abs(p0.fP[0])<50";
  48     //
  49     chainBudget->Draw(">>listB",cutptV+cutptI+cutncl+cutDr+cutDz,"entryList");
  50     TEntryList *elistB = (TEntryList*)gDirectory->Get("listB");
  51     chainBudget->SetEntryList(elistB);
  52
  53     chainBudget->SetAlias("dptrel","(pt0V-pt1V)/((pt0V+pt1V)*0.5)");
  54     chainBudget->SetAlias("dptInrel","(pt0In-pt1In)/((pt0In+pt1In)*0.5)");
  55     chainBudget->SetAlias("ptcorr","(pt0In-pt0V)/(pt0V)+(pt1V-pt1In)/(pt1In)");
  56 */
  57
  58
  59
  60 #include "Riostream.h"
  61 #include "TChain.h"
  62 #include "TTree.h"
  63 #include "TH1F.h"
  64 #include "TH2F.h"
  65 #include "TList.h"
  66 #include "TMath.h"
  67 #include "TCanvas.h"
  68 #include "TFile.h"
  69 #include "TF1.h"
  70 #include "THnSparse.h"
  71 #include "TDatabasePDG.h"
  72
  73 #include "AliTPCclusterMI.h"
  74 #include "AliTPCseed.h"
  75 #include "AliESDVertex.h"
  76 #include "AliESDEvent.h"
  77 #include "AliESDfriend.h"
  78 #include "AliESDInputHandler.h"
  79 #include "AliAnalysisManager.h"
  80
  81 #include "AliTracker.h"
  82 #include "AliMagF.h"
  83 #include "AliTPCCalROC.h"
  84
  85 #include "AliLog.h"
  86
  87 #include "AliTPCcalibCosmic.h"
  88 #include "TTreeStream.h"
  89 #include "AliTPCTracklet.h"
  90 //#include "AliESDcosmic.h"
  91
  92
  93 ClassImp(AliTPCcalibCosmic)
  94
  95
  96 AliTPCcalibCosmic::AliTPCcalibCosmic()
  97   :AliTPCcalibBase(),
  98    fHistNTracks(0),
  99    fClusters(0),
 100    fModules(0),
 101    fHistPt(0),
 102    fDeDx(0),
 103    fDeDxMIP(0),
 104    fMIPvalue(1),
 105    fCutMaxD(5),        // maximal distance in rfi ditection
 106    fCutMaxDz(40),      // maximal distance in z ditection
 107    fCutTheta(0.03),    // maximal distan theta
 108    fCutMinDir(-0.99)   // direction vector products
 109 {
 110   AliInfo("Default Constructor");
 111   for (Int_t ihis=0; ihis<6;ihis++){
 112     fHistoDelta[ihis]=0;
 113     fHistoPull[ihis]=0;
 114   }
 115   for (Int_t ihis=0; ihis<4;ihis++){
 116     fHistodEdxMax[ihis]    =0;
 117     fHistodEdxTot[ihis]    =0;
 118   }
 119 }
 120
 121
 122 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
 123   :AliTPCcalibBase(),
 124    fHistNTracks(0),
 125    fClusters(0),
 126    fModules(0),
 127    fHistPt(0),
 128    fDeDx(0),
 129    fDeDxMIP(0),
 130    fMIPvalue(1),
 131    fCutMaxD(5),        // maximal distance in rfi ditection
 132    fCutMaxDz(40),      // maximal distance in z ditection
 133    fCutTheta(0.03),    // maximal distan theta
 134    fCutMinDir(-0.99)   // direction vector products
 135 {
 136   SetName(name);
 137   SetTitle(title);
 138
 139   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
 140   fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
 141   fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
 142   fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
 143   fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
 144   BinLogX(fDeDx);
 145   fDeDxMIP =  new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
 146   Init();
 147   AliInfo("Non Default Constructor");
 148   //
 149 }
 150
 151 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
 152   //
 153   //
 154   //
 155   for (Int_t ihis=0; ihis<6;ihis++){
 156     delete fHistoDelta[ihis];
 157     delete fHistoPull[ihis];
 158   }
 159   for (Int_t ihis=0; ihis<4;ihis++){
 160     delete fHistodEdxTot[ihis];
 161     delete fHistodEdxMax[ihis];
 162   }
 163
 164   delete fHistNTracks;            //  histogram showing number of ESD tracks per event
 165   delete fClusters;               //  histogram showing the number of clusters per track
 166   delete fModules;                //  2d histogram of tracks which are propagated to the ACORDE scintillator array
 167   delete fHistPt;                 //  Pt histogram of reconstructed tracks
 168   delete fDeDx;                   //  dEdx spectrum showing the different particle types
 169   delete fDeDxMIP;                //  TPC signal close to the MIP region of muons 0.4 < p < 0.45 GeV
 170 }
 171
 172
 173 void AliTPCcalibCosmic::Init(){
 174   //
 175   // init component
 176   // Make performance histograms
 177   //
 178
 179   // tracking performance bins
 180   // 0 - delta of interest
 181   // 1 - min (track0, track1) number of clusters
 182   // 2 - R  - vertex radius
 183   // 3 - P1 - mean z
 184   // 4 - P2 - snp(phi)    at inner wall of TPC
 185   // 5 - P3 - tan(theta)  at inner wall of TPC
 186   // 6 - P4 - 1/pt mean
 187   // 7 - pt - pt mean
 188   // 8 - alpha
 189
 190   Double_t xminTrack[9], xmaxTrack[9];
 191   Int_t binsTrack[9];
 192   TString axisName[9];
 193   //
 194   binsTrack[0] =100;
 195   axisName[0]  ="#Delta";
 196   //
 197   binsTrack[1] =8;
 198   xminTrack[1] =80; xmaxTrack[1]=160;
 199   axisName[1]  ="N_{cl}";
 200   //
 201   binsTrack[2] =10;
 202   xminTrack[2] =0; xmaxTrack[2]=90;  //
 203   axisName[2]  ="dca_{r} (cm)";
 204   //
 205   binsTrack[3] =25;
 206   xminTrack[3] =-250; xmaxTrack[3]=250;  //
 207   axisName[3]  ="z (cm)";
 208   //
 209   binsTrack[4] =10;
 210   xminTrack[4] =-0.8; xmaxTrack[4]=0.8;  //
 211   axisName[4]  ="sin(#phi)";
 212   //
 213   binsTrack[5] =10;
 214   xminTrack[5] =-1; xmaxTrack[5]=1;  //
 215   axisName[5]  ="tan(#theta)";
 216   //
 217   binsTrack[6] =40;
 218   xminTrack[6] =-2; xmaxTrack[6]=2;  //
 219   axisName[6]  ="1/pt (1/GeV)";
 220   //
 221   binsTrack[7] =40;
 222   xminTrack[7] =0.2; xmaxTrack[7]=50;  //
 223   axisName[7]  ="pt (GeV)";
 224   //
 225   binsTrack[8] =18;
 226   xminTrack[8] =0; xmaxTrack[8]=TMath::Pi();  //
 227   axisName[8]  ="alpha";
 228   //
 229   // delta y
 230   xminTrack[0] =-1; xmaxTrack[0]=1;  //
 231   fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
 232   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 233   fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 234   //
 235   // delta z
 236   xminTrack[0] =-1; xmaxTrack[0]=1;  //
 237   fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
 238   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 239   fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 240   //
 241   // delta P2
 242   xminTrack[0] =-10; xmaxTrack[0]=10;  //
 243   fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
 244   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 245   fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 246   //
 247   // delta P3
 248   xminTrack[0] =-10; xmaxTrack[0]=10;  //
 249   fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
 250   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 251   fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 252   //
 253   // delta P4
 254   xminTrack[0] =-0.2; xmaxTrack[0]=0.2;  //
 255   fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", 9, binsTrack,xminTrack, xmaxTrack);
 256   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 257   fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 258
 259   //
 260   // delta Pt
 261   xminTrack[0] =-0.5; xmaxTrack[0]=0.5;  //
 262   fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", 9, binsTrack,xminTrack, xmaxTrack);
 263   xminTrack[0] =-5; xmaxTrack[0]=5;  //
 264   fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
 265   //
 266
 267   for (Int_t idedx=0;idedx<4;idedx++){
 268     xminTrack[0] =0.5; xmaxTrack[0]=1.5;  //
 269     binsTrack[1] =40;
 270     xminTrack[1] =10; xmaxTrack[1]=160;
 271
 272     fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
 273                                           Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
 274                                           9, binsTrack,xminTrack, xmaxTrack);
 275     fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
 276                                           Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
 277                                           9, binsTrack,xminTrack, xmaxTrack);
 278   }
 279
 280
 281
 282   for (Int_t ivar=0;ivar<6;ivar++){
 283     for (Int_t ivar2=0;ivar2<9;ivar2++){
 284       fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 285       fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 286       fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 287       fHistoPull[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 288       BinLogX(fHistoDelta[ivar],7);
 289       BinLogX(fHistoPull[ivar],7);
 290       if (ivar<4){
 291         fHistodEdxMax[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 292         fHistodEdxMax[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 293         fHistodEdxTot[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
 294         fHistodEdxTot[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
 295         BinLogX(fHistodEdxMax[ivar],7);
 296         BinLogX(fHistodEdxTot[ivar],7);
 297       }
 298     }
 299   }
 300 }
 301
 302
 303 void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){
 304   //
 305   //
 306   //
 307   for (Int_t ivar=0; ivar<6;ivar++){
 308     if (fHistoDelta[ivar] && cosmic->fHistoDelta[ivar]){
 309       fHistoDelta[ivar]->Add(cosmic->fHistoDelta[ivar]);
 310     }
 311     if (fHistoPull[ivar] && cosmic->fHistoPull[ivar]){
 312       fHistoPull[ivar]->Add(cosmic->fHistoPull[ivar]);
 313     }
 314   }
 315   for (Int_t ivar=0; ivar<4;ivar++){
 316     if (fHistodEdxMax[ivar] && cosmic->fHistodEdxMax[ivar]){
 317       fHistodEdxMax[ivar]->Add(cosmic->fHistodEdxMax[ivar]);
 318     }
 319     if (fHistodEdxTot[ivar] && cosmic->fHistodEdxTot[ivar]){
 320       fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]);
 321     }
 322   }
 323 }
 324
 325
 326
 327
 328 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
 329   //
 330   //
 331   //
 332   if (!event) {
 333     Printf("ERROR: ESD not available");
 334     return;
 335   }
 336   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 337   if (!esdFriend) {
 338    Printf("ERROR: esdFriend not available");
 339    return;
 340   }
 341
 342
 343   FindPairs(event); // nearly everything takes place in find pairs...
 344
 345   if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
 346   Int_t ntracks=event->GetNumberOfTracks();
 347   fHistNTracks->Fill(ntracks);
 348   if (ntracks==0) return;
 349  //  AliESDcosmic cosmicESD;
 350 //   TTreeSRedirector * cstream =  GetDebugStreamer();
 351 //   cosmicESD.SetDebugStreamer(cstream);
 352 //   cosmicESD.ProcessEvent(event);
 353 //   if (cstream) cosmicESD.DumpToTree();
 354
 355
 356 }
 357
 358
 359 void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam */*inner1*/, AliTPCseed *seed0,  AliTPCseed *seed1, const AliExternalTrackParam *param0Combined ){
 360   //
 361   // par0,par1       - parameter of tracks at DCA 0
 362   // inner0,inner1   - parameter of tracks at the TPC entrance
 363   // seed0, seed1    - detailed track information
 364   // param0Combined  - Use combined track parameters for binning
 365   //
 366   Int_t kMinCldEdx =20;
 367   Int_t ncl0 = seed0->GetNumberOfClusters();
 368   Int_t ncl1 = seed1->GetNumberOfClusters();
 369
 370   const Double_t kpullCut    = 10;
 371   Double_t x[9];
 372   Double_t xyz0[3];
 373   Double_t xyz1[3];
 374   par0->GetXYZ(xyz0);
 375   par1->GetXYZ(xyz1);
 376   Double_t radius0 = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
 377   Double_t radius1 = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
 378   inner0->GetXYZ(xyz0);
 379   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
 380   // bin parameters
 381   x[1] = TMath::Min(ncl0,ncl1);
 382   x[2] = (radius0+radius1)*0.5;
 383   x[3] = param0Combined->GetZ();
 384   x[4] = inner0->GetSnp();
 385   x[5] = param0Combined->GetTgl();
 386   x[6] = param0Combined->GetSigned1Pt();
 387   x[7] = param0Combined->Pt();
 388   x[8] = alpha;
 389   // deltas
 390   Double_t delta[6];
 391   Double_t sigma[6];
 392   delta[0] = (par0->GetY()+par1->GetY());
 393   delta[1] = (par0->GetZ()-par1->GetZ());
 394   delta[2] = (par0->GetAlpha()-par1->GetAlpha()-TMath::Pi());
 395   delta[3] = (par0->GetTgl()+par1->GetTgl());
 396   delta[4] = (par0->GetParameter()[4]+par1->GetParameter()[4]);
 397   delta[5] = (par0->Pt()-par1->Pt())/((par0->Pt()+par1->Pt())*0.5);
 398   //
 399   sigma[0] = TMath::Sqrt(par0->GetSigmaY2()+par1->GetSigmaY2());
 400   sigma[1] = TMath::Sqrt(par0->GetSigmaZ2()+par1->GetSigmaZ2());
 401   sigma[2] = TMath::Sqrt(par0->GetSigmaSnp2()+par1->GetSigmaSnp2());
 402   sigma[3] = TMath::Sqrt(par0->GetSigmaTgl2()+par1->GetSigmaTgl2());
 403   sigma[4] = TMath::Sqrt(par0->GetSigma1Pt2()+par1->GetSigma1Pt2());
 404   sigma[5] = sigma[4]*((par0->Pt()+par1->Pt())*0.5);
 405   //
 406   Bool_t isOK = kTRUE;
 407   for (Int_t ivar=0;ivar<6;ivar++){
 408     if (sigma[ivar]==0) isOK=kFALSE;
 409     x[0]= delta[ivar]/sigma[ivar];
 410     if (TMath::Abs(x[0])>kpullCut) isOK = kFALSE;
 411   }
 412   //
 413
 414   if (isOK) for (Int_t ivar=0;ivar<6;ivar++){
 415     x[0]= delta[ivar]/TMath::Sqrt(2);
 416     if (ivar==2 || ivar ==3) x[0]*=1000;
 417     fHistoDelta[ivar]->Fill(x);
 418     if (sigma[ivar]>0){
 419       x[0]= delta[ivar]/sigma[ivar];
 420       fHistoPull[ivar]->Fill(x);
 421     }
 422   }
 423
 424   //
 425   // Fill dedx performance
 426   //
 427   for (Int_t ipad=0; ipad<4;ipad++){
 428     //
 429     //
 430     //
 431     Int_t row0=0;
 432     Int_t row1=160;
 433     if (ipad==0) row1=63;
 434     if (ipad==1) {row0=63; row1=63+64;}
 435     if (ipad==2) {row0=128;}
 436     Int_t   nclUp       = TMath::Nint(seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
 437     Int_t   nclDown     = TMath::Nint(seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
 438     Int_t   minCl       = TMath::Min(nclUp,nclDown);
 439     if (minCl<kMinCldEdx) continue;
 440     x[1] = minCl;
 441     //
 442     Float_t dEdxTotUp   = seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
 443     Float_t dEdxTotDown = seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
 444     Float_t dEdxMaxUp   = seed0->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
 445     Float_t dEdxMaxDown = seed1->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
 446     //
 447     if (dEdxTotDown<=0) continue;
 448     if (dEdxMaxDown<=0) continue;
 449     x[0]=dEdxTotUp/dEdxTotDown;
 450     fHistodEdxTot[ipad]->Fill(x);
 451     x[0]=dEdxMaxUp/dEdxMaxDown;
 452     fHistodEdxMax[ipad]->Fill(x);
 453   }
 454
 455
 456
 457 }
 458
 459 void AliTPCcalibCosmic::MaterialBudgetDump(AliExternalTrackParam *const par0, AliExternalTrackParam *const par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam *inner1, AliTPCseed *const seed0,  AliTPCseed *const seed1, AliExternalTrackParam *const param0Combined, AliExternalTrackParam *const param1Combined){
 460   //
 461   // matrial budget AOD dump
 462   //
 463   // par0,par1       - parameter of tracks at DCA 0
 464   // inner0,inner1   - parameter of tracks at the TPC entrance
 465   // seed0, seed1    - detailed track information
 466   // param0Combined  - Use combined track parameters for binning
 467   // param1Combined  -
 468   Double_t p0In = inner0->GetP();
 469   Double_t p1In = inner1->GetP();
 470   Double_t p0V  = par0->GetP();
 471   Double_t p1V  = par1->GetP();
 472   //
 473   Double_t pt0In = inner0->Pt();
 474   Double_t pt1In = inner1->Pt();
 475   Double_t pt0V  = par0->Pt();
 476   Double_t pt1V  = par1->Pt();
 477   Int_t ncl0 = seed0->GetNumberOfClusters();
 478   Int_t ncl1 = seed1->GetNumberOfClusters();
 479   Int_t nclmin=TMath::Min(ncl0,ncl1);
 480   Double_t sign = (param0Combined->GetSigned1Pt()>0) ? 1:-1.;
 481   //
 482   TTreeSRedirector * pcstream =  GetDebugStreamer();
 483   if (pcstream){
 484     (*pcstream)<<"material"<<
 485       "run="<<fRun<<              //  run number
 486       "event="<<fEvent<<          //  event number
 487       "time="<<fTime<<            //  time stamp of event
 488       "trigger="<<fTrigger<<      //  trigger
 489       "triggerClass="<<&fTriggerClass<<      //  trigger
 490       "mag="<<fMagF<<             //  magnetic field
 491       "sign="<<sign<<             // sign of the track
 492       //
 493       "ncl0="<<ncl0<<
 494       "ncl1="<<ncl1<<
 495       "nclmin="<<nclmin<<
 496       //
 497       "p0In="<<p0In<<
 498       "p1In="<<p1In<<
 499       "p0V="<<p0V<<
 500       "p1V="<<p1V<<
 501       "pt0In="<<pt0In<<
 502       "pt1In="<<pt1In<<
 503       "pt0V="<<pt0V<<
 504       "pt1V="<<pt1V<<
 505       "p0.="<<par0<<
 506       "p1.="<<par1<<
 507       "up0.="<<param0Combined<<
 508       "up1.="<<param1Combined<<
 509       "\n";
 510   }
 511
 512 }
 513
 514 void AliTPCcalibCosmic::Analyze() {
 515
 516   fMIPvalue = CalculateMIPvalue(fDeDxMIP);
 517
 518   return;
 519
 520 }
 521
 522
 523
 524 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
 525   //
 526   // Find cosmic pairs
 527   //
 528   // Track0 is choosen in upper TPC part
 529   // Track1 is choosen in lower TPC part
 530   //
 531   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
 532   AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
 533   Int_t ntracks=event->GetNumberOfTracks();
 534   TObjArray  tpcSeeds(ntracks);
 535   if (ntracks==0) return;
 536   Double_t vtxx[3]={0,0,0};
 537   Double_t svtxx[3]={0.000001,0.000001,100.};
 538   AliESDVertex vtx(vtxx,svtxx);
 539   //
 540   //track loop
 541   //
 542   for (Int_t i=0;i<ntracks;++i) {
 543    AliESDtrack *track = event->GetTrack(i);
 544    fClusters->Fill(track->GetTPCNcls());
 545
 546    const AliExternalTrackParam * trackIn = track->GetInnerParam();
 547    const AliExternalTrackParam * trackOut = track->GetOuterParam();
 548    if (!trackIn) continue;
 549    if (!trackOut) continue;
 550    if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue;  // filter laser
 551
 552
 553    AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
 554    TObject *calibObject;
 555    AliTPCseed *seed = 0;
 556    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
 557      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
 558    }
 559    if (seed) tpcSeeds.AddAt(seed,i);
 560
 561    Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
 562    if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
 563      fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159));
 564      //
 565      if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159));
 566      //
 567      if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) {
 568        TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
 569        if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
 570        if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
 571      }
 572
 573    }
 574
 575   }
 576
 577   if (ntracks<2) return;
 578   //
 579   // Find pairs
 580   //
 581   for (Int_t i=0;i<ntracks;++i) {
 582     AliESDtrack *track0 = event->GetTrack(i);
 583     // track0 - choosen upper part
 584     if (!track0) continue;
 585     if (!track0->GetOuterParam()) continue;
 586     if (track0->GetOuterParam()->GetAlpha()<0) continue;
 587     Double_t dir0[3];
 588     track0->GetDirection(dir0);
 589     for (Int_t j=0;j<ntracks;++j) {
 590       if (i==j) continue;
 591       AliESDtrack *track1 = event->GetTrack(j);
 592       //track 1 lower part
 593       if (!track1) continue;
 594       if (!track1->GetOuterParam()) continue;
 595       if (track1->GetOuterParam()->GetAlpha()>0) continue;
 596       //
 597       Double_t dir1[3];
 598       track1->GetDirection(dir1);
 599
 600       AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
 601       AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
 602       if (! seed0) continue;
 603       if (! seed1) continue;
 604       Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159);
 605       Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159);
 606       //
 607       Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
 608       Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
 609       //
 610       Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159);
 611       Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159);
 612       //
 613       Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
 614       Float_t d0  = track0->GetLinearD(0,0);
 615       Float_t d1  = track1->GetLinearD(0,0);
 616       //
 617       // conservative cuts - convergence to be guarantied
 618       // applying before track propagation
 619       if (TMath::Abs(d0+d1)>fCutMaxD) continue;   // distance to the 0,0
 620       if (dir>fCutMinDir) continue;               // direction vector product
 621       Float_t bz = AliTracker::GetBz();
 622       Float_t dvertex0[2];   //distance to 0,0
 623       Float_t dvertex1[2];   //distance to 0,0
 624       track0->GetDZ(0,0,0,bz,dvertex0);
 625       track1->GetDZ(0,0,0,bz,dvertex1);
 626       if (TMath::Abs(dvertex0[1])>250) continue;
 627       if (TMath::Abs(dvertex1[1])>250) continue;
 628       //
 629       //
 630       //
 631       Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
 632       AliExternalTrackParam param0(*track0);
 633       AliExternalTrackParam param1(*track1);
 634       //
 635       // Propagate using Magnetic field and correct fo material budget
 636       //
 637       Double_t sign0=-1;
 638       Double_t sign1=1;
 639       Double_t maxsnp=0.90;
 640       AliTracker::PropagateTrackToBxByBz(&param0,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign0);
 641       AliTracker::PropagateTrackToBxByBz(&param1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign1);
 642       //
 643       // Propagate rest to the 0,0 DCA - z should be ignored
 644       //
 645       Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
 646       Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
 647       //
 648       param0.GetDZ(0,0,0,bz,dvertex0);
 649       param1.GetDZ(0,0,0,bz,dvertex1);
 650       if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
 651       //
 652       Double_t xyz0[3];//,pxyz0[3];
 653       Double_t xyz1[3];//,pxyz1[3];
 654       param0.GetXYZ(xyz0);
 655       param1.GetXYZ(xyz1);
 656       Bool_t isPair = IsPair(&param0,&param1);
 657       //
 658       if (isPair) FillAcordeHist(track0);
 659       //
 660       // combined track params
 661       //
 662       AliExternalTrackParam *par0U=MakeCombinedTrack(&param0,&param1);
 663       AliExternalTrackParam *par1U=MakeCombinedTrack(&param1,&param0);
 664
 665
 666       //
 667       if (fStreamLevel>0){
 668         TTreeSRedirector * cstream =  GetDebugStreamer();
 669         //printf("My stream=%p\n",(void*)cstream);
 670         AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
 671         AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
 672         AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
 673         AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
 674         Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
 675         Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
 676         Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
 677         Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
 678         //
 679         //
 680         //
 681         FillHistoPerformance(&param0, &param1, ip0, ip1, seed0, seed1,par0U);
 682         MaterialBudgetDump(&param0, &param1, ip0, ip1, seed0, seed1,par0U,par1U);
 683         if (cstream) {
 684           (*cstream) << "Track0" <<
 685             "run="<<fRun<<              //  run number
 686             "event="<<fEvent<<          //  event number
 687             "time="<<fTime<<            //  time stamp of event
 688             "trigger="<<fTrigger<<      //  trigger
 689             "triggerClass="<<&fTriggerClass<<      //  trigger
 690             "mag="<<fMagF<<             //  magnetic field
 691             "dir="<<dir<<               //  direction
 692             "OK="<<isPair<<             //  will be accepted
 693             "b0="<<b0<<                 //  propagate status
 694             "b1="<<b1<<                 //  propagate status
 695             "crossI="<<isCrossI<<       //  cross inner
 696             "crossO="<<isCrossO<<       //  cross outer
 697             //
 698             "Orig0.=" << track0 <<      //  original track  0
 699             "Orig1.=" << track1 <<      //  original track  1
 700             "Tr0.="<<&param0<<          //  track propagated to the DCA 0,0
 701             "Tr1.="<<&param1<<          //  track propagated to the DCA 0,0
 702             "Ip0.="<<ip0<<              //  inner param - upper
 703             "Ip1.="<<ip1<<              //  inner param - lower
 704             "Op0.="<<op0<<              //  outer param - upper
 705             "Op1.="<<op1<<              //  outer param - lower
 706             "Up0.="<<par0U<<           //  combined track 0
 707             "Up1.="<<par1U<<           //  combined track 1
 708             //
 709             "v00="<<dvertex0[0]<<       //  distance using kalman
 710             "v01="<<dvertex0[1]<<       //
 711             "v10="<<dvertex1[0]<<       //
 712             "v11="<<dvertex1[1]<<       //
 713             "d0="<<d0<<                 //  linear distance to 0,0
 714             "d1="<<d1<<                 //  linear distance to 0,0
 715             //
 716             //
 717             //
 718             "x00="<<xyz0[0]<<           // global position close to vertex
 719             "x01="<<xyz0[1]<<
 720             "x02="<<xyz0[2]<<
 721             //
 722             "x10="<<xyz1[0]<<           // global position close to vertex
 723             "x11="<<xyz1[1]<<
 724             "x12="<<xyz1[2]<<
 725             //
 726             "alpha0="<<alpha0<<
 727             "alpha1="<<alpha1<<
 728             "dir00="<<dir0[0]<<           // direction upper
 729             "dir01="<<dir0[1]<<
 730             "dir02="<<dir0[2]<<
 731             //
 732             "dir10="<<dir1[0]<<           // direction lower
 733             "dir11="<<dir1[1]<<
 734             "dir12="<<dir1[2]<<
 735             //
 736             //
 737             "Seed0.=" << seed0 <<       //  original seed 0
 738             "Seed1.=" << seed1 <<       //  original seed 1
 739             //
 740             "dedx0="<<dedx0<<           //  dedx0 - all
 741             "dedx1="<<dedx1<<           //  dedx1 - all
 742             //
 743             "dedx0I="<<dedx0I<<         //  dedx0 - inner ROC
 744             "dedx1I="<<dedx1I<<         //  dedx1 - inner ROC
 745             //
 746             "dedx0O="<<dedx0O<<         //  dedx0 - outer ROC
 747             "dedx1O="<<dedx1O<<         //  dedx1 - outer ROC
 748             "\n";
 749         }
 750       }
 751       delete par0U;
 752       delete par1U;
 753     }
 754   }
 755 }
 756
 757
 758
 759
 760 void  AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
 761
 762   // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
 763   if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
 764
 765   const Double_t acordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
 766   const Double_t roof = 210.5;       // distance from x =0 to end of magnet roof
 767
 768   Double_t r[3];
 769   upperTrack->GetXYZ(r);
 770   Double_t d[3];
 771   upperTrack->GetDirection(d);
 772   Double_t x,z;
 773   z = r[2] + (d[2]/d[1])*(acordePlane - r[1]);
 774   x = r[0] + (d[0]/d[1])*(acordePlane - r[1]);
 775
 776   if (x > roof) {
 777     x = r[0] + (d[0]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
 778     z = r[2] + (d[2]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
 779   }
 780   if (x < -roof) {
 781     x = r[0] + (d[0]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
 782     z = r[2] + (d[2]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
 783   }
 784
 785   fModules->Fill(z, x);
 786
 787 }
 788
 789
 790
 791 Long64_t AliTPCcalibCosmic::Merge(TCollection *const li) {
 792
 793   TIterator* iter = li->MakeIterator();
 794   AliTPCcalibCosmic* cal = 0;
 795
 796   while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
 797     if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
 798       //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
 799       return -1;
 800     }
 801
 802     fHistNTracks->Add(cal->GetHistNTracks());
 803     fClusters->Add(cal-> GetHistClusters());
 804     fModules->Add(cal->GetHistAcorde());
 805     fHistPt->Add(cal->GetHistPt());
 806     fDeDx->Add(cal->GetHistDeDx());
 807     fDeDxMIP->Add(cal->GetHistMIP());
 808     Add(cal);
 809   }
 810   return 0;
 811
 812 }
 813
 814
 815 Bool_t  AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
 816   //
 817   //
 818   /*
 819   // 0. Same direction - OPOSITE  - cutDir +cutT
 820   TCut cutDir("cutDir","dir<-0.99")
 821   // 1.
 822   TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
 823   //
 824   // 2. The same rphi
 825   TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
 826   //
 827   //
 828   //
 829   TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
 830   // 1/Pt diff cut
 831   */
 832   const Double_t *p0 = tr0->GetParameter();
 833   const Double_t *p1 = tr1->GetParameter();
 834   if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
 835   if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
 836   if (TMath::Abs(p0[0]+p1[0])>fCutMaxD)  return kFALSE;
 837
 838   Double_t d0[3], d1[3];
 839   tr0->GetDirection(d0);
 840   tr1->GetDirection(d1);
 841   if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
 842   //
 843   return kTRUE;
 844 }
 845
 846
 847
 848 Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
 849
 850   TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
 851   funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
 852                                 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
 853   hist->Fit(funcDoubleGaus);
 854   Double_t aMIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
 855
 856   delete funcDoubleGaus;
 857
 858   return aMIPvalue;
 859
 860 }
 861
 862
 863
 864
 865 void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
 866   //
 867   // Not implemented yet
 868   //
 869   return;
 870
 871 }
 872
 873
 874 void AliTPCcalibCosmic::BinLogX(THnSparse *const h, Int_t axisDim) {
 875
 876   // Method for the correct logarithmic binning of histograms
 877
 878   TAxis *axis = h->GetAxis(axisDim);
 879   int bins = axis->GetNbins();
 880
 881   Double_t from = axis->GetXmin();
 882   Double_t to = axis->GetXmax();
 883   Double_t *newBins = new Double_t[bins + 1];
 884
 885   newBins[0] = from;
 886   Double_t factor = pow(to/from, 1./bins);
 887
 888   for (int i = 1; i <= bins; i++) {
 889    newBins[i] = factor * newBins[i-1];
 890   }
 891   axis->Set(bins, newBins);
 892   delete newBins;
 893
 894 }
 895
 896
 897 void AliTPCcalibCosmic::BinLogX(TH1 *const h) {
 898
 899   // Method for the correct logarithmic binning of histograms
 900
 901   TAxis *axis = h->GetXaxis();
 902   int bins = axis->GetNbins();
 903
 904   Double_t from = axis->GetXmin();
 905   Double_t to = axis->GetXmax();
 906   Double_t *newBins = new Double_t[bins + 1];
 907
 908   newBins[0] = from;
 909   Double_t factor = pow(to/from, 1./bins);
 910
 911   for (int i = 1; i <= bins; i++) {
 912    newBins[i] = factor * newBins[i-1];
 913   }
 914   axis->Set(bins, newBins);
 915   delete newBins;
 916
 917 }
 918
 919
 920 AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
 921   //
 922   //
 923   //
 924   AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
 925   par1R->Rotate(track0->GetAlpha());
 926   par1R->PropagateTo(track0->GetX(),AliTracker::GetBz());
 927   //
 928   //
 929   Double_t * param = (Double_t*)par1R->GetParameter();
 930   Double_t * covar = (Double_t*)par1R->GetCovariance();
 931
 932   param[0]*=1;  //OK
 933   param[1]*=1;  //OK
 934   param[2]*=1;  //?
 935   param[3]*=-1; //OK
 936   param[4]*=-1; //OK
 937   //
 938   covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
 939   //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
 940   covar[13]*=-1.;
 941   return par1R;
 942 }
 943
 944 AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
 945   //
 946   // Make combined track
 947   //
 948   //
 949   AliExternalTrackParam * par1T = MakeTrack(track0,track1);
 950   AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
 951   //
 952   UpdateTrack(*par0U,*par1T);
 953   delete par1T;
 954   return par0U;
 955 }
 956
 957
 958 void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
 959   //
 960   // Update track 1 with track 2
 961   //
 962   //
 963   //
 964   TMatrixD vecXk(5,1);    // X vector
 965   TMatrixD covXk(5,5);    // X covariance
 966   TMatrixD matHk(5,5);    // vector to mesurement
 967   TMatrixD measR(5,5);    // measurement error
 968   TMatrixD vecZk(5,1);    // measurement
 969   //
 970   TMatrixD vecYk(5,1);    // Innovation or measurement residual
 971   TMatrixD matHkT(5,5);
 972   TMatrixD matSk(5,5);    // Innovation (or residual) covariance
 973   TMatrixD matKk(5,5);    // Optimal Kalman gain
 974   TMatrixD mat1(5,5);     // update covariance matrix
 975   TMatrixD covXk2(5,5);   //
 976   TMatrixD covOut(5,5);
 977   //
 978   Double_t *param1=(Double_t*) track1.GetParameter();
 979   Double_t *covar1=(Double_t*) track1.GetCovariance();
 980   Double_t *param2=(Double_t*) track2.GetParameter();
 981   Double_t *covar2=(Double_t*) track2.GetCovariance();
 982   //
 983   // copy data to the matrix
 984   for (Int_t ipar=0; ipar<5; ipar++){
 985     for (Int_t jpar=0; jpar<5; jpar++){
 986       covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
 987       measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
 988       matHk(ipar,jpar)=0;
 989       mat1(ipar,jpar)=0;
 990     }
 991     vecXk(ipar,0) = param1[ipar];
 992     vecZk(ipar,0) = param2[ipar];
 993     matHk(ipar,ipar)=1;
 994     mat1(ipar,ipar)=0;
 995   }
 996   //
 997   //
 998   //
 999   //
1000   //
1001   vecYk = vecZk-matHk*vecXk;                 // Innovation or measurement residual
1002   matHkT=matHk.T(); matHk.T();
1003   matSk = (matHk*(covXk*matHkT))+measR;      // Innovation (or residual) covariance
1004   matSk.Invert();
1005   matKk = (covXk*matHkT)*matSk;              //  Optimal Kalman gain
1006   vecXk += matKk*vecYk;                      //  updated vector
1007   covXk2 = (mat1-(matKk*matHk));
1008   covOut =  covXk2*covXk;
1009   //
1010   //
1011   //
1012   // copy from matrix to parameters
1013   if (0) {
1014     vecXk.Print();
1015     vecZk.Print();
1016     //
1017     measR.Print();
1018     covXk.Print();
1019     covOut.Print();
1020     //
1021     track1.Print();
1022     track2.Print();
1023   }
1024
1025   for (Int_t ipar=0; ipar<5; ipar++){
1026     param1[ipar]= vecXk(ipar,0) ;
1027     for (Int_t jpar=0; jpar<5; jpar++){
1028       covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
1029     }
1030   }
1031 }
1032
1033
1034
1035