HLT/ITS/trackingSAP/AliITSSAPTracker.cxx

   1 #include "AliITSSAPTracker.h"
   2 #include "AliITSSAPLayer.h"
   3 #include "AliITSRecPoint.h"
   4 #include "AliGeomManager.h"
   5 #include "AliVParticle.h"
   6 #include "AliSymMatrix.h"
   7 //
   8 #include "AliRunLoader.h"
   9 #include "AliHeader.h"
  10 #include "AliGenEventHeader.h"
  11 #include "AliStack.h"
  12 #include <TParticle.h>
  13 #include <TParticlePDG.h>
  14 #include <TFile.h>
  15
  16
  17 ClassImp(AliITSSAPTracker)
  18
  19 const Float_t AliITSSAPTracker::fgkZSpanITS[AliITSSAPTracker::kMaxLrITS] =
  20 { 36. ,14.1,14.1,  38., 22.2,29.7, 51.   ,43.1,48.9};
  21
  22 const Float_t AliITSSAPTracker::fgkRLayITS[AliITSSAPTracker::kMaxLrITS] =
  23   { 2.94, 3.9,7.6, 11.04, 15.0,23.9, 29.44 ,38.0,43.0};
  24
  25 const Float_t AliITSSAPTracker::fgkRSpanITS[AliITSSAPTracker::kMaxLrITS] = // half span in R
  26   { 0.04, 0.5,0.5, 0.5, 0.8, 0.8, 0.5 ,0.6,0.6};
  27
  28 const Int_t    AliITSSAPTracker::fgkPassivLrITS[AliITSSAPTracker::kNLrPassive] =
  29   {AliITSSAPTracker::kLrBeamPime,AliITSSAPTracker::kLrShield1,AliITSSAPTracker::kLrShield2};
  30
  31 const Int_t    AliITSSAPTracker::fgkActiveLrITS[AliITSSAPTracker::kNLrActive] =
  32   {AliITSSAPTracker::kLrSPD1,AliITSSAPTracker::kLrSPD2,
  33    AliITSSAPTracker::kLrSDD1,AliITSSAPTracker::kLrSDD2,
  34    AliITSSAPTracker::kLrSSD1,AliITSSAPTracker::kLrSSD2};
  35
  36 const Int_t    AliITSSAPTracker::fgkLr2Active[AliITSSAPTracker::kMaxLrITS] = // conversion to active lr.
  37   {-1, 0, 1, -1, 2, 3, -1, 4, 5};
  38
  39 const Float_t AliITSSAPTracker::fgkRhoLITS[AliITSSAPTracker::kMaxLrITS] = {
  40   0.162802, 0.321960,0.354588, 0.274995, 0.193789,0.198168, 0.435372, 0.195828,0.226940};
  41
  42 const Float_t AliITSSAPTracker::fgkX2X0ITS[AliITSSAPTracker::kMaxLrITS] = {
  43   0.002757, 0.011660,0.012614, 0.006488, 0.007714,0.007916, 0.012689, 0.007849,0.009128};
  44
  45
  46 const Double_t AliITSSAPTracker::fgkClSystYErr2[AliITSSAPTracker::kNLrActive] =
  47   {0.0010*0.0010, 0.0030*0.0030, 0.0500*0.0500, 0.0500*0.0500, 0.0020*0.0020, 0.0020*0.0020};
  48
  49 const Double_t AliITSSAPTracker::fgkClSystZErr2[AliITSSAPTracker::kNLrActive] =
  50   {0.0050*0.0050, 0.0050*0.0050, 0.0050*0.0050, 0.0050*0.0050, 0.1000*0.1000, 0.1000*0.1000};
  51
  52
  53 const Int_t    AliITSSAPTracker::fgkLrDefBins[AliITSSAPTracker::kNLrActive][2] = // n bins in z, phi
  54   { {20,20}, {20,20}, {20,20}, {20,20}, {20,20}, {20,20} };
  55
  56 const Float_t AliITSSAPTracker::fgkDefMass = 0.14;
  57 const Int_t   AliITSSAPTracker::fgkDummyLabel = -3141593;
  58
  59 #ifdef _TIMING_
  60 const char* AliITSSAPTracker::fgkSWNames[AliITSSAPTracker::kNSW] = {
  61   "Total"
  62   ,"Tracklets"
  63   ,"Tracks"
  64   ,"Vertex"
  65 };
  66 #endif
  67
  68
  69 //______________________________________________
  70 AliITSSAPTracker::AliITSSAPTracker() :
  71   fBlacklist(0)
  72   ,fPhiShift(0.0045)
  73   ,fSigThetaTracklet(0.025)
  74   ,fSigPhiTracklet(0.08)
  75   ,fChi2CutTracklet(1.5)
  76   ,fPhiShiftSc(0.)
  77   ,fDThetaTrackletSc(0)
  78   ,fDPhiTrackletSc(0)
  79   ,fBz(5.0)
  80   ,fDPhiTol(0.)
  81   ,fDThSig2Inv(0.)
  82   ,fDPhSig2Inv(0.)
  83   //
  84   ,fMinPt(0.3)
  85   ,fCurvMax(0)
  86   ,fZSPD2CutMin(1e9)
  87   ,fZSPD2CutMax(-1e9)
  88   ,fMaxChi2Tr2Cl(40)
  89   ,fAddErr2YspdVtx(0.02*0.02)
  90   ,fAddErr2ZspdVtx(0.04*0.04)
  91   //
  92   ,fMissChi2Penalty(3)
  93   ,fMaxMissedLayers(1)
  94   ,fNTracks(0)
  95   ,fFitVertex(kTRUE)
  96   //
  97   ,fSPDVertex(0)
  98 #ifdef _CONTROLH_
  99   ,fHTrackletMC(0),fHTrackletAll(0),fHTrackletFake(0),fHTrackMC(0),fHTrackAll(0),fHTrackFake(0)
 100   ,fHVtxDiffXY(0)
 101   ,fHVtxDiffXMlt(0),fHVtxDiffYMlt(0),fHVtxDiffZMlt(0)
 102   ,fHVtxPullXMlt(0),fHVtxPullYMlt(0),fHVtxPullZMlt(0)
 103   ,fHVtxMCSPDDiffXY(0)
 104   ,fHVtxMCSPDDiffXMlt(0),fHVtxMCSPDDiffYMlt(0),fHVtxMCSPDDiffZMlt(0)
 105   ,fHVtxMCSPDPullXMlt(0),fHVtxMCSPDPullYMlt(0),fHVtxMCSPDPullZMlt(0)
 106   ,fHChi2NDFvsPT(0),fHChi2vsNC(0)
 107   ,fHVtxMltRef(0),fHVtxOKMlt(0),fHVtxDiffZ(0),fHVtxMCSPDDiffZ(0)
 108 #endif
 109 {
 110   // def. c-tor
 111   for (int i=kNLrActive;i--;) fLayers[i] = 0;
 112 }
 113
 114 //______________________________________________
 115 AliITSSAPTracker::~AliITSSAPTracker()
 116 {
 117   // d-tor
 118   for (int i=0;i<kNLrActive;i++) delete fLayers[i];
 119 }
 120
 121 //______________________________________________
 122 void AliITSSAPTracker::Init()
 123 {
 124   // init tracker
 125   //
 126   if (!AliGeomManager::GetGeometry()) {
 127     AliGeomManager::LoadGeometry("geometry.root");
 128     AliGeomManager::ApplyAlignObjsFromCDB("ITS");
 129   }
 130   //
 131   for (int i=0;i<kNLrActive;i++) {
 132     int iAct = fgkActiveLrITS[i];
 133     fLayers[i] = new AliITSSAPLayer(i,fgkZSpanITS[iAct]+1,fgkLrDefBins[i][0],fgkLrDefBins[i][1]);
 134     fSkipLayer[i] = kFALSE;
 135     fNSigma2[i] = 7*7;
 136     fYToler2[i] = 0.2*0.2;
 137     fZToler2[i] = 0.2*0.2;
 138     fChi2TotCut[i] = 0;
 139   }
 140   fChi2TotCut[1] = 40; // 2 cl+vtx -> NDF=1
 141   fChi2TotCut[2] = 40;
 142   fChi2TotCut[3] = 30;
 143   fChi2TotCut[4] = 35;
 144   fChi2TotCut[5] = 40;
 145   //
 146   fMissChi2Penalty = 3;
 147   fMaxMissedLayers = 1;
 148   //
 149   // auxialary precalculated variables
 150   if (fChi2CutTracklet<0.1) fChi2CutTracklet = 0.1;
 151   double scl = TMath::Sqrt(fChi2CutTracklet);
 152   fDThetaTrackletSc = fSigThetaTracklet*scl;
 153   fDPhiTrackletSc   = fSigPhiTracklet*scl;
 154   //
 155   fDThSig2Inv = 1./(fSigThetaTracklet*fSigThetaTracklet);
 156   fDPhSig2Inv = 1./(fSigPhiTracklet*fSigPhiTracklet);
 157   //
 158   fBlacklist = new TBits(100*100);
 159   //
 160 #ifdef _TIMING_
 161   for (int i=kNSW;i--;) {
 162     fSW[i].Stop();
 163     fSW[i].Reset();
 164   }
 165 #endif
 166   //
 167 #ifdef _CONTROLH_
 168   BookHistos();
 169 #endif
 170 }
 171
 172 //______________________________________________
 173 void AliITSSAPTracker::ProcessEvent()
 174 {
 175   // do full reconstruction
 176 #ifdef _TIMING_
 177   fSW[kSWTotal].Start(0);
 178   fSW[kSWTracklets].Start(0);
 179 #endif
 180   //
 181   fNTracks = 0;
 182   FindTracklets();
 183   //
 184 #ifdef _TIMING_
 185   fSW[kSWTracklets].Stop();
 186   fSW[kSWTracks].Start(0);
 187 #endif
 188   //
 189   Tracklets2Tracks();
 190   RefitInward();
 191 #ifdef _TIMING_
 192   fSW[kSWTracks].Stop();
 193   fSW[kSWVertex].Start(0);
 194 #endif
 195   if (fFitVertex) {
 196     if (FitTrackVertex()) {
 197 #ifdef _DEBUG_
 198       printf("FittedVertex: "); fTrackVertex.Print();
 199       printf("SPD   Vertex: "); fSPDVertex->Print();
 200 #endif
 201     }
 202   }
 203   //
 204   //
 205 #ifdef _TIMING_
 206   fSW[kSWVertex].Stop();
 207   fSW[kSWTotal].Stop();
 208   PrintTiming();
 209 #endif
 210   //
 211 #ifdef _CONTROLH_
 212   FillRecoStat();
 213 #endif
 214 }
 215
 216
 217 //______________________________________________
 218 void AliITSSAPTracker::Clear(Option_t*)
 219 {
 220   // reset event info
 221   ClearTracklets();
 222   ClearTracks();
 223   for (int i=kNLrActive;i--;) {
 224     fNClusters[i] = 0;
 225     if (fLayers[i]) fLayers[i]->Clear();
 226   }
 227 }
 228
 229 //______________________________________________
 230 void AliITSSAPTracker::ClearTracklets()
 231 {
 232   // reset tracklets info
 233   fSPD2Discard.clear();
 234   fTracklets.clear();
 235   fSPD1Tracklet.clear();
 236   if (fBlacklist) fBlacklist->ResetAllBits();
 237 }
 238
 239
 240 //______________________________________________
 241 void AliITSSAPTracker::AddCluster(AliITSRecPoint* cl)
 242 {
 243   // add cluster to corresponding layer
 244   if (!cl->Misalign()) AliWarning("Can't misalign this cluster !");
 245   fLayers[cl->GetLayer()]->AddCluster(cl);
 246 }
 247
 248 //______________________________________________
 249 Bool_t AliITSSAPTracker::FindTracklets()
 250 {
 251   // find SPD tracklets
 252   //
 253   if (!fSPDVertex) {
 254     //    AliInfo("No SPD vertex set");
 255     return kFALSE;
 256   }
 257   float rv2 = fSPDVertex->GetX()*fSPDVertex->GetX()+fSPDVertex->GetY()*fSPDVertex->GetY();
 258   if (rv2>0.25*fgkRLayITS[kLrBeamPime]*fgkRLayITS[kLrBeamPime]) {
 259     //    AliInfo("SPD vertex is too far from beam line");
 260     fSPDVertex->Print();
 261     return kFALSE;
 262   }
 263   fPhiShiftSc = fPhiShift*TMath::Abs(fBz/5.0);
 264   fDPhiTol = fDPhiTrackletSc + fPhiShiftSc;
 265   //
 266   AliITSSAPLayer &spdL1 = *fLayers[kALrSPD1];
 267   AliITSSAPLayer &spdL2 = *fLayers[kALrSPD2];
 268   spdL1.SortClusters(fSPDVertex);
 269   spdL2.SortClusters(fSPDVertex);
 270   fNClusters[0] = spdL1.GetNClusters();
 271   fNClusters[1] = spdL2.GetNClusters();
 272   //
 273   if (fNClusters[0]<1 || fNClusters[1]<1) return kFALSE;
 274   //
 275   fSPD2Discard.resize(fNClusters[1]);
 276   fSPD1Tracklet.resize(fNClusters[0]);
 277   //
 278   fBlacklist->SetBitNumber(TMath::Max(fNClusters[0]*fNClusters[1],10000),kFALSE); // to reserve the space
 279   //
 280   int nfound;
 281   do {
 282     nfound = 0;
 283     for (int icl2=fNClusters[1];icl2--;) if (!fSPD2Discard[icl2]) nfound += AssociateClusterOfL2(icl2);
 284   } while(nfound);
 285   //
 286   for (int itr=GetNTracklets();itr--;) CookLabel(fTracklets[itr]);
 287   //
 288   return kTRUE;
 289 }
 290
 291 //______________________________________________
 292 Int_t AliITSSAPTracker::AssociateClusterOfL2(int icl2)
 293 {
 294   // find SPD1 cluster matching to SPD2 cluster icl2
 295   AliITSSAPLayer &spdL1 = *fLayers[kALrSPD1];
 296   AliITSSAPLayer &spdL2 = *fLayers[kALrSPD2];
 297   AliITSSAPLayer::ClsInfo* cli2 = spdL2.GetClusterInfo(icl2);
 298   // expected z at SPD1
 299   float zV = fSPDVertex->GetZ();
 300   float z2 = cli2->z - zV;
 301   float tg2Inv = z2/cli2->r;
 302   float dzt = (1.+tg2Inv*tg2Inv)*fDThetaTrackletSc;
 303   float dz = dzt*fgkRLayITS[kLrSPD1] + TMath::Abs(tg2Inv)*fgkRSpanITS[kLrSPD1]; // uncertainty from dTheta and from Layer1 R spread
 304   float zL1 = zV + tg2Inv*fgkRLayITS[kLrSPD1]; // center of expected Z1
 305   int nsel1 = spdL1.SelectClusters(zL1-dz,zL1+dz, cli2->phi-fDPhiTol,cli2->phi+fDPhiTol);
 306   if (!nsel1) {
 307     fSPD2Discard[icl2] = true;
 308     return 0; // no candidates
 309   }
 310   float chiBest = 9999;
 311   SPDtracklet_t trk;
 312   trk.id1 = -1;
 313   int icl1,nCand=0;
 314   while ( (icl1=spdL1.GetNextClusterInfoID())!=-1) {  // loop over matching clusters of lr1
 315     if (IsBlacklisted(icl1,icl2)) continue;
 316     AliITSSAPLayer::ClsInfo* cli1 = spdL1.GetClusterInfo(icl1);
 317     float z1 = cli1->z - zV;
 318     float tg1Inv = z1/cli1->r;
 319     //
 320     float dTheta = (tg2Inv-tg1Inv)/(1.+tg1Inv*tg1Inv);        // fast check on theta
 321     if (TMath::Abs(dTheta)>fDThetaTrackletSc) continue;
 322     //
 323     float dPhi = cli1->phi - cli2->phi;                       // fast check on phi
 324     if (dPhi>TMath::Pi()) dPhi = TMath::TwoPi()-dPhi;
 325     double dPhiS = TMath::Abs(dPhi)-fPhiShiftSc;
 326     if (TMath::Abs(dPhiS)>fDPhiTrackletSc) continue;
 327     //
 328     float chi2 = dTheta*dTheta*fDThSig2Inv + dPhiS*dPhiS*fDPhSig2Inv; // check final chi2
 329     if (chi2>1.) continue;
 330     nCand++;
 331     if (chi2>chiBest) continue;
 332     // check if cl1 is already associated with better
 333     trk.id1 = icl1;
 334     trk.id2 = icl2;
 335     trk.dtht = dTheta;
 336     trk.dphi = dPhi;
 337     trk.chi2 = chiBest = chi2;
 338   }
 339   //
 340   if (trk.id1!=-1) { // check if there is no better icl1 candidate for icl2
 341     int oldId = fSPD1Tracklet[trk.id1];
 342     if (!oldId) { // store new tracklet
 343       fTracklets.push_back(trk);
 344       fSPD1Tracklet[trk.id1] = fTracklets.size(); // refer from clusters to tracklet (id+1)
 345       Blacklist(trk.id1,trk.id2);
 346       return 1;
 347     }
 348     SPDtracklet_t& oldTrk = (SPDtracklet_t&)fSPD1Tracklet[--oldId];
 349     if (oldTrk.chi2 < trk.chi2) { // previous is better
 350       Blacklist(trk.id1,trk.id2);  // shall we blacklist new combination?
 351       if (nCand==1)  fSPD2Discard[icl2] = true; // there was just 1 candidate and it is discarded
 352       return 0;
 353     }
 354     oldTrk = trk;         // new combination is better, overwrite it with new one
 355     Blacklist(trk.id1,trk.id2);
 356     return 1;
 357   }
 358   //
 359   fSPD2Discard[icl2] = true; // no chance to find partner for this cluster
 360   return 0;
 361   //
 362 }
 363
 364
 365 //______________________________________________
 366 void AliITSSAPTracker::Tracklets2Tracks()
 367 {
 368   // try to extend tracklets to outer layers
 369   int nTrk = GetNTracklets();
 370   if (!nTrk) return;
 371   //
 372   CalcAuxTracking(); // RS??? do we need to repeat this?
 373   //
 374   for (int ila=kALrSDD1;ila<kNLrActive;ila++) {
 375     if (fSkipLayer[ila]) continue;
 376     fLayers[ila]->SortClusters(0);
 377     fNClusters[ila] = fLayers[ila]->GetNClusters();
 378   }
 379   //
 380   fTracks.resize(nTrk);
 381
 382   //
 383   for (int itr=0;itr<nTrk;itr++) {
 384     SPDtracklet_t& trlet = fTracklets[itr];
 385     //
 386 #ifdef _DEBUG_
 387     printf("TestTracklet %d\t|",itr);
 388     int stat = GetTrackletMCTruth(trlet);
 389     //
 390     int nmiss=0;
 391     for (int i=2;i<kNLrActive;i++) {
 392       printf("%c", (stat&(0x1<<i)) ? '*':'-');
 393       if (!(stat&(0x1<<i))) nmiss++;
 394     }
 395     printf("|\n");
 396     PrintTracklet(itr);
 397 #endif
 398     //
 399     float zspd2 = fLayers[kALrSPD2]->GetClusterInfo(trlet.id2)->z;
 400     if (zspd2<fZSPD2CutMin || zspd2>fZSPD2CutMax) continue;
 401     ITStrack_t &track = fTracks[fNTracks];
 402     if (!CreateTrack(track, trlet)) continue;
 403     track.trackletID = itr;
 404     Bool_t res;
 405 #ifdef _DEBUG_
 406     double xyz[3];
 407     track.paramOut.GetXYZAt(0,fBz,xyz);
 408     printf("process track pt:%f XYZ: %+.4f %+.4f %+.4f\n",track.paramOut.Pt(),xyz[0],xyz[1],xyz[2]);
 409 #endif
 410     for (int lrID=kLrShield1;lrID<kMaxLrITS;lrID++) {
 411       res = FollowToLayer(track,lrID) && IsAcceptableTrack(track);
 412       if (!res) break;
 413     }
 414 #ifdef _DEBUG_
 415     printf("%s:%d\n",res ? "OK" : "Fail",nmiss<=fMaxMissedLayers);
 416 #endif
 417     if (!res) continue;
 418     track.paramOut.ResetBit(kInvalidBit); // flag that outward fit succeeded
 419     CookLabel(track);
 420     fNTracks++;
 421     //
 422   }
 423 }
 424
 425 //______________________________________________
 426 Bool_t AliITSSAPTracker::IsAcceptableTrack(const AliITSSAPTracker::ITStrack_t& track) const
 427 {
 428   // check if the track is acceptable
 429   return kTRUE;
 430 }
 431
 432 //______________________________________________
 433 void AliITSSAPTracker::PrintTrack(const AliITSSAPTracker::ITStrack_t& track) const
 434 {
 435   // print track info
 436   printf("Chi2 = %f for %d clusters. Tracklet %d\n",track.chi2,track.ncl,track.trackletID);
 437   //
 438   for (int ilr=0;ilr<kNLrActive;ilr++) {
 439     if (track.clID[ilr]<0) continue;
 440     AliITSRecPoint* cl = fLayers[ilr]->GetClusterSorted(track.clID[ilr]);
 441     printf("L%d #%4d ",ilr,track.clID[ilr]);
 442     for (int i=0;i<3;i++) printf("%d ",cl->GetLabel(i)); printf("\n");
 443   }
 444   track.paramOut.Print();
 445   track.paramInw.Print();
 446 }
 447
 448 //______________________________________________
 449 void AliITSSAPTracker::PrintTracklets() const
 450 {
 451   // print traklets info
 452   int ntr = fTracklets.size();
 453   printf("NTracklets: %d\n",ntr);
 454   printf("Nspd1: %4d Nspd2: %4d, Ntracklets: %d\n",fNClusters[0],fNClusters[1],ntr);
 455   for (int itr=0;itr<ntr;itr++) PrintTracklet(itr);
 456   //
 457 }
 458
 459 //______________________________________________
 460 void AliITSSAPTracker::PrintTracklet(Int_t itr) const
 461 {
 462   // print single tracklet
 463   const SPDtracklet_t* trk = &fTracklets[itr];
 464   AliITSRecPoint* cl1 = fLayers[kALrSPD1]->GetClusterSorted(trk->id1);
 465   AliITSRecPoint* cl2 = fLayers[kALrSPD2]->GetClusterSorted(trk->id2);
 466   AliITSSAPLayer::ClsInfo_t* cli0 = fLayers[kALrSPD1]->GetClusterInfo(trk->id1);
 467   printf("#%3d Phi:%+.3f Eta:%+.3f Dphi:%+.3f Dtht:%+.3f Chi2:%.3f | Lbl:",
 468          itr,cli0->phi,
 469          -TMath::Log(TMath::Tan(TMath::ATan2(cli0->r,cli0->z-fSPDVertex->GetZ())/2.)),
 470          trk->dphi,trk->dtht,trk->chi2);
 471   int lab=-1,lb = -1;
 472   for (int i=0;i<3;i++) if ( (lb=cl1->GetLabel(i))>=0 ) {if (lab<0)lab=lb; printf(" %5d",lb);} printf("|");
 473   for (int i=0;i<3;i++) if ( (lb=cl2->GetLabel(i))>=0 ) printf(" %5d",lb);
 474   printf("| ->%d\n",trk->label);
 475   lab = TMath::Abs(trk->label);
 476   //
 477   AliStack* stack = 0;
 478   AliRunLoader* rl = AliRunLoader::Instance();
 479   if (lab>=0 && rl && (stack=rl->Stack())) {
 480     TParticle* mctr = stack->Particle(lab);
 481     if (mctr) {
 482       TParticlePDG* mctrPDG = mctr->GetPDG();
 483       if (mctrPDG) {
 484         double qpt = mctrPDG->Charge()>0 ? mctr->Pt() : -mctr->Pt();
 485         printf("MCTrack: Prim:%d Vxyz: {%+.4f %+.4f %+.4f} 1/pt: %.3f tgl: %.3f\n",
 486                stack->IsPhysicalPrimary(lab),
 487                mctr->Vx(),mctr->Vy(),mctr->Vz(),
 488                TMath::Abs(qpt)>0 ? 1./qpt : 9999., TMath::Tan(TMath::Pi()/2. - mctr->Theta()));
 489       }
 490     }
 491   }
 492 }
 493
 494
 495 //______________________________________________
 496 void AliITSSAPTracker::PrintTracks() const
 497 {
 498   // print tracks info
 499   printf("NTracks: %d\n",fNTracks);
 500   for (int itr=0;itr<fNTracks;itr++) PrintTrack(fTracks[itr]);
 501   //
 502 }
 503
 504
 505 //______________________________________________
 506 void AliITSSAPTracker::CalcAuxTracking()
 507 {
 508   // precalculate auxilarry variables for tracking
 509   //
 510   // largest track curvature to search
 511   const double ztolerEdge = 1.0;
 512   fCurvMax = TMath::Abs(fBz*kB2C/fMinPt);
 513   double thMin =-1e9;
 514   double thMax = 1e9;
 515   for (int ilA=kNLrActive-1;ilA>kALrSPD2;ilA--) {
 516     if (!IsObligatoryLayer(ilA)) continue;
 517     int ilr=fgkActiveLrITS[ilA];
 518     double r   = fgkRLayITS[ilr] - fgkRSpanITS[ilr];
 519     double dz = fgkZSpanITS[ilr]+ztolerEdge+fDThetaTrackletSc*r;
 520     double ri  = 1./r;
 521     double tmin= (-dz-fSPDVertex->GetZ())*ri;
 522     double tmax= ( dz-fSPDVertex->GetZ())*ri;
 523     if (tmin>thMin) thMin = tmin;
 524     if (tmax<thMax) thMax = tmax;
 525   }
 526   double r = fgkRLayITS[kLrSPD2] + fgkRSpanITS[kLrSPD2];
 527   fZSPD2CutMin = fSPDVertex->GetZ()+thMin*r; // min Z of SPD2 in tracklet to consider tracking
 528   fZSPD2CutMax = fSPDVertex->GetZ()+thMax*r; // max Z of SPD2 in tracklet to consider tracking
 529   //
 530 }
 531
 532 //______________________________________________
 533 Bool_t AliITSSAPTracker::CreateTrack(AliITSSAPTracker::ITStrack_t& track,
 534                                      AliITSSAPTracker::SPDtracklet_t& trlet)
 535 {
 536   // create track seed from tracklet
 537   // init track
 538   track.label = trlet.label;
 539   //
 540   AliITSSAPLayer::ClsInfo_t *cli1=fLayers[kALrSPD1]->GetClusterInfo(trlet.id1);
 541   AliITSSAPLayer::ClsInfo_t *cli2=fLayers[kALrSPD2]->GetClusterInfo(trlet.id2);
 542   AliITSRecPoint *cl1=fLayers[kALrSPD1]->GetClusterUnSorted(cli1->index);
 543   AliITSRecPoint *cl2=fLayers[kALrSPD2]->GetClusterUnSorted(cli2->index);
 544   int det1 = cl1->GetVolumeId()-fLayers[kALrSPD1]->GetVIDOffset();
 545   int det2 = cl2->GetVolumeId()-fLayers[kALrSPD2]->GetVIDOffset();
 546   AliITSSAPLayer::ITSDetInfo_t& detInfo1 = fLayers[kALrSPD1]->GetDetInfo(det1);
 547   AliITSSAPLayer::ITSDetInfo_t& detInfo2 = fLayers[kALrSPD2]->GetDetInfo(det2);
 548   //
 549   // crude momentun estimate
 550   float dx=cli1->x-cli2->x,dy=cli1->y-cli2->y,d=TMath::Sqrt(dx*dx+dy*dy);
 551   float qptInv = fBz ? 2*TMath::Sin(cli2->phi-cli1->phi)/d/fBz/kB2C : 0; // positive particle goes anticlockwise in B+
 552   //
 553   // we initialize the seed in the tracking frame of 1st detector
 554   float xv= fSPDVertex->GetX()*detInfo1.cosTF + fSPDVertex->GetY()*detInfo1.sinTF;
 555   float yv=-fSPDVertex->GetX()*detInfo1.sinTF + fSPDVertex->GetY()*detInfo1.cosTF;
 556   float zv= fSPDVertex->GetZ();
 557   float par[5] = {yv, zv, (float)TMath::Sin(cli1->phi-detInfo1.phiTF), (cli1->z-zv)/cli1->r, qptInv};
 558   double covVtx[6];
 559   fSPDVertex->GetCovarianceMatrix(covVtx);
 560   float cov[15] = {float(covVtx[0]+covVtx[2] + fAddErr2YspdVtx),
 561                    0, float(covVtx[5] + fAddErr2ZspdVtx),
 562                    0,0,1,
 563                    0,0,0,1,
 564                    0,0,0,0,100*100};
 565   AliExternalTrackParam& param = track.paramOut;
 566   param.Set(xv, detInfo1.phiTF, par, cov);
 567   track.chi2 = 0;   // chi2 at 1st two point is 0
 568   // go to 1st layer, ignoring the MS (errors are anyway not defined)
 569   if (!param.PropagateTo(detInfo1.xTF+cl1->GetX(), fBz)) return kFALSE;
 570   Double_t cpar0[2]={ cl1->GetY(), cl1->GetZ()};
 571   Double_t ccov0[3]={ cl1->GetSigmaY2() + GetClSystYErr2(kALrSPD1), 0., cl1->GetSigmaZ2() + GetClSystZErr2(kALrSPD1)};
 572   if (!param.Update(cpar0,ccov0)) return kFALSE;
 573   if (!param.CorrectForMeanMaterial(fgkX2X0ITS[kLrSPD1],-fgkRhoLITS[kLrSPD1],fgkDefMass)) return kFALSE;
 574   // go to 2nd layer
 575   if (!param.Rotate(detInfo2.phiTF)) return kFALSE;
 576   if (!param.PropagateTo(detInfo2.xTF+cl2->GetX(), fBz)) return kFALSE;
 577   Double_t cpar1[2]={ cl2->GetY(), cl2->GetZ()};
 578   Double_t ccov1[3]={ cl2->GetSigmaY2() + GetClSystYErr2(kALrSPD2), 0., cl2->GetSigmaZ2() + GetClSystZErr2(kALrSPD2)};
 579   track.chi2 += param.GetPredictedChi2(cpar1,ccov1);
 580   if (!param.Update(cpar1,ccov1)) return kFALSE;
 581 #ifdef _CONTROLH_
 582   FillTrackingControlHistos(1,track.label,&param,cpar1,ccov1,cl2);
 583 #endif
 584   //
 585   track.clID[0] = trlet.id1;
 586   track.clID[1] = trlet.id2;
 587   track.clID[2] = track.clID[3] = track.clID[4] = track.clID[5] = -1;
 588   track.ncl = 2;
 589   track.nmiss=0;
 590   //
 591   param.SetBit(kInvalidBit); // flag that track is not yer refitted outward
 592   track.paramOut.SetBit(kInvalidBit); // flag that track was not refitter inward
 593   return kTRUE;
 594 }
 595
 596 //______________________________________________
 597 Bool_t AliITSSAPTracker::CrossPassiveLayer(AliExternalTrackParam& param, Int_t lrID)
 598 {
 599   // cross the layer, applying mat. corrections
 600   double xStart=param.GetX();
 601   double xToGo = GetXatLabRLin(param,fgkRLayITS[lrID]);
 602   if (xToGo<0 || !param.PropagateTo(xToGo,fBz)) return kFALSE;
 603   double x2x0=fgkX2X0ITS[lrID],xrho=fgkRhoLITS[lrID];
 604   if (xStart<xToGo) xrho = -xrho; // inward propagation
 605   return param.CorrectForMeanMaterial(x2x0,xrho,fgkDefMass,kFALSE);
 606 //
 607 }
 608
 609 //______________________________________________
 610 Bool_t AliITSSAPTracker::FollowToLayer(AliITSSAPTracker::ITStrack_t& track, Int_t lrID)
 611 {
 612   // take track to given layer, searching hits if needed and applying mat. corrections
 613   int lrIDA = fgkLr2Active[lrID]; // active layer ID
 614   if (lrIDA<0 || fSkipLayer[lrIDA]) return CrossPassiveLayer(track.paramOut,lrID);
 615   //
 616   AliExternalTrackParam trCopy(track.paramOut);
 617   double xToGo = GetXatLabRLin(trCopy,fgkRLayITS[lrID]); // aproximate X at lrID
 618   if (!trCopy.PropagateTo(xToGo,fBz)) return kFALSE;
 619   double xyz[3];
 620   trCopy.GetXYZ(xyz);
 621   double phi=TMath::ATan2(xyz[1],xyz[0]),z=trCopy.GetZ();
 622   // we need track errors in the plane nearly tangential to crossing point
 623   if (!trCopy.Rotate(phi)) return kFALSE;
 624   double dphi = TMath::Sqrt(trCopy.GetSigmaY2()*fNSigma2[lrIDA]+fYToler2[lrIDA])/fgkRLayITS[lrID];
 625   double dz   = TMath::Sqrt(trCopy.GetSigmaZ2()*fNSigma2[lrIDA]+fZToler2[lrIDA]);
 626   AliITSSAPLayer* lrA = fLayers[lrIDA];
 627   int nCl = lrA->SelectClusters(z-dz,z+dz,phi-dphi,phi+dphi);
 628   Bool_t updDone = kFALSE;
 629   //
 630 #ifdef _DEBUG_
 631   printf("at Lr%d, Ncl:%d ",lrIDA,nCl);
 632   trCopy.Print();
 633 #endif
 634   //
 635   if (nCl) {
 636     int icl,iclBest=-1;
 637     double chi2Best = fMaxChi2Tr2Cl;
 638     AliITSRecPoint* bestCl = 0;
 639     AliExternalTrackParam bestTr;
 640     //
 641 #ifdef _DEBUG_
 642     int iclt=0;
 643 #endif
 644     while ( (icl=lrA->GetNextClusterInfoID())!=-1) {
 645       AliITSSAPLayer::ClsInfo_t *cli = lrA->GetClusterInfo(icl);
 646       AliITSRecPoint *cl=lrA->GetClusterUnSorted(cli->index);
 647       int detId = cl->GetVolumeId()-lrA->GetVIDOffset();
 648       AliITSSAPLayer::ITSDetInfo_t& detInfo = lrA->GetDetInfo(detId);
 649       trCopy = track.paramOut;
 650       if (!trCopy.Propagate(detInfo.phiTF, detInfo.xTF+cl->GetX(), fBz)) continue;
 651       double cpar[2]={ cl->GetY(), cl->GetZ()};
 652       double ccov[3]={ cl->GetSigmaY2() + GetClSystYErr2(lrIDA) , 0., cl->GetSigmaZ2() + GetClSystZErr2(lrIDA)};
 653       double chi2cl = trCopy.GetPredictedChi2(cpar,ccov);
 654       //
 655 #ifdef _DEBUG_
 656       float clXYZ[3]; cl->GetGlobalXYZ(clXYZ);
 657       double trXYZ[3]; trCopy.GetXYZ(trXYZ);
 658       Float_t xCl, alphaCl;
 659       cl->GetXAlphaRefPlane(xCl,alphaCl);
 660       //
 661       printf("cl%d Chi2:%.2f Dyz: %+e %+e Err: %e %e %e |Lb:",iclt++,chi2cl,
 662              cl->GetY()-trCopy.GetY(),cl->GetZ()-trCopy.GetZ(),
 663              TMath::Sqrt(ccov[0]),ccov[1],TMath::Sqrt(ccov[2])); //TMP
 664       for (int j=0;j<3;j++) if (cl->GetLabel(j)>=0) printf(" %d",cl->GetLabel(j)); printf("\n");
 665       printf("CL: X:%.4f Alp:%+.4f XYZ: %+.4f %+.4f %+.4f\n",xCl,alphaCl,clXYZ[0],clXYZ[1],clXYZ[2]);
 666       printf("TR: X:%.4f Alp:%+.4f XYZ: %+.4f %+.4f %+.4f\n",detInfo.xTF,detInfo.phiTF,trXYZ[0],trXYZ[1],trXYZ[2]);
 667       trCopy.Print();
 668 #endif
 669       //
 670 #ifdef _CONTROLH_
 671       FillTrackingControlHistos(lrIDA,track.label,&trCopy,cpar,ccov,cl);
 672 #endif
 673       //
 674       if (chi2cl>fMaxChi2Tr2Cl) continue;
 675       //    SaveCandidate(lrIDA,trCopy,chi2cl,icl);  // RS: do we need this?
 676       if (chi2cl>chi2Best) continue;
 677       chi2Best = chi2cl;
 678       iclBest = icl;
 679       bestCl = cl;
 680       bestTr = trCopy;
 681       if (nCl==1) { // in absence of competitors, do the fit on spot
 682         if (!bestTr.Update(cpar,ccov)) return kFALSE;
 683         updDone = kTRUE;
 684       }
 685     }
 686 #ifdef _DEBUG_
 687     printf("Lr%d -> %f\n",lrIDA,chi2Best);
 688 #endif
 689     //
 690     if (bestCl) {
 691       if (!updDone) {
 692         double cpar[2]={ bestCl->GetY(), bestCl->GetZ()};
 693         double ccov[3]={ bestCl->GetSigmaY2(), 0., bestCl->GetSigmaZ2()}; // RS: add syst errors
 694         if (!bestTr.Update(cpar,ccov)) return kFALSE;
 695         updDone = kTRUE;
 696       }
 697       track.paramOut = bestTr;
 698       track.clID[lrIDA] = iclBest;
 699       track.ncl++;
 700       track.chi2 += chi2Best;
 701     }
 702   }
 703   //
 704   if (!updDone) {
 705     if (++track.nmiss > fMaxMissedLayers)  return kFALSE;
 706     track.paramOut = trCopy;
 707     track.chi2 += fMissChi2Penalty;
 708   }
 709   //
 710 #ifdef _CONTROLH_
 711   int ndf = 2*track.ncl-5;
 712   if (ndf>0) {
 713     fHChi2vsNC->Fill(track.ncl,track.chi2);
 714     if (lrID==kNLrActive-1) fHChi2NDFvsPT->Fill(track.paramOut.Pt(),track.chi2/ndf);
 715   }
 716 #endif
 717   if (track.chi2 > GetChi2TotCut(track.ncl+1)) return kFALSE;
 718   //
 719   return track.paramOut.CorrectForMeanMaterial(fgkX2X0ITS[lrID],-fgkRhoLITS[lrID],fgkDefMass,kFALSE);
 720   //
 721 }
 722
 723 //______________________________________________
 724 void AliITSSAPTracker::CookLabel(AliITSSAPTracker::ITStrack_t& track)
 725 {
 726   // cook mc label for the track
 727   track.label = fgkDummyLabel;
 728   if (!track.ncl) return;
 729   const int kMaxLbPerCl = 3;
 730   int lbID[kNLrActive*6],lbStat[kNLrActive*6];
 731   Int_t nLab=0;
 732   for (int i=kNLrActive;i--;) {
 733     int clid = track.clID[i];
 734     if (clid<0) continue;
 735     AliITSRecPoint* cl = fLayers[i]->GetClusterSorted(clid);
 736     for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
 737       int trLb = cl->GetLabel(imc);
 738       if (trLb<0) break;
 739       // search this mc track in already accounted ones
 740       int iLab;
 741       for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
 742       if (iLab<nLab) lbStat[iLab]++;
 743       else {
 744         lbID[nLab] = trLb;
 745         lbStat[nLab++] = 1;
 746       }
 747     } // loop over given cluster's labels
 748   } // loop over all clusters
 749   //
 750   if (nLab) {
 751     int maxLab=0;
 752     for (int ilb=nLab;ilb--;) if (lbStat[maxLab]<lbStat[ilb]) maxLab=ilb;
 753     track.label = lbStat[maxLab]==track.ncl ? lbID[maxLab] : -lbID[maxLab];
 754   }
 755   //
 756 }
 757
 758 //______________________________________________
 759 void AliITSSAPTracker::CookLabel(AliITSSAPTracker::SPDtracklet_t& tracklet)
 760 {
 761   // cook mc label for the tracklet
 762   tracklet.label = fgkDummyLabel;
 763   const int kMaxLbPerCl = 3;
 764   int lbID[kNLrActive*6],lbStat[kNLrActive*6];
 765   Int_t nLab=0;
 766   for (int i=2;i--;) {
 767     int clid = i ? tracklet.id2 : tracklet.id1;
 768     AliITSRecPoint* cl = fLayers[i]->GetClusterSorted(clid);
 769     for (int imc=0;imc<kMaxLbPerCl;imc++) { // labels within single cluster
 770       int trLb = cl->GetLabel(imc);
 771       if (trLb<0) break;
 772       // search this mc track in already accounted ones
 773       int iLab;
 774       for (iLab=0;iLab<nLab;iLab++) if (lbID[iLab]==trLb) break;
 775       if (iLab<nLab) lbStat[iLab]++;
 776       else {
 777         lbID[nLab] = trLb;
 778         lbStat[nLab++] = 1;
 779       }
 780     } // loop over given cluster's labels
 781   } // loop over all clusters
 782   //
 783   if (nLab) {
 784     int maxLab=0;
 785     for (int ilb=nLab;ilb--;) if (lbStat[maxLab]<lbStat[ilb]) maxLab=ilb;
 786     tracklet.label = lbStat[maxLab]==2 ? lbID[maxLab] : -lbID[maxLab];
 787   }
 788   //
 789 }
 790
 791 //______________________________________________
 792 Double_t AliITSSAPTracker::GetXatLabRLin(AliExternalTrackParam& track, double r)
 793 {
 794   // X of track circle intersection in current tracking frame, neglecting the curvature
 795   // Solution of equation (x+d)^2+(y+b*d)^2 - r^2, where x,y are current coordinates of
 796   // track and d=X-x0. b = tg(phi)
 797   //double sn=tr.GetSnp();
 798   double sn=track.GetSnp();
 799   if (TMath::Abs(sn)>kAlmost1) return -999;
 800   double x=track.GetX(), y=track.GetY();
 801   double cs2=(1.-sn)*(1.+sn), tg=sn/TMath::Sqrt(cs2);
 802   double t0=x+tg*y, t1=x*x+y*y-r*r, det=t0*t0-t1/cs2;
 803   if (det<0) return -999; // does not touch circle
 804   det = TMath::Sqrt(det);
 805   return x+(det-t0)*cs2;
 806   //
 807 }
 808
 809 //______________________________________________
 810 Int_t AliITSSAPTracker::GetTrackletMCTruth(AliITSSAPTracker::SPDtracklet_t& trlet) const
 811 {
 812   int status = 0;
 813   AliITSSAPLayer::ClsInfo_t *cli1=fLayers[kALrSPD1]->GetClusterInfo(trlet.id1);
 814   AliITSSAPLayer::ClsInfo_t *cli2=fLayers[kALrSPD2]->GetClusterInfo(trlet.id2);
 815   AliITSRecPoint *cl1=fLayers[kALrSPD1]->GetClusterUnSorted(cli1->index);
 816   AliITSRecPoint *cl2=fLayers[kALrSPD2]->GetClusterUnSorted(cli2->index);
 817   //
 818   int lab = -1;
 819   //
 820   for (int i=0;i<3;i++) {
 821     int lb1 = cl1->GetLabel(i);
 822     if (lb1<0) continue;
 823     for (int j=0;j<3;j++) {
 824       int lb2 = cl2->GetLabel(i);
 825       if (lb2<0) break;
 826       if (lb1==lb2) {lab = lb1; break;}
 827     }
 828     if (lab>=0) break;
 829   }
 830   if (lab<0) return 0;
 831   //
 832   for (int ila=kALrSDD1;ila<kNLrActive;ila++) {
 833     for (int icl=fNClusters[ila];icl--;) {
 834       AliITSRecPoint *cl=fLayers[ila]->GetClusterUnSorted(icl);
 835       for (int i=0;i<3;i++) {
 836         if (cl->GetLabel(i)<0) break;
 837         if (cl->GetLabel(i)==lab) {status |= 0x1<<ila; break;}
 838       }
 839       if (status & (0x1<<ila)) break;
 840     }
 841   }
 842   return status;
 843 }
 844
 845 //______________________________________________
 846 Bool_t AliITSSAPTracker::RefitInward(int itr)
 847 {
 848   // refit track inward with material correction
 849   ITStrack_t &track = fTracks[itr];
 850   AliExternalTrackParam &trout = track.paramOut;
 851   if (trout.TestBit(kInvalidBit)) return kFALSE;
 852   AliExternalTrackParam &trin  = track.paramInw;
 853   trin = trout;
 854   int ilA = kNLrActive;
 855   for (;ilA--;) {                    // find outermost layer with cluster
 856     if (track.clID[ilA]<0) continue;
 857     break;
 858   }
 859   int ilStart = fgkActiveLrITS[ilA]; // corresponding total lr id
 860   AliITSSAPLayer* lrA = fLayers[ilA];
 861   AliITSRecPoint *cl=lrA->GetClusterSorted(track.clID[ilA]);
 862   AliITSSAPLayer::ITSDetInfo_t& detInfo = lrA->GetDetInfo(cl->GetVolumeId()-lrA->GetVIDOffset());
 863   if (!trin.RotateParamOnly(detInfo.phiTF)) return kFALSE;
 864   if (!trin.PropagateParamOnlyTo(detInfo.xTF+cl->GetX(), fBz)) return kFALSE;
 865   // init with outer cluster y,z and slopes, q/pt of outward track
 866   double par[5] = {cl->GetY(), cl->GetZ(), trin.GetSnp(), trin.GetTgl(), trin.GetSigned1Pt()};
 867   double cov[15] = {cl->GetSigmaY2() + GetClSystYErr2(kALrSPD1),
 868                    0., cl->GetSigmaZ2() + GetClSystZErr2(kALrSPD1),
 869                    0,0,1,
 870                    0,0,0,1,
 871                    0,0,0,0,100*100};
 872   trin.Set(double(detInfo.xTF+cl->GetX()),double(detInfo.phiTF), par, cov);
 873   // !!! no material correction is needed: errors are not defined yer
 874   //
 875   for (int ilr=ilStart;ilr--;) {
 876     //
 877     if ( (ilA=fgkLr2Active[ilr])<0 || track.clID[ilA]<0) { // either passive layer or no cluster
 878       if (CrossPassiveLayer(trin,ilr)) continue;
 879       else return kFALSE;
 880     }
 881     // there is a cluster, need to update
 882     lrA = fLayers[ilA];
 883     cl = lrA->GetClusterSorted(track.clID[ilA]);
 884     AliITSSAPLayer::ITSDetInfo_t& detInfo1 = lrA->GetDetInfo(cl->GetVolumeId()-lrA->GetVIDOffset());
 885     if (!trin.Propagate(detInfo1.phiTF, detInfo1.xTF+cl->GetX(), fBz)) return kFALSE;
 886     double cpar[2]={ cl->GetY(), cl->GetZ()};
 887     double ccov[3]={ cl->GetSigmaY2() + GetClSystYErr2(ilA) , 0., cl->GetSigmaZ2() + GetClSystZErr2(ilA)};
 888     if (!trin.Update(cpar,ccov)) return kFALSE;
 889     //
 890     // correct for layer materials
 891     if (!trin.CorrectForMeanMaterial(fgkX2X0ITS[ilr],fgkRhoLITS[ilr],fgkDefMass,kFALSE)) return kFALSE;
 892     //
 893   }
 894   //
 895   // now go to PCA to vertex
 896   //double dca[2],dcaCov[3];
 897   if (!trin.PropagateToDCA(fSPDVertex,fBz,fgkRLayITS[kLrBeamPime])) return kFALSE; //,dca,dcaCov);
 898   //
 899   trin.ResetBit(kInvalidBit); // flag that inward fit succeeded
 900   return kTRUE;
 901   //
 902 }
 903
 904 //______________________________________________
 905 void AliITSSAPTracker::RefitInward()
 906 {
 907   // refit tracks inward with material correction
 908   for (int itr=fNTracks;itr--;) {
 909     if (!RefitInward(itr)) {
 910 #ifdef _DEBUG_
 911       printf("RefitInward failed for track %d\n",itr);
 912       PrintTrack(fTracks[itr]);
 913 #endif
 914     }
 915   }
 916   //
 917 }
 918
 919
 920 //______________________________________________
 921 Bool_t AliITSSAPTracker::FitTrackVertex()
 922 {
 923   // Fit the vertexTracks. The inner tracks must be already propagated to the SPD vertex.
 924   // In this case straight line extrapolation can be used
 925   //
 926   fTrackVertex.SetNContributors(0); // invalidate
 927   //
 928   if (fNTracks<3) return kFALSE;
 929   double cxx=0,cxy=0,cxz=0,cx0=0,cyy=0,cyz=0,cy0=0,czz=0,cz0=0;
 930   //
 931   int ntAcc = 0;
 932   for (int itr=fNTracks;itr--;) {
 933     //
 934     AliExternalTrackParam& trc = fTracks[itr].paramInw;
 935     if (trc.TestBit(kInvalidBit)) continue; // the track is invalidated, skip
 936     double *param = (double*)trc.GetParameter();
 937     double *covar = (double*)trc.GetCovariance();
 938     //
 939     double  x0=trc.GetX();
 940     double &y0=param[0];
 941     double &z0=param[1];
 942     double sn=param[2];
 943     double cs2=(1.-sn)*(1.+sn);
 944     if (cs2<kAlmost0) continue;
 945     double cs=TMath::Sqrt(cs2), tgp=sn/cs, tgl=trc.GetTgl()/cs;
 946     // assume straight track equation Y=y0+tgp*X, Z=z0+tgl*X in tracking frame
 947     //
 948     double alp = trc.GetAlpha();
 949     sn = TMath::Sin(alp); // parameters for rotation of vertex to
 950     cs = TMath::Cos(alp); // tracking frame
 951     //
 952     double &syy=covar[0], &syz=covar[1], &szz=covar[2];
 953     double detI = syy*szz - syz*syz;
 954     if (TMath::Abs(detI)<kAlmost0) return kFALSE;
 955     detI = 1./detI;
 956     double syyI = szz*detI;
 957     double szzI = syy*detI;
 958     double syzI =-syz*detI;
 959     //
 960     double tmpSP = sn*tgp;
 961     double tmpCP = cs*tgp;
 962     double tmpSC = sn+tmpCP;
 963     double tmpCS =-cs+tmpSP;
 964     double tmpCL = cs*tgl;
 965     double tmpSL = sn*tgl;
 966     double tmpYXP = y0-tgp*x0;
 967     double tmpZXL = z0-tgl*x0;
 968     //
 969     double tmpCLzz = tmpCL*szzI;
 970     double tmpSLzz = tmpSL*szzI;
 971     double tmpSCyz = tmpSC*syzI;
 972     double tmpCSyz = tmpCS*syzI;
 973     double tmpCSyy = tmpCS*syyI;
 974     double tmpSCyy = tmpSC*syyI;
 975     double tmpSLyz = tmpSL*syzI;
 976     double tmpCLyz = tmpCL*syzI;
 977     //
 978     cxx += tmpCL*(tmpCLzz+tmpSCyz+tmpSCyz)+tmpSC*tmpSCyy;          // dchi^2/dx/dx
 979     cxy += tmpCL*(tmpSLzz+tmpCSyz)+tmpSL*tmpSCyz+tmpSC*tmpCSyy;    // dchi^2/dx/dy
 980     cxz += -sn*syzI-tmpCLzz-tmpCP*syzI;                            // dchi^2/dx/dz
 981     cx0 += -(tmpCLyz+tmpSCyy)*tmpYXP-(tmpCLzz+tmpSCyz)*tmpZXL;     // RHS
 982     //
 983     //double cyx
 984     cyy += tmpSL*(tmpSLzz+tmpCSyz+tmpCSyz)+tmpCS*tmpCSyy;          // dchi^2/dy/dy
 985     cyz += -(tmpCSyz+tmpSLzz);                                     // dchi^2/dy/dz
 986     cy0 += -tmpYXP*(tmpCSyy+tmpSLyz)-tmpZXL*(tmpCSyz+tmpSLzz);     // RHS
 987     //
 988     //double czx
 989     //double czy
 990     czz += szzI;                                                    // dchi^2/dz/dz
 991     cz0 += tmpZXL*szzI+tmpYXP*syzI;                                 // RHS
 992     //
 993     ntAcc++;
 994   }
 995   //
 996   double vec[3] = {cx0,cy0,cz0};
 997   AliSymMatrix mat(3);
 998   mat(0,0) = cxx;
 999   mat(0,1) = cxy;
1000   mat(0,2) = cxz;
1001   mat(1,1) = cyy;
1002   mat(1,2) = cyz;
1003   mat(2,2) = czz;
1004
1005   //-------------------------TMP>>>
1006   AliRunLoader* rl = AliRunLoader::Instance();
1007   AliHeader* hd = 0;
1008   AliGenEventHeader* hdmc=0;
1009   TArrayF vtxMC(3);
1010   if (rl && (hd=rl->GetHeader()) && (hdmc=hd->GenEventHeader())) {
1011     hdmc->PrimaryVertex(vtxMC);
1012   }
1013   //-------------------------TMP<<<
1014 #ifdef _DEBUG_
1015   printf("MatBefore: \n"); mat.Print("d");
1016 #endif
1017   if (mat.SolveChol(vec,kTRUE)) {
1018 #ifdef _DEBUG_
1019     printf("MatAfter : \n"); mat.Print("d");
1020 #endif
1021     //
1022     double vtCov[6] = {mat(0,0),mat(0,1),mat(1,1),mat(0,2),mat(1,2),mat(2,2)};
1023     fTrackVertex.SetXYZ(vec);
1024     fTrackVertex.SetCovarianceMatrix(vtCov);
1025     fTrackVertex.SetNContributors(ntAcc);
1026     //
1027     // calculate explicitly chi2
1028     double chiTRC = 0;
1029     double chiSPD = 0;
1030     //
1031     for (int itr=fNTracks;itr--;) {
1032       AliExternalTrackParam& trc = fTracks[itr].paramInw;
1033       if (trc.TestBit(kInvalidBit)) continue; // the track is invalidated, skip
1034       AliExternalTrackParam trT(trc);
1035       AliExternalTrackParam trS(trc);
1036       double dz[2],covdum[3],*covt;
1037       trT.PropagateToDCA(&fTrackVertex,fBz,10,dz,covdum);
1038       covt = (double*)trT.GetCovariance();
1039       double detI = covt[0]*covt[2] - covt[1]*covt[1];
1040       detI = 1./detI;
1041       double syyI = covt[2]*detI;
1042       double szzI = covt[0]*detI;
1043       double syzI =-covt[1]*detI;
1044       chiTRC += dz[0]*dz[0]*syyI + dz[1]*dz[1]*szzI + 2*dz[0]*dz[1]*syzI;
1045       //
1046       trS.PropagateToDCA(fSPDVertex,fBz,10,dz,covdum);
1047       covt = (double*)trT.GetCovariance();
1048       detI = covt[0]*covt[2] - covt[1]*covt[1];
1049       detI = 1./detI;
1050       syyI = covt[2]*detI;
1051       szzI = covt[0]*detI;
1052       syzI =-covt[1]*detI;
1053       chiSPD += dz[0]*dz[0]*syyI + dz[1]*dz[1]*szzI + 2*dz[0]*dz[1]*syzI;
1054     }
1055 #ifdef _DEBUG_
1056     printf("VTFIT %f %f %f %d %8.2f %8.2f   %.4f %.4f %.4f   %.4f %.4f %.4f\n",
1057            vtxMC[0],vtxMC[1],vtxMC[2],
1058            ntAcc,chiTRC,chiSPD,
1059            fTrackVertex.GetX(),fTrackVertex.GetY(),fTrackVertex.GetZ(),
1060            fSPDVertex->GetX(),fSPDVertex->GetY(),fSPDVertex->GetZ());
1061 #endif
1062     //
1063     return kTRUE;
1064   }
1065   //
1066   return kFALSE;
1067 }
1068
1069 #ifdef _CONTROLH_
1070 //______________________________________________
1071 void AliITSSAPTracker::FillRecoStat()
1072 {
1073   // fill data for preformance study
1074   //
1075   AliStack* stack = 0;
1076   AliRunLoader* rl = AliRunLoader::Instance();
1077   if (!rl || !(stack=rl->Stack())) return;
1078   //
1079   TBits patternMC;
1080   enum {kIsPrim=kNLrActive,kValidTracklet,kValidTrack,kRecDone,kBitPerTrack};
1081   int nTrkMC = stack->GetNtrack();
1082   patternMC.SetBitNumber(nTrkMC*kBitPerTrack,0);
1083   //
1084   // fill MC track patterns
1085   for (int ilr=kNLrActive;ilr--;) {
1086     AliITSSAPLayer *lr = fLayers[ilr];
1087     int ncl = lr->GetNClusters();
1088     for (int icl=ncl;icl--;) {
1089       AliITSRecPoint* cl = lr->GetClusterUnSorted(icl);
1090       for (int j=0;j<3;j++) {
1091         int lb = cl->GetLabel(j);
1092         if (lb<0 || lb>=nTrkMC) break;
1093         patternMC.SetBitNumber(lb*kBitPerTrack+ilr,kTRUE);
1094       }
1095     }
1096   }
1097   // set reconstructability
1098   for (int itr=nTrkMC;itr--;) {
1099     int bitoffs = itr*kBitPerTrack;
1100     Bool_t isPrim = stack->IsPhysicalPrimary(itr);
1101     patternMC.SetBitNumber(bitoffs+kIsPrim,isPrim);
1102     if (patternMC.TestBitNumber(bitoffs+kALrSPD1) && patternMC.TestBitNumber(bitoffs+kALrSPD2)) {
1103       patternMC.SetBitNumber(bitoffs+kValidTracklet,kTRUE);
1104       //
1105       TParticle* mctr = stack->Particle(itr);
1106       fHTrackletMC->Fill(mctr->Pt(),isPrim);
1107       // check outer layers reconstructability
1108       int nmiss = 0;
1109       for (int il=kALrSDD1;il<=kALrSSD2;il++) if (!fSkipLayer[il] && !patternMC.TestBitNumber(bitoffs+il)) nmiss++;
1110       if (nmiss<=fMaxMissedLayers) {
1111         patternMC.SetBitNumber(bitoffs+kValidTrack);
1112         fHTrackMC->Fill(mctr->Pt(),isPrim);
1113       }
1114     }
1115   }
1116   //
1117   int nTrk = GetNTracklets();
1118   if (!nTrk) return;
1119   for (int itr=0;itr<nTrk;itr++) {
1120     SPDtracklet_t& trlet = fTracklets[itr];
1121     //    PrintTracklet(itr);
1122     //
1123     int lbl = trlet.label;
1124     if (lbl==fgkDummyLabel) continue;
1125     int lblA = TMath::Abs(lbl);
1126     int bitoffs = lblA*kBitPerTrack;
1127     Bool_t isPrim = patternMC.TestBitNumber(bitoffs+kIsPrim);
1128     TParticle* mctr = stack->Particle(lblA);
1129     double pt = mctr->Pt();
1130     fHTrackletAll->Fill(pt,isPrim);
1131     if (lbl<0) fHTrackletFake->Fill(pt,isPrim);
1132   }
1133   //
1134   nTrk = GetNTracks();
1135   for (int itr=0;itr<nTrk;itr++) {
1136     ITStrack_t &track = fTracks[itr];
1137     CookLabel(track);
1138     //
1139     int lbl = track.label;
1140     if (lbl==fgkDummyLabel) continue;
1141     int lblA = TMath::Abs(lbl);
1142     int bitoffs = lblA*kBitPerTrack;
1143     Bool_t isPrim = patternMC.TestBitNumber(bitoffs+kIsPrim);
1144     TParticle* mctr = stack->Particle(lblA);
1145     double pt = mctr->Pt();
1146     Bool_t clone = patternMC.TestBitNumber(bitoffs+kRecDone); // was the track already reconstructed?
1147     float bn = isPrim ? (clone ? 2:1):(clone ? -1:0);  // fill clones in over/underflow
1148     fHTrackAll->Fill(pt,bn);
1149     patternMC.SetBitNumber(bitoffs+kRecDone);
1150     if (lbl<0) fHTrackFake->Fill(pt,bn);
1151   }
1152   //
1153   AliHeader* hd = rl->GetHeader();
1154   AliGenEventHeader* hdmc;
1155   TArrayF vtxMC;
1156   if (hd && (hdmc=hd->GenEventHeader())) hdmc->PrimaryVertex(vtxMC);
1157   //
1158   nTrk = GetNTracklets();
1159   fHVtxMltRef->Fill(nTrk);
1160   if (fTrackVertex.GetStatus()==1) {
1161     if (hdmc) {
1162       double dx = vtxMC[0]-fTrackVertex.GetX();
1163       double dy = vtxMC[1]-fTrackVertex.GetY();
1164       double dz = vtxMC[2]-fTrackVertex.GetZ();
1165       fHVtxDiffXY->Fill(dx,dy);
1166       fHVtxDiffZ->Fill(dz);
1167       fHVtxDiffXMlt->Fill(nTrk, dx);
1168       fHVtxDiffYMlt->Fill(nTrk, dy);
1169       fHVtxDiffZMlt->Fill(nTrk, dz);
1170       //
1171       double sig[3];
1172       fTrackVertex.GetSigmaXYZ(sig);
1173       if (sig[0]>0) fHVtxPullXMlt->Fill(nTrk, dx/sig[0]);
1174       if (sig[1]>0) fHVtxPullYMlt->Fill(nTrk, dy/sig[1]);
1175       if (sig[2]>0) fHVtxPullZMlt->Fill(nTrk, dz/sig[2]);
1176     }
1177     fHVtxOKMlt->Fill(nTrk);
1178   }
1179   //
1180   if (fSPDVertex->GetStatus()==1 && hdmc) {
1181     double dx = vtxMC[0]-fSPDVertex->GetX();
1182     double dy = vtxMC[1]-fSPDVertex->GetY();
1183     double dz = vtxMC[2]-fSPDVertex->GetZ();
1184     fHVtxMCSPDDiffXY->Fill(dx,dy);
1185     fHVtxMCSPDDiffZ->Fill(dz);
1186     fHVtxMCSPDDiffXMlt->Fill(nTrk, dx);
1187     fHVtxMCSPDDiffYMlt->Fill(nTrk, dy);
1188     fHVtxMCSPDDiffZMlt->Fill(nTrk, dz);
1189     //
1190     double sig[3];
1191     fSPDVertex->GetSigmaXYZ(sig);
1192     if (sig[0]>0) fHVtxMCSPDPullXMlt->Fill(nTrk, dx/sig[0]);
1193     if (sig[1]>0) fHVtxMCSPDPullYMlt->Fill(nTrk, dy/sig[1]);
1194     if (sig[2]>0) fHVtxMCSPDPullZMlt->Fill(nTrk, dz/sig[2]);
1195     //
1196   }
1197   //
1198 }
1199
1200 //______________________________________________
1201 void AliITSSAPTracker::BookHistos()
1202 {
1203   // book control histos
1204   const int kNBinMlt=20, kNBPt=15, kNBDiffVtx=50, kNResBins=250,kNPullBins=50,kNChiClBins=50,kNBPullVtx=50;
1205   const double kMinMlt=1,kMaxMlt=5000,kMinPt=0.01,kMaxPt=3, kMaxDiffVtx=0.05, kMaxResidYZ=2.5,kMaxPullYZ=10,kChiClMax=100,kMaxPullVtx=10;
1206   //
1207   double* axLogPt  = DefLogAx(kMinPt,kMaxPt,kNBPt);
1208   double* axLogMlt = DefLogAx(kMinMlt,kMaxMlt,kNBinMlt);
1209
1210   for (int ilr=0;ilr<kNLrActive;ilr++) {
1211     //
1212     // ----------------- These are histos to be filled during tracking
1213     // PropagateBack and RefitInward will be stored among the histos of 1st pass
1214     //
1215     int ilrS = ilr*10;
1216     TString ttl = Form("residY%d",ilr);
1217     TH2F* hdy = new TH2F(ttl.Data(),ttl.Data(),kNBPt,axLogPt,kNResBins,-kMaxResidYZ,kMaxResidYZ);
1218     fArrHisto.AddAtAndExpand(hdy,ilrS+kHResidY);
1219     hdy->SetDirectory(0);
1220     //
1221     ttl = Form("residYPull%d",ilr);
1222     TH2F* hdyp = new TH2F(ttl.Data(),ttl.Data(),kNBPt,axLogPt,kNPullBins,-kMaxPullYZ,kMaxPullYZ);
1223     fArrHisto.AddAtAndExpand(hdyp,ilrS+kHPullY);
1224     hdyp->SetDirectory(0);
1225     //
1226     ttl = Form("residZ%d",ilr);
1227     TH2F* hdz = new TH2F(ttl.Data(),ttl.Data(),kNBPt,axLogPt,kNResBins,-kMaxResidYZ,kMaxResidYZ);
1228     fArrHisto.AddAtAndExpand(hdz,ilrS+kHResidZ);
1229     hdz->SetDirectory(0);
1230     //
1231     ttl = Form("residZPull%d",ilr);
1232     TH2F* hdzp = new TH2F(ttl.Data(),ttl.Data(),kNBPt,axLogPt,kNPullBins,-kMaxPullYZ,kMaxPullYZ);
1233     hdzp->SetDirectory(0);
1234     fArrHisto.AddAtAndExpand(hdzp,ilrS+kHPullZ);
1235     //
1236     ttl = Form("chi2Cl%d",ilr);
1237     TH2F* hchi = new TH2F(ttl.Data(),ttl.Data(),kNBPt,axLogPt, kNChiClBins,0.,kChiClMax);
1238     hchi->SetDirectory(0);
1239     fArrHisto.AddAtAndExpand(hchi,ilrS+kHChi2Cl);
1240   } // loop over layers
1241   //
1242   fHChi2NDFvsPT = new TH2F("chi2ndfPT","chi2/ndf total vs pt",kNBPt,axLogPt, kNChiClBins,0.,kChiClMax);
1243   fArrHisto.AddLast(fHChi2NDFvsPT);
1244   fHChi2NDFvsPT->SetDirectory(0);
1245   //
1246   fHChi2vsNC = new TH2F("chi2NC","chi2 total vs NCl",kNLrActive-2,2.5,kNLrActive+0.5, kNChiClBins,0.,kChiClMax);
1247   fArrHisto.AddLast(fHChi2vsNC);
1248   fHChi2vsNC->SetDirectory(0);
1249
1250   // SPDvertex vs MC
1251   fHVtxMCSPDDiffXY = new TH2F("vtxMCSPDDiffXY","vtxMC-vtxSPD XY",kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx,
1252                             kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1253   fArrHisto.AddLast(fHVtxMCSPDDiffXY);
1254   fHVtxMCSPDDiffXY->SetDirectory(0);
1255   //
1256   fHVtxMCSPDDiffZ = new TH1F("vtxMCSPDDiffZ","vtxMC-vtxSPD Z",kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1257   fArrHisto.AddLast(fHVtxMCSPDDiffZ);
1258   fHVtxMCSPDDiffZ->SetDirectory(0);
1259   //
1260   fHVtxMCSPDDiffXMlt = new TH2F("VtxMCSPDDiffXMlt","vX_{MC}-vX_{SPD} vs mlt",kNBinMlt,axLogMlt,
1261                            kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1262   fArrHisto.AddLast(fHVtxMCSPDDiffXMlt);
1263   fHVtxMCSPDDiffXMlt->SetDirectory(0);
1264   //
1265   fHVtxMCSPDDiffYMlt = new TH2F("VtxMCSPDDiffYMlt","vY_{MC}-vY_{rec} vs mlt",kNBinMlt,axLogMlt,
1266                            kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1267   fArrHisto.AddLast(fHVtxMCSPDDiffYMlt);
1268   fHVtxMCSPDDiffYMlt->SetDirectory(0);
1269   //
1270   fHVtxMCSPDDiffZMlt = new TH2F("VtxMCSPDDiffZMlt","vZ_{MC}-vZ_{rec} vs mlt",kNBinMlt,axLogMlt,
1271                            kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1272   fArrHisto.AddLast(fHVtxMCSPDDiffZMlt);
1273   fHVtxMCSPDDiffZMlt->SetDirectory(0);
1274   //
1275   //
1276   fHVtxMCSPDPullXMlt = new TH2F("VtxMCSPDPullXMlt","Pull vX_{MC}-vX_{SPD} vs mlt",kNBinMlt,axLogMlt,
1277                            kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1278   fArrHisto.AddLast(fHVtxMCSPDPullXMlt);
1279   fHVtxMCSPDPullXMlt->SetDirectory(0);
1280   //
1281   fHVtxMCSPDPullYMlt = new TH2F("VtxMCSPDPullYMlt","Pull vY_{MC}-vY_{rec} vs mlt",kNBinMlt,axLogMlt,
1282                            kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1283   fArrHisto.AddLast(fHVtxMCSPDPullYMlt);
1284   fHVtxMCSPDPullYMlt->SetDirectory(0);
1285   //
1286   fHVtxMCSPDPullZMlt = new TH2F("VtxMCSPDPullZMlt","Pull vZ_{MC}-vZ_{rec} vs mlt",kNBinMlt,axLogMlt,
1287                            kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1288   fArrHisto.AddLast(fHVtxMCSPDPullZMlt);
1289   fHVtxMCSPDPullZMlt->SetDirectory(0);
1290    //
1291   fHTrackletMC = new TH2F("MCRefTracklet","MCRef Tracklet",kNBPt,axLogPt, 2, -0.5, 1.5);
1292   fHTrackletMC->SetXTitle("p_{T}");
1293   fHTrackletMC->GetYaxis()->SetBinLabel(1,"sec");
1294   fHTrackletMC->GetYaxis()->SetBinLabel(2,"prim");
1295   fArrHisto.AddLast(fHTrackletMC);
1296   fHTrackletMC->SetDirectory(0);
1297   //
1298   fHTrackletAll = new TH2F("TrackletAll","Tracklet All rec",kNBPt,axLogPt, 2, -0.5, 1.5);
1299   fHTrackletAll->SetXTitle("p_{T}");
1300   fHTrackletAll->GetYaxis()->SetBinLabel(1,"sec");
1301   fHTrackletAll->GetYaxis()->SetBinLabel(2,"prim");
1302   fArrHisto.AddLast(fHTrackletAll);
1303   fHTrackletAll->SetDirectory(0);
1304   //
1305   fHTrackletFake = new TH2F("TrackletFake","Tracklet Fake rec",kNBPt,axLogPt, 2, -0.5, 1.5);
1306   fHTrackletFake->SetXTitle("p_{T}");
1307   fHTrackletFake->GetYaxis()->SetBinLabel(1,"sec");
1308   fHTrackletFake->GetYaxis()->SetBinLabel(2,"prim");
1309   fArrHisto.AddLast(fHTrackletFake);
1310   fHTrackletFake->SetDirectory(0);
1311   //
1312   fHTrackMC = new TH2F("MCRefTrack","MCRef Track",kNBPt,axLogPt, 2, -0.5, 1.5);
1313   fHTrackMC->SetXTitle("p_{T}");
1314   fHTrackMC->GetYaxis()->SetBinLabel(1,"sec");
1315   fHTrackMC->GetYaxis()->SetBinLabel(2,"prim");
1316   fArrHisto.AddLast(fHTrackMC);
1317   fHTrackMC->SetDirectory(0);
1318   //
1319   fHTrackAll = new TH2F("TrackAll","Track All rec",kNBPt,axLogPt, 2, -0.5, 1.5);
1320   fHTrackAll->SetXTitle("p_{T}");
1321   fHTrackAll->GetYaxis()->SetBinLabel(1,"sec");
1322   fHTrackAll->GetYaxis()->SetBinLabel(2,"prim");
1323   fArrHisto.AddLast(fHTrackAll);
1324   fHTrackAll->SetDirectory(0);
1325   //
1326   fHTrackFake = new TH2F("TrackFake","Track Fake rec",kNBPt,axLogPt, 2, -0.5, 1.5);
1327   fHTrackFake->SetXTitle("p_{T}");
1328   fHTrackFake->GetYaxis()->SetBinLabel(1,"sec");
1329   fHTrackFake->GetYaxis()->SetBinLabel(2,"prim");
1330   fArrHisto.AddLast(fHTrackFake);
1331   fHTrackFake->SetDirectory(0);
1332   //
1333   fHVtxMltRef = new TH1F("vtxRef","vtxRef",kNBinMlt,axLogMlt);
1334   fHVtxMltRef->SetXTitle("sqrt(Ntracklets)");
1335   fArrHisto.AddLast(fHVtxMltRef);
1336   fHVtxMltRef->SetDirectory(0);
1337   //
1338   fHVtxOKMlt = new TH1F("vtxOK","vtxOK",kNBinMlt,axLogMlt);
1339   fHVtxOKMlt->SetXTitle("sqrt(Ntracklets)");
1340   fArrHisto.AddLast(fHVtxOKMlt);
1341   fHVtxOKMlt->SetDirectory(0);
1342   //
1343   fHVtxDiffXY = new TH2F("vtxMCDiffXY","vtxMC-vtxRec XY",kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx,
1344                          kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1345   fArrHisto.AddLast(fHVtxDiffXY);
1346   fHVtxDiffXY->SetDirectory(0);
1347   //
1348   fHVtxDiffZ = new TH1F("vtxMCDiffZ","vtxMC-vtxRec Z",kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1349   fArrHisto.AddLast(fHVtxDiffZ);
1350   fHVtxDiffZ->SetDirectory(0);
1351   //
1352   fHVtxDiffXMlt = new TH2F("VtxDiffXMlt","vX_{MC}-vX_{rec} vs mlt",kNBinMlt,axLogMlt,
1353                            kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1354   fArrHisto.AddLast(fHVtxDiffXMlt);
1355   fHVtxDiffXMlt->SetDirectory(0);
1356   //
1357   fHVtxDiffYMlt = new TH2F("VtxDiffYMlt","vY_{MC}-vY_{rec} vs mlt",kNBinMlt,axLogMlt,
1358                            kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1359   fArrHisto.AddLast(fHVtxDiffYMlt);
1360   fHVtxDiffYMlt->SetDirectory(0);
1361   //
1362   fHVtxDiffZMlt = new TH2F("VtxDiffZMlt","vZ_{MC}-vZ_{rec} vs mlt",kNBinMlt,axLogMlt,
1363                            kNBDiffVtx,-kMaxDiffVtx,kMaxDiffVtx);
1364   fArrHisto.AddLast(fHVtxDiffZMlt);
1365   fHVtxDiffZMlt->SetDirectory(0);
1366   //
1367   fHVtxPullXMlt = new TH2F("VtxPullXMlt","Pull vX_{MC}-vX_{rec} vs mlt",kNBinMlt,axLogMlt,
1368                            kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1369   fArrHisto.AddLast(fHVtxPullXMlt);
1370   fHVtxPullXMlt->SetDirectory(0);
1371   //
1372   fHVtxPullYMlt = new TH2F("VtxPullYMlt","Pull vY_{MC}-vY_{rec} vs mlt",kNBinMlt,axLogMlt,
1373                            kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1374   fArrHisto.AddLast(fHVtxPullYMlt);
1375   fHVtxPullYMlt->SetDirectory(0);
1376   //
1377   fHVtxPullZMlt = new TH2F("VtxPullZMlt","Pull vZ_{MC}-vZ_{rec} vs mlt",kNBinMlt,axLogMlt,
1378                            kNBPullVtx,-kMaxPullVtx,kMaxPullVtx);
1379   fArrHisto.AddLast(fHVtxPullZMlt);
1380   fHVtxPullZMlt->SetDirectory(0);
1381   //
1382 }
1383
1384 //______________________________________________
1385 void AliITSSAPTracker::SaveHistos(const char* outFName)
1386 {
1387   // save control histos
1388   TString fnms = outFName;
1389   if (fnms.IsNull()) fnms = "XXXITSTrackerControlH.root";
1390   TFile* fl = TFile::Open(fnms.Data(),"recreate");
1391   if (!fl) {
1392     printf("Failed to open output file %s\n",fnms.Data());
1393     return;
1394   }
1395   fArrHisto.Write();
1396   fl->Close();
1397   delete fl;
1398   printf("Stored control histos in %s\n",fnms.Data());
1399   //
1400 }
1401
1402 //______________________________________________
1403 void AliITSSAPTracker::FillTrackingControlHistos(int lrID,int lbl,const AliExternalTrackParam* track,
1404                                                  const double cpar[2],const double ccov[3],
1405                                                  const AliITSRecPoint* bestCl)
1406 {
1407   // fill control histos for tracking for correct matches
1408   Bool_t corr = kFALSE;
1409   for (int i=0;i<3;i++) if (bestCl->GetLabel(i)==lbl) {corr=kTRUE; break;}
1410   if (!corr) return;
1411   double pt = track->Pt();
1412   double dy = cpar[0]-track->GetY();
1413   double dz = cpar[1]-track->GetZ();
1414   double sgy = TMath::Sqrt(ccov[0]+track->GetSigmaY2());
1415   double sgz = TMath::Sqrt(ccov[2]+track->GetSigmaZ2());
1416   int lrIDS = lrID*10;
1417   ((TH2F*)fArrHisto[lrIDS+kHResidY])->Fill(pt,dy);
1418   ((TH2F*)fArrHisto[lrIDS+kHPullY])->Fill(pt,dy/sgy);
1419   ((TH2F*)fArrHisto[lrIDS+kHResidZ])->Fill(pt,dz);
1420   ((TH2F*)fArrHisto[lrIDS+kHPullZ])->Fill(pt,dz/sgz);
1421   ((TH2F*)fArrHisto[lrIDS+kHChi2Cl])->Fill(pt,track->GetPredictedChi2(cpar,ccov));
1422   //
1423 }
1424
1425 //______________________________________________
1426 Double_t* AliITSSAPTracker::DefLogAx(double xMn,double xMx, int nbin)
1427 {
1428   // get array for log axis
1429   if (xMn<=0 || xMx<=xMn || nbin<2) {
1430     printf("Wrong axis request: %f %f %d\n",xMn,xMx,nbin);
1431     return 0;
1432   }
1433   double dx = log(xMx/xMn)/nbin;
1434   double *xax = new Double_t[nbin+1];
1435   for (int i=0;i<=nbin;i++) xax[i]= xMn*exp(dx*i);
1436   return xax;
1437 }
1438
1439 #endif
1440 //
1441 #ifdef _TIMING_
1442 //______________________________________________
1443 void AliITSSAPTracker::PrintTiming()
1444 {
1445   // print timing info
1446   for (int i=0;i<kNSW;i++) {printf("%-10s:\t",fgkSWNames[i]); fSW[i].Print();}
1447 }
1448 #endif