Package for fast semi-analitical MC with realistic geometry
authorshahoian <ruben.shahoyan@cern.ch>
Thu, 11 Dec 2014 19:25:14 +0000 (20:25 +0100)
committershahoian <ruben.shahoyan@cern.ch>
Thu, 11 Dec 2014 19:25:14 +0000 (20:25 +0100)
ITS/UPGRADE/FT2/FT2.cxx [new file with mode: 0644]
ITS/UPGRADE/FT2/FT2.h [new file with mode: 0644]
ITS/UPGRADE/FT2/README [new file with mode: 0644]
ITS/UPGRADE/FT2/test.C [new file with mode: 0644]
ITS/UPGRADE/FT2/testF.C [new file with mode: 0644]

diff --git a/ITS/UPGRADE/FT2/FT2.cxx b/ITS/UPGRADE/FT2/FT2.cxx
new file mode 100644 (file)
index 0000000..aa6231e
--- /dev/null
@@ -0,0 +1,1243 @@
+#include "FT2.h"
+#include "AliCDBManager.h"
+#include "AliCDBEntry.h"
+#include "AliGRPManager.h"
+#include "AliLog.h"
+#include "AliITSUReconstructor.h"
+#include "AliITSURecoDet.h"
+#include <TGeoGlobalMagField.h>
+#include "AliMagF.h"
+#include "AliGeomManager.h"
+#include "AliITSUGeomTGeo.h"
+#include "AliTrackerBase.h"
+#include "AliITSUTrackerGlo.h"
+#include "AliITSURecoSens.h"
+#include "AliITSURecoLayer.h"
+#include "AliITSsegmentation.h"
+#include "AliESDVertex.h"
+#include <TParticle.h>
+#include <TDatabasePDG.h>
+#include <TRandom.h>
+#include <../TPC/Base/AliTPCcalibDB.h>
+#include <../TPC/Base/AliTPCParam.h>
+#include "AliPIDResponse.h"
+#include "AliDetectorPID.h"
+
+
+using namespace AliITSUAux;
+
+
+ClassImp(FT2)
+
+float FT2::fgMaxStepTGeo = 1.0;
+
+//________________________________________________
+FT2::FT2() :
+fITSRec(0)
+  ,fITS(0)
+  ,fUsePIDForTracking(0)
+  ,fIsTPC(kFALSE)
+  ,fPIDResponse(0)
+  ,fTPCSignalElectron(0)
+  ,fTPCSignalMuon(0)
+  ,fTPCSignalPion(0)
+  ,fTPCSignalKaon(0)
+  ,fTPCSignalProton(0)
+  ,fTPCSectorEdge(0)
+  ,fMaxSnpTPC(0.95)
+  ,fTPCLayers()
+  ,fTPCHitLr()
+  ,fNTPCHits(0)
+  ,fMinTPCHits(60)
+  ,fMinITSLrHit(4)
+  ,fProbe()
+  ,fKalmanOutward(0)
+  ,fUseKalmanOut(kTRUE)
+  //,fProbeMass(0.14)
+  ,fBz(0.5)
+  ,fSimMat(kFALSE)
+  ,fCurrITSLr(-1)
+  ,fNClTPC(0)
+  ,fNClITS(0)
+  ,fNClITSFakes(0)
+  ,fITSPatternFake(0)  
+  ,fITSPattern(0)
+  ,fChi2TPC(0)
+  ,fChi2ITS(0)
+  ,fSigYITS(3.14e-4) // 5e-4
+  ,fSigZITS(3.38e-4) // 5e-4
+  ,fNITSLrHit(0)
+  ,fdNdY(-1)
+{
+  //
+  AliInfo("Default Constructor");
+
+  memset(fITSSensHit,0,kMaxITSLr*2*sizeof(AliITSURecoSens*));
+}
+
+//________________________________________________
+FT2::~FT2()
+{
+  //destroy
+  AliInfo("Destroy");
+  delete[] fKalmanOutward;
+  delete fITSRec;
+}
+
+//________________________________________________
+void FT2::InitTPCSignalFiles(const char *TPCsignalElectrons,const char *TPCsignalMuons,const char *TPCsignalPions,const char *TPCsignalKaons,const char *TPCsignalProtons)
+{
+  AliInfo("Setting Files");
+       
+  fTPCSignalElectron = TFile::Open(TPCsignalElectrons);
+  if(fTPCSignalElectron->IsZombie()) AliFatal("Problem with opening TPC signal file for electrons - File not available!");
+       
+  fTPCSignalMuon = TFile::Open(TPCsignalMuons);
+  if(fTPCSignalMuon->IsZombie()) AliFatal("Problem with opening TPC signal file for muons - File not available!");
+       
+  fTPCSignalPion = TFile::Open(TPCsignalPions);
+  if(fTPCSignalPion->IsZombie()) AliFatal("Problem with opening TPC signal file for pions - File not available!");
+       
+  fTPCSignalKaon = TFile::Open(TPCsignalKaons);
+  if(fTPCSignalKaon->IsZombie()) AliFatal("Problem with opening TPC signal file for kaons - File not available!");
+       
+  fTPCSignalProton = TFile::Open(TPCsignalProtons);
+  if(fTPCSignalProton->IsZombie()) AliFatal("Problem with opening TPC signal file for protons - File not available!");
+}
+//________________________________________________
+void FT2::InitDetector(Bool_t addTPC, Float_t sigYTPC,Float_t sigZTPC,Float_t effTPC,Float_t scEdge)
+{
+  AliInfo("Initializing Detector");
+
+  // read full geometry and initialize qauzi-layers for TPC if requested
+  AliCDBManager* man = AliCDBManager::Instance();
+  if (!man->IsDefaultStorageSet() || man->GetRun()<0) AliFatal("CDB manager is not configured");
+  //
+  if ( !TGeoGlobalMagField::Instance()->GetField() ) {
+    AliWarning("Magnetic field is not initialized, loading from GRP");
+    AliGRPManager grpMan;
+    if(!grpMan.ReadGRPEntry()) AliFatal("Cannot get GRP entry");
+    if(!grpMan.SetMagField()) AliFatal("Problem with magnetic field setup");
+  }
+       
+  fBz = AliTrackerBase::GetBz();
+  //
+  AliGeomManager::LoadGeometry("geometry.root");
+  AliCDBEntry* ent = man->Get("ITS/Calib/RecoParam");
+  AliITSURecoParam* par = (AliITSURecoParam*)((TObjArray*)ent->GetObject())->At(1); // need just to initialize the interface
+  fITSRec = new AliITSUReconstructor();
+  fITSRec->SetRecoParam(par);
+  fITSRec->Init();
+  //
+  fITS = fITSRec->GetITSInterface();
+  if (addTPC) AddTPC(sigYTPC,sigZTPC,effTPC,scEdge);
+  int nLrITS = fITS->GetNLayersActive();
+  fKalmanOutward = new AliExternalTrackParam[nLrITS];
+  //
+  AliTPCcalibDB * calib = AliTPCcalibDB::Instance();
+  const AliMagF * field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
+  calib->SetExBField(field);
+  calib->SetRun(138871);
+       
+  //
+}
+//____________________________________________________
+void FT2::InitTPCPIDResponse()
+{
+  AliInfo("Initializing TPC PID Response");
+       
+  fPIDResponse = new AliPIDResponse(kTRUE);
+  fPIDResponse->GetTPCResponse().SetUseDatabase(kFALSE);
+  AliTPCParam* param = AliTPCcalibDB::Instance()->GetParameters();
+  if (param){
+    TVectorD *paramBB = param->GetBetheBlochParameters();
+    if (paramBB){
+      fPIDResponse->GetTPCResponse().SetBetheBlochParameters((*paramBB)(0),(*paramBB)(1),(*paramBB)(2),(*paramBB)(3),(*paramBB)(4));
+      AliInfo(Form("Setting BB parameters from OCDB (AliTPCParam): %.2g, %.2g, %.2g, %.2g, %.2g",(*paramBB)(0),(*paramBB)(1),(*paramBB)(2),(*paramBB)(3),(*paramBB)(4)));
+    }
+    else {
+      AliError("Couldn't get BB parameters from OCDB, the old default values will be used instead");
+    }
+  }
+  else {
+    AliError("Couldn't get TPC parameters");
+  }
+}
+//____________________________________________________
+void FT2::InitEnvLocal()
+{
+  AliInfo("Initializing Local Environment");
+
+  // init with environment of local ITSU generation
+  AliCDBManager* man = AliCDBManager::Instance();
+  man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+
+  man->SetSpecificStorage("GRP/*/*","alien://folder=/alice/data/2010/OCDB");
+  man->SetSpecificStorage("ITS/*/*", "alien://folder=/alice/simulation/LS1_upgrade/Ideal");
+  man->SetSpecificStorage("TPC/*/*", "alien://folder=/alice/simulation/2008/v4-15-Release/Ideal/");
+  //   man->SetSpecificStorage("GRP/GRP/Data","alien://folder=/alice/data/2010/OCDB");
+  //   man->SetSpecificStorage("ITS/Align/Data", "alien://folder=/alice/simulation/LS1_upgrade/Ideal");
+  //   man->SetSpecificStorage("ITS/Calib/RecoParam", "alien://folder=/alice/simulation/LS1_upgrade/Ideal");
+  man->SetRun(138871);
+  man->Print("");
+  /*
+    AliCDBManager* man = AliCDBManager::Instance();
+    man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
+    man->SetSpecificStorage("GRP/GRP/Data",
+    Form("local://%s",gSystem->pwd()));
+    man->SetSpecificStorage("ITS/Align/Data",
+    Form("local://%s",gSystem->pwd()));
+    man->SetSpecificStorage("ITS/Calib/RecoParam",
+    Form("local://%s",gSystem->pwd()));
+    man->SetRun(0);*/
+  //
+}
+
+//____________________________________________________
+void FT2::AddTPC(Float_t sigY, Float_t sigZ, Float_t eff,Float_t scEdge)
+{
+  // add TPC mock-up
+  //
+  const Float_t kRadLPerRow = 0.000036;
+  //
+  const Float_t kTPCInnerRadialPitch  =    0.75 ;    // cm
+  const Float_t kTPCMiddleRadialPitch =    1.0  ;    // cm
+  const Float_t kTPCOuterRadialPitch  =    1.5  ;    // cm
+  const Int_t   kTPCInnerRows            =   63 ;
+  const Int_t   kTPCMiddleRows           =   64  ;
+  const Int_t   kTPCOuterRows            =   32  ;
+  const Int_t   kTPCRows       =   (kTPCInnerRows + kTPCMiddleRows + kTPCOuterRows) ;
+  const Float_t kTPCRowOneRadius         =   85.2  ;    // cm
+  const Float_t kTPCRow64Radius          =  135.1  ;    // cm
+  const Float_t kTPCRow128Radius         =  199.2  ;    // cm
+  //
+  if (fIsTPC) {
+    printf("TPC was already added\n");
+    return;
+  }
+  fIsTPC = kTRUE;
+  fTPCSectorEdge = scEdge;
+  //
+  for ( Int_t k = 0 ; k < kTPCRows ; k++ ) {
+    //
+    Float_t rowRadius =0;
+    if (k<kTPCInnerRows) rowRadius =  kTPCRowOneRadius + k*kTPCInnerRadialPitch ;
+    else if ( k>=kTPCInnerRows && k<(kTPCInnerRows+kTPCMiddleRows) )
+      rowRadius =  kTPCRow64Radius + (k-kTPCInnerRows+1)*kTPCMiddleRadialPitch ;
+    else if (k>=(kTPCInnerRows+kTPCMiddleRows) && k<kTPCRows )
+      rowRadius = kTPCRow128Radius + (k-kTPCInnerRows-kTPCMiddleRows+1)*kTPCOuterRadialPitch ;
+    AddTPCLayer(rowRadius,kRadLPerRow,sigY,sigZ,eff);
+  }
+  //
+  fTPCHitLr.resize(fTPCLayers.size());
+}
+
+//____________________________________________________
+void FT2::AddTPCLayer(Float_t x, Float_t x2x0,Float_t sigY, Float_t sigZ, Float_t eff)
+{
+  // add single TPC layer
+  fTPCLayers.push_back(FT2TPCLayer(x,x2x0,sigY,sigZ,eff));
+}
+
+//____________________________________________________
+void FT2::PrintLayout()
+{
+  // print setup
+  if (fITS) fITS->Print("lr");
+  if (fIsTPC) {
+    printf("TPC inactive sector edge: %.3f\n",fTPCSectorEdge);
+    printf("TPC \t  R   \tx2x0\tsgRPhi\t sigZ  \t Eff  \n");
+    for (int ilr=0;ilr<(int)fTPCLayers.size();ilr++) {
+      FT2TPCLayer& lr = fTPCLayers[ilr];
+      printf("%3d\t%6.2f\t%1.4f\t",ilr,lr.x,lr.x2x0);
+      //
+      if (lr.isDead) printf("      \t");
+      else printf("%6.4f\t",lr.rphiRes);
+      //
+      if (lr.isDead) printf("      \t");
+      else printf("%6.4f\t",lr.zRes);
+      //
+      if (lr.isDead) printf("      \t");
+      else printf("%6.4f\t",lr.eff);
+      //
+      printf("\n");
+    }
+  }
+}
+
+//____________________________________________________
+Bool_t FT2::ProcessTrack(TParticle* part, AliVertex* vtx)
+{
+  // track the particle throught the setup; it must be provided in the production point
+  // then relate to vtx
+  //
+  static AliESDVertex vtx0;
+  static Bool_t first = kTRUE;
+  if (first) {
+    double vd[6] = {0};
+    vtx0.SetXYZ(&vd[0]);
+    vtx0.SetCovarianceMatrix(&vd[0]);
+    first = kFALSE;
+  }
+  if (!InitProbe(part)) {
+#if DEBUG
+    printf("Initialization failed\n");
+    part->Print();
+    fProbe.Print();
+#endif
+    return kFALSE;
+  }
+  PrepareProbe();
+  //
+  // check if there is something to reconstruct
+  if ( (fIsTPC && fNTPCHits<fMinTPCHits) || (fNITSLrHit<fMinITSLrHit)) {
+#if DEBUG
+    printf("Track hit requirements are not satisfied: Hit ITS Layers: %d (need %d) TPCHits: %d (need %d)\n",
+          fNITSLrHit,fMinITSLrHit, fNTPCHits, fIsTPC ? fMinTPCHits:0);
+    fProbe.Print();
+#endif
+    return kFALSE;
+  }
+       
+  ResetCovMat(&fProbe);
+  //
+  if (!ReconstructProbe()) {    
+#if DEBUG
+    printf("Track reconstruction failed\n");
+    fProbe.Print();
+#endif
+    return kFALSE;
+  }
+  //
+  // go to innermost radius of ITS (including beam pipe)
+  if (!PropagateToR(fITS->GetRMin(),-1, kTRUE, kFALSE, kTRUE)) {
+#if DEBUG
+    printf("Track propagatation to ITS RMin=%f failed\n",fITS->GetRMin());
+    fProbe.Print();
+#endif  
+    return kFALSE; // don't exit on faile propagation, may simply not reach this point
+  }
+  //
+  AliVertex* vtuse = vtx ? vtx : (AliVertex*)&vtx0; // if no vertex provided, relate to 0,0,0
+  //
+  fDCACov[0] = fDCACov[1] = fDCACov[2] = 0;
+  Bool_t res = fProbe.PropagateToDCA(vtuse,fBz, fITS->GetRMin(), fDCA, fDCACov);
+
+#if DEBUG
+  if (!res) {
+    printf("Track propagation to vertex failed\n");
+    vtuse->Print();
+    fProbe.Print();
+  }
+#endif  
+  //
+  /*
+    printf("Tracking done: NclITS: %d (chi2ITS=%.3f) NclTPC: %3d (chi2TPC=%.3f)\n",fNClITS,fChi2ITS,fNClTPC,fChi2TPC);
+    fProbe.Print();
+    //
+    if (res) {
+    printf("DCA: {%e %e} DCACov: {%e %e %e}\n",fDCA[0],fDCA[1],fDCACov[0],fDCACov[1],fDCACov[2]);
+    }
+    else {
+    printf("Failed to propagate to vertex\n");
+    }
+  */
+  //
+  return res;
+}
+
+//____________________________________________________
+Bool_t FT2::InitProbe(TParticle* part)
+{
+  // create probe (code copied from AliESDtrack::AliESDtrack(TParticle * part) :
+  //
+  Double_t xref;
+  Double_t alpha;
+  Double_t param[5];
+  Double_t covar[15] = {0};
+  //
+  fNTPCHits = fNClTPC = fNClITS = fNClITSFakes = fNITSLrHit = fITSPatternFake = fITSPattern = 0;
+  fChi2TPC = fChi2ITS = 0;
+  fDCA[0] = fDCA[1] = fDCACov[0] = fDCACov[1] = fDCACov[2] = 0.;
+  //
+  // Calculate alpha: the rotation angle of the corresponding local system (TPC sector)
+  alpha = part->Phi()*180./TMath::Pi();
+  if (alpha<0) alpha+= 360.;
+  if (alpha>360) alpha -= 360.;
+  //
+  Int_t sector = (Int_t)(alpha/20.);
+  alpha = 10. + 20.*sector;
+  alpha /= 180;
+  alpha *= TMath::Pi();
+  //
+  // Get the vertex of origin and the momentum
+  TVector3 ver(part->Vx(),part->Vy(),part->Vz());
+  TVector3 mom(part->Px(),part->Py(),part->Pz());
+       
+  // Rotate to the local coordinate system (TPC sector)
+  ver.RotateZ(-alpha);
+  mom.RotateZ(-alpha);
+       
+  // X of the referense plane
+  xref = ver.X();
+       
+  Int_t pdgCode = part->GetPdgCode();
+  TParticlePDG* pdgp = TDatabasePDG::Instance()->GetParticle(pdgCode);
+  Double_t charge = pdgp->Charge();
+  fProbe.fProbeMass = pdgp->Mass();
+  fProbe.fAbsPdgCode = TMath::Abs(pdgCode);
+  //
+  param[0] = ver.Y();
+  param[1] = ver.Z();
+  param[2] = TMath::Sin(mom.Phi());
+  param[3] = mom.Pz()/mom.Pt();
+  param[4] = TMath::Sign(1/mom.Pt(),charge);
+  //   AliInfo(Form("Probe Mass is: %f ",fProbe.fProbeMass));
+  // Set AliExternalTrackParam
+  fProbe.Set(xref, alpha, param, covar);
+  ResetCovMat(&fProbe);
+  return kTRUE;
+}
+
+//____________________________________________________
+Bool_t FT2::PrepareProbe()
+{
+  // propagate the probe to max allowed R of the setup
+  int nlrITS = fITS->GetNLayers(); // including passive layers
+  //
+  double xyzTmp[3];
+  //
+  for (int ilr=0;ilr<nlrITS;ilr++) {
+    AliITSURecoLayer* lr = fITS->GetLayer(ilr);
+    if (lr->IsPassive()) { // cylindric passive layer, just go to this layer
+      if (!PropagateToR(lr->GetRMax(),1,kFALSE,fSimMat,fSimMat)) return kFALSE;
+    }
+    else { // active layer, need to simulate the hit positions
+      if (!PassActiveITSLayer(lr)) return kFALSE;
+    }
+  }
+  //
+  fNTPCHits = 0;
+  if (fIsTPC) {
+    const double kTanSectH = 1.76326980708464975e-01; // tangent of  10 degree (sector half opening)
+    double xyz[3];
+    fProbe.GetXYZ(xyz);
+    Double_t alphan = TMath::ATan2(xyz[1], xyz[0]);
+    if (!fProbe.RotateParamOnly(alphan)) {
+      return kFALSE; // start in the frame with X pointing to track position
+    }
+    if (!PropagateToR(fTPCLayers[0].x-0.1,1,kFALSE,fSimMat,fSimMat)) {
+      return kFALSE; // reach inner layer
+    }
+    //
+    int sector = -1;
+    for (int ilr=0;ilr<(int)fTPCLayers.size();ilr++) {
+#if DEBUG>5
+      printf("At TPC lr %d\n",ilr);
+      fProbe.Print();
+#endif
+      FT2TPCLayer_t &tpcLr = fTPCLayers[ilr];
+      tpcLr.hitSect = -1;
+      fProbe.GetXYZ(xyzTmp);
+      double phi = TMath::ATan2(xyzTmp[1],xyzTmp[0]);
+      BringTo02Pi(phi);
+      Int_t sectNew = (Int_t)(phi*TMath::RadToDeg()/20.);
+      if (sector!=sectNew) {
+       sector = sectNew;
+       phi = 10. + 20.*sector;
+       phi /= 180;
+       phi *= TMath::Pi();
+       if (!fProbe.RotateParamOnly(phi)) {
+#if DEBUG
+         printf("TPC:Failed to rotate to alpha %.4f (sc %d)\n",phi,sector);
+         fProbe.Print();
+#endif
+         return kFALSE; // go to sector frame
+       }
+      }
+      if (!fProbe.PropagateParamOnlyTo(tpcLr.x, fBz)) {
+#if DEBUG
+       printf("TPC:Failed to go to X: %.4f\n",tpcLr.x);
+       fProbe.Print();
+#endif
+       return kFALSE; // no materials inside TPC
+      }
+      if (TMath::Abs(fProbe.GetZ())>kMaxZTPC) break; // exit from the TPC
+      double maxY = kTanSectH*tpcLr.x - fTPCSectorEdge; // max allowed Y in the sector
+      if (TMath::Abs(fProbe.GetY())>maxY) {
+#if DEBUG>3
+       printf("TPC: No hit in dead zone: %f\n",maxY);
+       fProbe.Print();
+#endif
+       continue; // in dead zone
+      }
+      //
+      // if track Snp > limit, stop propagation
+      if (TMath::Abs(fProbe.GetSnp())>fMaxSnpTPC) {
+#if DEBUG>2
+       printf("TPC: Max snp %f reached\n",fMaxSnpTPC);
+       fProbe.Print();
+#endif
+       return kFALSE;
+      }
+                       
+      // register hit
+#if DEBUG>5
+      printf("Register TPC hit at sect:%d, Y:%.4f Z:%.4f\n",sector,fProbe.GetY(),fProbe.GetZ());
+      fProbe.Print();
+#endif
+      if (tpcLr.isDead || (tpcLr.eff<1 && gRandom->Rndm()>tpcLr.eff)) continue;
+      //
+      tpcLr.hitSect = sector;
+      tpcLr.hitY = fProbe.GetY();
+      tpcLr.hitZ = fProbe.GetZ();
+      fTPCHitLr[fNTPCHits++] = ilr;
+      //
+    }
+  }
+  //
+  return kTRUE;
+}
+
+
+//________________________________________________________________________
+Bool_t FT2::ReconstructProbe()
+{
+  // reconstruct the probe
+       
+  //
+  int sect = -1;
+  fChi2TPC = 0;
+  fChi2ITS = 0;
+  fNClITS = 0;
+  fNClITSFakes = 0;
+  fITSPatternFake = 0;
+  fITSPattern = 0;
+  fNClTPC = 0;
+  fCurrITSLr = -1;     
+  //
+  if (fNTPCHits) {
+    // TPC tracking
+#if DEBUG>5
+    AliInfo(Form("Number of clusters generated on the pad rows: %i",fNTPCHits));
+#endif
+    for (int ih=fNTPCHits;ih--;) {
+#if DEBUG>5
+      AliInfo(Form("Entering TPC cluster loop: %i",ih));
+#endif
+      if(ih==0 && fUsePIDForTracking){
+       double signalTPC, p = fProbe.P(); fProbe.fTPCmomentum = p;
+
+       TH1* hdedx;
+       AliPID::EParticleType typePID=AliPID::EParticleType(2);
+
+       if(fProbe.fAbsPdgCode==11)
+         {
+           int pbin = ((TH1*)fTPCSignalElectron->Get("hMomentumAxis"))->FindBin(p);
+           hdedx = (TH1*)fTPCSignalElectron->Get(Form("electronDEDX%i",pbin));
+           typePID=AliPID::EParticleType(0);
+         }
+       else if(fProbe.fAbsPdgCode==13)
+         {
+           int pbin = ((TH1*)fTPCSignalMuon->Get("hMomentumAxis"))->FindBin(p);
+           hdedx = (TH1*)fTPCSignalMuon->Get(Form("muonDEDX%i",pbin));
+           typePID=AliPID::EParticleType(1);
+         }
+       else if(fProbe.fAbsPdgCode==211)
+         {
+           int pbin = ((TH1*)fTPCSignalPion->Get("hMomentumAxis"))->FindBin(p);
+           hdedx = (TH1*)fTPCSignalPion->Get(Form("pionDEDX%i",pbin));
+           typePID=AliPID::EParticleType(2);
+         }
+       else if(fProbe.fAbsPdgCode==321)
+         {
+           int pbin = ((TH1*)fTPCSignalKaon->Get("hMomentumAxis"))->FindBin(p);
+           hdedx = (TH1*)fTPCSignalKaon->Get(Form("kaonDEDX%i",pbin));
+           typePID=AliPID::EParticleType(3);
+         }
+       else if(fProbe.fAbsPdgCode==2212)
+         {
+           int pbin = ((TH1*)fTPCSignalProton->Get("hMomentumAxis"))->FindBin(p);
+           hdedx = (TH1*)fTPCSignalProton->Get(Form("protonDEDX%i",pbin));
+           typePID=AliPID::EParticleType(4);
+         }
+       else {hdedx=0;AliFatal("PDC Code %4.f cannot be treated in the code - This case should not be possible here!");}
+#if DEBUG>5
+       AliInfo(Form("PDG code of Probe: %4.f",fProbe.fAbsPdgCode));
+       AliInfo(Form("Momentum of Probe: %f",fProbe.fTPCmomentum));
+#endif
+       if(hdedx->Integral()>0){ signalTPC = ((TH1*)hdedx)->GetRandom();} //30
+       else{
+                                       
+         Double_t mean = fPIDResponse->GetTPCResponse().GetExpectedSignal(&fProbe,typePID,AliTPCPIDResponse::kdEdxDefault,kFALSE,kFALSE);
+         fProbe.fTPCSignalN = (UShort_t)mean;
+         Double_t sigma = fPIDResponse->GetTPCResponse().GetExpectedSigma(&fProbe,typePID,AliTPCPIDResponse::kdEdxDefault,kFALSE,kFALSE);
+         signalTPC = (40./50)*gRandom->Gaus(mean,sigma); // 40. or 33. ?
+#if DEBUG>5
+         AliInfo(Form("TPC Signal for %i from scaling %f",typePID,signalTPC));
+#endif
+       }
+       fProbe.fTPCSignal = signalTPC;
+       fProbe.fTPCSignalN = (UShort_t)signalTPC;
+
+       Double_t prob[AliPID::kSPECIES]={0.};
+       GetComputeTPCProbability(&fProbe,AliPID::kSPECIES,prob);
+                               
+       Float_t max=0.,min=1.e9;
+       Int_t pid=-1;
+       for (Int_t i=0; i<AliPID::kSPECIES; ++i) {
+         if (prob[i]>max) {pid=i; max=prob[i];}
+         if (prob[i]<min) min=prob[i];
+       }
+                               
+       if (pid>AliPID::kSPECIES-1 || (min>=max)) pid = AliPID::kPion;
+                               
+#if DEBUG>5
+       for(Int_t i=0;i<AliPID::kSPECIES;i++){
+         AliInfo(Form("Pid Prob : %f",prob[i]));
+       }
+#endif
+       Int_t pdgCode = 0;
+       if(pid==0) pdgCode = 11;
+       if(pid==1) pdgCode = 13;
+       if(pid==2) pdgCode = 211;
+       if(pid==3) pdgCode = 321;
+       if(pid==4) pdgCode = 2212;
+       if(pdgCode==0) AliFatal("This should not have happened");
+#if DEBUG>5
+       AliInfo(Form("Probe mass for tracking: %f",fProbe.fProbeMass));
+#endif
+       if(pid==0) fProbe.fAbsPdgCodeForTracking = 211;                 // e(11)        --> track as pion due to bug in fMC
+       if(pid==1) fProbe.fAbsPdgCodeForTracking = 211;                 // mu(13)       --> track as pion due to bug in fMC
+       if(pid==2) fProbe.fAbsPdgCodeForTracking = 211;                 // pi(211)
+       if(pid==3) fProbe.fAbsPdgCodeForTracking = 321;                 // ka(321)
+       if(pid==4) fProbe.fAbsPdgCodeForTracking = 2212;                // pr(2212)
+       if(fProbe.fAbsPdgCodeForTracking==0) AliFatal("This should not have happened");
+       TParticlePDG* pdgp = TDatabasePDG::Instance()->GetParticle(fProbe.fAbsPdgCodeForTracking);
+       fProbe.fProbeMass = pdgp->Mass();
+#if DEBUG>5
+       AliInfo(Form("PID for tracking: %i with probe mass %f",pid,fProbe.fProbeMass));
+#endif
+      }
+                       
+      FT2TPCLayer_t &tpcLr = fTPCLayers[fTPCHitLr[ih]];
+      if (tpcLr.hitSect!=sect) { // rotate to hit sector
+       sect = tpcLr.hitSect;
+       if (!fProbe.Rotate( TMath::DegToRad()*(10.+20.*sect))) return kFALSE;
+      }
+      // go to the padrow
+      if (!fProbe.PropagateTo(tpcLr.x,fBz)) return kFALSE;
+      if (tpcLr.x2x0>0 && !fProbe.CorrectForMeanMaterial(tpcLr.x2x0,0,fProbe.fProbeMass) ) return kFALSE;
+      //
+      if (tpcLr.isDead) continue;
+      double chi = UpdateKalman(&fProbe,tpcLr.hitY,tpcLr.hitZ,tpcLr.rphiRes,tpcLr.zRes,kTRUE,kFALSE);
+      if (chi<0) return kFALSE;
+      fChi2TPC += chi;
+      fNClTPC++;
+    }
+    // go to ITS/TPC matching R, accounting for TGeo materials
+    if (!PropagateToR(fITS->GetRITSTPCRef(),-1,kTRUE, kFALSE, kTRUE)) return kFALSE;
+  }
+  //
+  // 
+#if DEBUG
+  printf("Outward kalman results:\n");
+  for (int ilr=0;ilr<fITS->GetNLayersActive();ilr++) {
+    if (GetKalmanOut(ilr).GetSigmaY2()<=0) continue;
+    printf("Extrap to Lr %d : ",ilr); GetKalmanOut(ilr).Print();
+  }
+#endif
+  //
+  // ITS tracking
+  for (int ilr=fITS->GetNLayersActive();ilr--;) {
+    fCurrITSLr = ilr;
+#if DEBUG>5
+    AliInfoF("Entering ITS cluster loop on lr%d: %d clusters",ilr,fNITSHits[ilr]);
+#endif
+    if (!fNITSHits[ilr]) continue;
+    // for 2 hits (overlap) chose only 1
+    int iht = fNITSHits[ilr]==1 ? 0 : (gRandom->Rndm()>0.5 ? 0:1);
+    AliITSURecoSens* sens = fITSSensHit[ilr][iht];
+    if (!fProbe.Rotate(sens->GetPhiTF())) {
+#if DEBUG
+      printf("Failed to rotate to sensor phi %f at lr %d\n",sens->GetPhiTF(),ilr); 
+      fProbe.Print();
+      return kFALSE;
+#endif
+    }
+    if (!PropagateToX(sens->GetXTF(),-1,kTRUE, kFALSE, kTRUE)) return kFALSE; // account for materials
+    double chi = UpdateKalman(&fProbe,fITSHitYZ[ilr][iht][0],fITSHitYZ[ilr][iht][1],fSigYITS,fSigZITS,kTRUE,kTRUE);
+ #if DEBUG>5
+    AliInfoF("Updated at ITS Lr%d, chi2:%f NclITS:%d Fakes: %d",ilr,chi,fNClITS,fNClITSFakes);
+#endif   
+    if (chi<0) return kFALSE;
+    fChi2ITS += chi;
+    fNClITS++;
+    fITSPattern |= 0x1<<ilr;
+  }
+  //
+  return kTRUE;
+}
+
+//________________________________________________________________________
+Double_t FT2::UpdateKalman(AliExternalTrackParam* trc, double y,double z,double sigY,double sigZ,Bool_t randomize, Bool_t fake)
+{
+  // calman update, return chi2 increment or -1 on failure
+  //
+  if (randomize) {
+    double ry,rz;
+    gRandom->Rannor(ry,rz);
+    y += sigY*ry;
+    z += sigZ*rz;
+  }
+  double meas[2] = {y,z};
+  double measErr2[3] = {sigY*sigY,0,sigZ*sigZ};
+  //
+  // do we want to simulate fakes?
+  if (fake && fdNdY>0) {
+    double err2a,err2b;   
+    double trCov[3],trPos[2];
+    if (fCurrITSLr>=0 && fKalmanOutward[fCurrITSLr].GetSigmaY2()>0) { // weight with outward calman track
+      AliExternalTrackParam trJoint = fKalmanOutward[fCurrITSLr];
+      if (TMath::Abs(trJoint.GetAlpha()-trc->GetAlpha())>1e-6 && !trJoint.Rotate(trc->GetAlpha())) {
+#if DEBUG
+       AliInfoF("Failed to rotate outward Kalman track to alpha=%f",trc->GetAlpha());
+       trJoint.Print();
+       return -1;
+#endif 
+      }
+      if (TMath::Abs(trJoint.GetX()-trc->GetX())>1e-4 && !trJoint.PropagateTo(trc->GetX(),fBz)) {
+#if DEBUG
+       AliInfoF("Failed to propagate outward Kalman track to X=%f",trc->GetX());
+       trJoint.Print();
+       return -1;
+#endif         
+      }
+      double inwPos[2] = {trc->GetY(),trc->GetZ()};
+      double inwCov[3] = {trc->GetSigmaY2(),trc->GetSigmaZY(),trc->GetSigmaZ2()};
+      if (!trJoint.Update(inwPos,inwCov)) {
+#if DEBUG
+       AliInfo("Failed to update outward Kalman track by inward one");
+       trJoint.Print();
+       trc->Print();
+       return -1;
+#endif
+      }
+#if DEBUG//>5
+      printf("Errors on Lr %d: \n",fCurrITSLr);
+      printf("Inward  :  {%+e,%+e,%+e} {%+e,%+e}\n",trc->GetSigmaY2(),trc->GetSigmaZY(),trc->GetSigmaZ2(), trc->GetY(),trc->GetZ());
+      printf("Outward :  {%+e,%+e,%+e} {%+e,%+e}\n",fKalmanOutward[fCurrITSLr].GetSigmaY2(),
+            fKalmanOutward[fCurrITSLr].GetSigmaZY(),fKalmanOutward[fCurrITSLr].GetSigmaZ2(),
+            fKalmanOutward[fCurrITSLr].GetY(),fKalmanOutward[fCurrITSLr].GetZ());
+      printf("Weighted:  {%+e,%+e,%+e} {%+e,%+e}\n",trJoint.GetSigmaY2(),trJoint.GetSigmaZY(),trJoint.GetSigmaZ2(),
+            trJoint.GetY(),trJoint.GetZ());
+#endif
+      trPos[0] = trJoint.GetY();
+      trPos[1] = trJoint.GetZ();
+      trCov[0] = trJoint.GetSigmaY2();
+      trCov[1] = trJoint.GetSigmaZY();      
+      trCov[2] = trJoint.GetSigmaZ2();
+    }
+    else {
+      trPos[0] = trc->GetY();
+      trPos[1] = trc->GetZ();
+      trCov[0] = trc->GetSigmaY2();
+      trCov[1] = trc->GetSigmaZY();
+      trCov[2] = trc->GetSigmaZ2();
+    }
+    //
+    if (!DiagonalizeErrors(trCov,err2a,err2b)) {
+#if DEBUG
+      printf("Failed to diagonalize erros\n"); trc->Print();
+#endif
+      return -1;
+    }
+    double xyz[3]; 
+    trc->GetXYZ(xyz);
+    double rho = HitDensity(xyz[0]*xyz[0]+xyz[1]*xyz[1],trc->GetTgl());
+    //
+    // prob. of good hit (http://rnc.lbl.gov/~wieman/HitFinding2D.htm)
+    Double_t sx2, sy2, probFake;
+    sx2 = 2 * TMath::Pi()*err2a*rho;
+    sy2 = 2 * TMath::Pi()*err2b*rho;
+    probFake = 1.-TMath::Sqrt(1./((1+sx2)*(1+sy2)));
+#if DEBUG
+    printf("DiagErr: %e %e Rho:%e PFake: %e\n",err2a,err2b, rho, probFake);
+#endif
+    if (gRandom->Rndm()<probFake) { // need to fake the hit
+      BiasAsFake(meas, trPos, trCov);
+      if (fCurrITSLr>=0) {
+       fITSPatternFake |= 0x1<<fCurrITSLr;
+       fNClITSFakes++;
+      }
+      /*
+       printf("rho:%f errs: %f %f -> %f -> %d r=%f\n",rho,err2a,err2b,probFake,fNClITSFakes, 
+       TMath::Sqrt(xyz[0]*xyz[0]+xyz[1]*xyz[1]));
+      */
+    }
+  }
+  //
+  double chi2 = trc->GetPredictedChi2(meas,measErr2);
+#if DEBUG>5
+  if (fake) {
+    printf("Updating {%e %e}/{%e %e} {%e %e %e} Ncl:%d NclF:%d -> %f\n", meas[0],meas[1],y,z,
+          measErr2[0],measErr2[1],measErr2[2],fNClITS,fNClITSFakes,chi2); 
+    trc->Print();
+  }
+#endif
+  //
+  if (!trc->Update(meas,measErr2)) {
+#if DEBUG
+    printf("Failed to Update {%e %e}/{%e %e} {%e %e %e} Ncl:%d NclF:%d\n", meas[0],meas[1],y,z,
+          measErr2[0],measErr2[1],measErr2[2],fNClITS,fNClITSFakes); 
+    trc->Print();
+#endif
+     return -1;
+  }
+  return chi2;
+}
+
+//________________________________________________________________________
+void FT2::ResetCovMat(AliExternalTrackParam* trc)
+{
+  // assign huge errors to probe at max.radius
+  enum {kY,kZ,kSnp,kTgl,kPtI};              // track parameter aliases
+  enum {kY2,kYZ,kZ2,kYSnp,kZSnp,kSnp2,kYTgl,kZTgl,kSnpTgl,kTgl2,kYPtI,kZPtI,kSnpPtI,kTglPtI,kPtI2}; // cov.matrix aliases
+  const double kLargeErr2Coord = 5*5;
+  const double kLargeErr2Dir = 0.7*0.7;
+  const double kLargeErr2PtI = 30.5*30.5;
+  double *trPars = (double*)trc->GetParameter();
+  double *trCov  = (double*)trc->GetCovariance();
+  //
+  for (int ic=15;ic--;) trCov[ic] = 0.;
+  trCov[kY2]   = trCov[kZ2]   = kLargeErr2Coord;
+  trCov[kSnp2] = trCov[kTgl2] = kLargeErr2Dir;
+  trCov[kPtI2] = kLargeErr2PtI*trPars[kPtI]*trPars[kPtI];
+  trc->CheckCovariance();
+  //
+}
+
+
+//________________________________________________________________________
+Bool_t FT2::PropagateToR(double r, int dir,
+                        Bool_t propErr, // propagate errors
+                        Bool_t simMat,  // simulate material effects
+                        Bool_t useTGeo)
+{
+  // propagate track to given R (not X!)
+  double xTgt;
+  if (!fProbe.GetXatLabR(r, xTgt, fBz, dir)) {
+#if DEBUG
+    printf("PropagateToR failed in GetXatLabR(%.3f,%.3f, %.3f, %d)\n", r,xTgt, fBz, dir);
+    fProbe.Print();
+#endif
+    return kFALSE;
+  }
+  double dx = xTgt-fProbe.GetX();
+  if (dir*dx<0) return kTRUE; // probe is already above given R
+  return PropagateToX(xTgt,dir,propErr,simMat,useTGeo);
+}
+
+//________________________________________________________________________
+Bool_t FT2::PropagateToX(double xTgt, int dir,
+                        Bool_t propErr, // propagate errors
+                        Bool_t simMat,  // simulate material effects
+                        Bool_t useTGeo)
+{
+  // propagate track to given X
+  //
+  Bool_t doKalmanOut = kFALSE;
+  if (simMat && dir>0 && fUseKalmanOut && xTgt<fITS->GetRMax()) {
+    doKalmanOut = kTRUE;
+  }
+  double maxStep = useTGeo ? fgMaxStepTGeo : 1e6;
+  //
+  Double_t xyz0[3],xyz1[3],param[7];
+  //
+  if (useTGeo) {
+    fProbe.GetXYZ(xyz0);
+  }
+  while ( dir*(xTgt-fProbe.GetX()) > kTrackToler) {
+    Double_t step = dir*TMath::Min(TMath::Abs(xTgt-fProbe.GetX()), maxStep);
+    Double_t x    = fProbe.GetX()+step;
+    Bool_t res = (propErr||doKalmanOut) ? fProbe.PropagateTo(x,fBz) : fProbe.PropagateParamOnlyTo(x,fBz);
+    if (!res)  {
+#if DEBUG
+      printf("Failed Prop PropagateToX(%.1f,%d,%d,%d,%d) doKalmanOut=%d",xTgt,dir,propErr,simMat,useTGeo,doKalmanOut);
+      fProbe.Print();
+#endif
+      return kFALSE;
+    }
+    if (useTGeo) {
+      fProbe.GetXYZ(xyz1);
+      AliTrackerBase::MeanMaterialBudget(xyz0,xyz1,param);
+      Double_t xrho=param[0]*param[4], xx0=param[1];
+      if (dir>0) xrho = -xrho;
+      if (simMat) {
+       if (!ApplyMSEloss(xx0,xrho)) {
+#if DEBUG
+         printf("Failed ApplyMSEloss(%f,%f) in PropagateToX(%.1f,%d,%d,%d,%d)  doKalmanOut=%d",
+                xx0,xrho,xTgt,dir,propErr,simMat,useTGeo,doKalmanOut);
+         fProbe.Print();
+#endif
+         return kFALSE;
+       }
+       if (doKalmanOut) { // special mode for outward kalman: we don't want to affect params by materials since
+         // it is already done, onle cov.matrix need to be updated
+         AliExternalTrackParam trTmp = fProbe;
+         if (!trTmp.CorrectForMeanMaterial(xx0,xrho,fProbe.fProbeMass,kTRUE)) {
+#if DEBUG
+           printf("Failed MatCorr(%f %f) in PropagateToX(%.1f,%d,%d,%d,%d) for doKalmanOut",xx0,xrho,xTgt,dir,propErr,simMat,useTGeo);
+             trTmp.Print();
+#endif
+         }
+         memcpy((double*)fProbe.GetCovariance(),trTmp.GetCovariance(),15*sizeof(double)); // restore parameter
+       }
+      }
+      else if (!fProbe.CorrectForMeanMaterial(xx0,xrho,fProbe.fProbeMass,kTRUE)) {
+#if DEBUG
+       printf("Failed CorrMeanMat(%f,%f) in PropagateToX(%.1f,%d,%d,%d,%d)",xx0,xrho,xTgt,dir,propErr,simMat,useTGeo);
+       fProbe.Print();
+#endif
+       return kFALSE;
+      }
+    }
+    memcpy(xyz0,xyz1,3*sizeof(double));
+    //
+  }
+  //
+  return kTRUE;
+}
+
+//______________________________________________________________
+Bool_t FT2::ApplyMSEloss(double x2X0, double xrho)
+{
+  // simulate random modification of track params due to the MS
+  if (x2X0<=0) return kTRUE;
+  //  printf("BeforeMAT: "); fProbe.Print();
+  double alpha = fProbe.GetAlpha(); // store original alpha
+  double mass2 = fProbe.fProbeMass*fProbe.fProbeMass;
+  //
+  double snp = fProbe.GetSnp();
+  double dip = fProbe.GetTgl();
+  Double_t angle=TMath::Sqrt((1.+ dip*dip)/((1-snp)*(1.+snp)));
+  x2X0 *= angle;
+  //
+  static double covCorr[15],covDum[21]={0};
+  static double mom[3],pos[3];
+  double *cov = (double*) fProbe.GetCovariance();
+  memcpy(covCorr,cov,15*sizeof(double));
+  fProbe.GetXYZ(pos);
+  fProbe.GetPxPyPz(mom);
+  double pt2 = mom[0]*mom[0]+mom[1]*mom[1];
+  double pt = TMath::Sqrt(pt2);
+  double ptot2 = pt2 + mom[2]*mom[2];
+  double ptot  = TMath::Sqrt(ptot2);
+  double beta = ptot/TMath::Sqrt(ptot2 + mass2);
+  double sigth = TMath::Sqrt(x2X0)*0.014/(ptot*beta);
+  //
+  // a la geant
+  double phiSC = gRandom->Rndm()*TMath::Pi();
+  double thtSC = gRandom->Gaus(0,1.4142*sigth);
+  //  printf("MS phi: %+.5f tht: %+.5f\n",phiSC,thtSC);
+  double sn = TMath::Sin(thtSC);
+  double dx = sn*TMath::Sin(phiSC);
+  double dy = sn*TMath::Cos(phiSC);
+  double dz = TMath::Cos(thtSC);
+  double v[3];
+  //  printf("Before: %+.3e %+.3e %+.3e | MS: %+.3e %+.3e\n",mom[0],mom[1],mom[2],thtSC,phiSC);
+  for (int i=3;i--;) mom[i] /= ptot;
+  double vmm = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
+  if (!IsZero(pt)) {
+    double pd1 = mom[0]/vmm;
+    double pd2 = mom[1]/vmm;
+    v[0] = pd1*mom[2]*dx - pd2*dy + mom[0]*dz;
+    v[1] = pd2*mom[2]*dx + pd1*dy + mom[1]*dz;
+    v[2] = -vmm*dx                + mom[2]*dz;
+  }
+  else {
+    v[0] = dx;
+    v[1] = dy;
+    v[2] = dz*TMath::Sign(1.,mom[2]);
+  }
+  //
+  // account for eloss
+  Double_t etot=TMath::Sqrt(ptot2 + mass2);
+  Double_t dE=fProbe.BetheBlochSolid(ptot/fProbe.fProbeMass)*xrho;
+  //  printf("X:%e E=%e dE=%e (xrho:%e)-> %e | ptot=%e M2 = %e\n",fProbe.GetX(),etot,dE,xrho,etot+dE,ptot,mass2);
+  if ( TMath::Abs(dE) > 0.3*ptot ) {
+#if DEBUG>1
+    printf("StopEloss ");fProbe.Print();
+#endif
+    return kFALSE;
+  } //30% energy loss is too much!
+  etot += dE;
+  if (etot<fProbe.fProbeMass+1./kRidiculous) {
+#if DEBUG>1
+    printf("StopEloss ");fProbe.Print();
+#endif
+    return kFALSE;
+  } // stopped
+  ptot = TMath::Sqrt(etot*etot - mass2);
+  //
+  double nrm = TMath::Sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
+  //  printf("before :%+e %+e %+e  || %+e %+e %+e %+e\n",mom[0],mom[1],mom[2],  sigth, x2X0, pt, beta);
+  //  fProbe.Print();
+  // direction cosines -> p
+  for (int i=3;i--;) mom[i] = ptot*v[i]/nrm;
+  //  printf("After : %+.3e %+.3e %+.3e\n",mom[0],mom[1],mom[2]);
+  fProbe.Set(pos,mom,covDum,fProbe.Charge());
+  //
+  fProbe.RotateParamOnly(alpha);
+  memcpy(cov,covCorr,15*sizeof(double));
+       
+  //  printf("AfterMAT: "); fProbe.Print();
+  return kTRUE;
+  //
+}
+
+//__________________________________________________________________
+Bool_t FT2::GetRoadWidth(AliITSURecoLayer* lrA,double *pr)
+{
+  static AliExternalTrackParam sc;   // seed copy for manipulations
+#if DEBUG>5
+  printf(">>RW at Lr%d\n",lrA->GetActiveID());
+#endif
+  sc = fProbe;
+  double xt;
+  if (!sc.GetXatLabR(lrA->GetRMin(),xt,fBz,1)) return kFALSE;
+  if (!sc.PropagateParamOnlyTo(xt,fBz)) return kFALSE;
+  sc.GetXYZ(&pr[AliITSUTrackerGlo::kTrXIn]);
+  pr[AliITSUTrackerGlo::kTrPhiIn] = TMath::ATan2(pr[AliITSUTrackerGlo::kTrYIn],pr[AliITSUTrackerGlo::kTrXIn]);
+  if (!sc.RotateParamOnly(pr[AliITSUTrackerGlo::kTrPhiIn])) return kFALSE; // go to the frame of the entry point into the layer
+  BringTo02Pi(pr[AliITSUTrackerGlo::kTrPhiIn]);
+  double dr  = lrA->GetDR();                              // approximate X dist at the inner radius
+  if (!sc.GetXYZAt(sc.GetX()-dr, fBz, pr + AliITSUTrackerGlo::kTrXOut)) {
+    // special case: track does not reach inner radius, might be tangential
+    double r = sc.GetD(0,0,fBz);
+    double x;
+    if (!sc.GetXatLabR(r,x,fBz,-1)) return kFALSE;
+    dr = Abs(sc.GetX() - x);
+    if (!sc.GetXYZAt(x, fBz, pr + AliITSUTrackerGlo::kTrXOut)) return kFALSE;
+  }
+  //
+  pr[AliITSUTrackerGlo::kTrPhiOut] = ATan2(pr[AliITSUTrackerGlo::kTrYOut],pr[AliITSUTrackerGlo::kTrXOut]);
+  BringTo02Pi(pr[AliITSUTrackerGlo::kTrPhiOut]);
+  double sgy = 1e-6; // dummy spread
+  double sgz = 1e-6; // dummy spread
+  double phi0  = MeanPhiSmall(pr[AliITSUTrackerGlo::kTrPhiOut],pr[AliITSUTrackerGlo::kTrPhiIn]);
+  double dphi0 = DeltaPhiSmall(pr[AliITSUTrackerGlo::kTrPhiOut],pr[AliITSUTrackerGlo::kTrPhiIn]);
+  //
+  pr[AliITSUTrackerGlo::kTrPhi0] = phi0;
+  pr[AliITSUTrackerGlo::kTrZ0]   = 0.5*(pr[AliITSUTrackerGlo::kTrZOut]+pr[AliITSUTrackerGlo::kTrZIn]);
+  dphi0 += sgy/lrA->GetR();
+  pr[AliITSUTrackerGlo::kTrDPhi] =  dphi0<PiOver2() ? dphi0 : PiOver2();
+  pr[AliITSUTrackerGlo::kTrDZ]   = 0.5*Abs(pr[AliITSUTrackerGlo::kTrZOut]-pr[AliITSUTrackerGlo::kTrZIn])   + sgz;
+  //
+#if DEBUG>5
+  printf("<<RW at Lr%d: phi0: %e+-%e z0: %e+-%e\n",lrA->GetActiveID(),pr[AliITSUTrackerGlo::kTrPhi0],
+        pr[AliITSUTrackerGlo::kTrDPhi],pr[AliITSUTrackerGlo::kTrZ0],pr[AliITSUTrackerGlo::kTrDZ]);
+#endif
+  return kTRUE;
+  //
+}
+
+//__________________________________________________________
+Bool_t FT2::PassActiveITSLayer(AliITSURecoLayer* lr)
+{
+  // pass the layer and register the hits
+  Double_t trImpData[AliITSUTrackerGlo::kNTrImpData];
+  AliITSURecoSens *hitSens[AliITSURecoLayer::kMaxSensMatching];
+  AliITSUGeomTGeo* gm = fITS->GetGeom();
+  //
+  int lrAID = lr->GetActiveID();
+  ((double*)fKalmanOutward[lrAID].GetCovariance())[0] = -1; // invalidate outward track
+  //
+  if (!GetRoadWidth(lr, trImpData)) return kFALSE;
+  fNITSHits[lrAID] = 0;
+  const AliITSsegmentation* segm = gm->GetSegmentation(lrAID);
+  int nsens = lr->FindSensors(&trImpData[AliITSUTrackerGlo::kTrPhi0], hitSens);
+  int idxHit[2] = {0,1};
+  //
+#if DEBUG>5
+  printf("Will test %d sensors on lr %d\n",nsens,lr->GetActiveID());
+  for (int isn=0;isn<nsens;isn++) hitSens[isn]->Print();
+  fProbe.Print();
+#endif
+  //
+  if (nsens>1) { // in case of overlaping sensors need to traverse them in increasing R order
+    double r2hit[AliITSURecoLayer::kMaxSensMatching];
+    AliExternalTrackParam tmpt;
+    for (int isn=0;isn<nsens;isn++) {
+      AliITSURecoSens* sens =  hitSens[isn]; // estimate hit point, we need them ordered in R
+      tmpt = fProbe;
+      if (!tmpt.RotateParamOnly(sens->GetPhiTF())) {
+#if DEBUG
+       printf("testFailed to rotate for sens %d ",isn);sens->Print();
+       tmpt.Print();
+#endif
+       return kFALSE;
+      }
+      double x = sens->GetXTF(), y;
+      if (!tmpt.GetYAt(x,fBz,y)) {
+#if DEBUG
+       printf("tesdFailed to do GetYAt for sensor %d ",isn);sens->Print();
+       tmpt.Print();
+#endif
+       return kFALSE;
+      }
+      r2hit[isn] = x*x+y*y;
+    }
+    if (r2hit[0]>r2hit[1]) {idxHit[0]=1;idxHit[1]=0;} // only 2 hits per layer are possible
+    nsens=2;
+  }
+  //
+  for (int isn=0;isn<nsens;isn++) {
+    AliITSURecoSens* sens =  hitSens[idxHit[isn]];
+    Bool_t res = fUseKalmanOut ? fProbe.Rotate(sens->GetPhiTF()) : fProbe.RotateParamOnly(sens->GetPhiTF());
+    if (!res) {
+#if DEBUG
+      printf("Failed to rotate for sens %d ",isn); sens->Print();
+      fProbe.Print();
+#endif
+      return kFALSE;
+    }
+    if (!PropagateToX(sens->GetXTF(),1,kFALSE,fSimMat,fSimMat)) {
+      return kFALSE;
+    }
+    const TGeoHMatrix* mt = gm->GetMatrixT2L(sens->GetID());
+    double xyzL[3],xyzT[3] = {fProbe.GetX(),fProbe.GetY(),fProbe.GetZ()};
+    mt->LocalToMaster(xyzT,xyzL);
+    int ix,iz;
+    if (!segm->LocalToDet(xyzL[0],xyzL[2],ix,iz)) {
+#if DEBUG>5
+      double xyzt[3];
+      double ph = TMath::ATan2(xyzt[1],xyzt[0]);
+      BringTo02Pi(ph);
+      printf("NoHit with XZloc %+.4f %+.4f, phi %.4f\n",xyzL[0],xyzL[2],ph);
+      fProbe.Print();
+      sens->Print();
+      segm->Print();
+#endif
+      continue; // not in the active zone
+    }
+    // register hit
+    int nh = fNITSHits[lrAID];
+    fITSSensHit[lrAID][nh] = sens;
+    fITSHitYZ[lrAID][nh][0] = fProbe.GetY();
+    fITSHitYZ[lrAID][nh][1] = fProbe.GetZ();
+    fNITSHits[lrAID]++;
+    if (fNITSHits[lrAID]) fNITSLrHit++;
+#if DEBUG>5
+    printf("Register hit %d at lr%d\n",fNITSHits[lrAID],lr->GetActiveID());
+    fProbe.Print();
+#endif
+    // emulate outward Kalman track
+    if (fUseKalmanOut) {
+      AliExternalTrackParam& trKO = fKalmanOutward[lrAID];
+      if (((double*)trKO.GetCovariance())[0]<0) {
+       trKO = fProbe;
+       AliExternalTrackParam trKOUp = fProbe;
+       double chi = UpdateKalman(&trKOUp,fITSHitYZ[lrAID][nh][0],fITSHitYZ[lrAID][nh][1],
+                                 fSigYITS,fSigZITS,kTRUE,kFALSE);
+       if (chi<0) {
+         printf("Failed on emulating outward Kalman on lr %d\n",lrAID);
+         trKOUp.Print();
+         return kFALSE;
+       }
+       // assign updated error matrix to fProbe
+       memcpy((double*)fProbe.GetCovariance(),trKOUp.GetCovariance(),15*sizeof(double));
+      } // update for outward Kalman
+    }
+  }
+  //
+  return kTRUE;
+}
+//__________________________________________________________
+AliPIDResponse::EDetPidStatus FT2::GetComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
+{
+  // Similar as AliPIDResponse
+  // Compute PID response for the TPC
+  //
+  // set flat distribution (no decision)
+  for (Int_t k=0; k<nSpecies; k++) p[k]=1./nSpecies;
+       
+       
+  Double_t dedx=track->GetTPCsignal();
+  Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
+       
+  Double_t bethe = 0.;
+  Double_t sigma = 0.;
+       
+  for (Int_t j=0; j<nSpecies; j++) {
+    AliPID::EParticleType type=AliPID::EParticleType(j);
+        
+    bethe=fPIDResponse->GetTPCResponse().GetExpectedSignal(track,type,AliTPCPIDResponse::kdEdxDefault,kFALSE, kFALSE);
+    sigma=fPIDResponse->GetTPCResponse().GetExpectedSigma(track, type,AliTPCPIDResponse::kdEdxDefault,kFALSE, kFALSE);
+#if DEBUG>5
+    AliInfo(Form("Species Hypothesis %i - dEdx: %f - Bethe :%f +/- %f",j,dedx,bethe,sigma));
+#endif
+    if (TMath::Abs(dedx-bethe) > 5.*sigma) {
+      p[j]=TMath::Exp(-0.5*5.*5.)/sigma;
+    } else {
+      p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
+      mismatch=kFALSE;
+    }
+  }
+  if (mismatch){
+    for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+  }
+  return AliPIDResponse::kDetPidOk;
+}
+
+//__________________________________________________________
+Bool_t FT2::DiagonalizeErrors(const double* cov, double &sy2d, double &sz2d) 
+{
+  // diagonalize cov. matrix
+  double dd = cov[0]-cov[2]; 
+  dd = TMath::Sqrt(dd*dd + 4.*cov[1]*cov[1]);
+  double sd = cov[0]+cov[2];
+  sy2d = 0.5*(sd - dd);
+  sz2d = 0.5*(sd + dd);
+  if (sy2d<=0 || sz2d<=0) return kFALSE;
+  return kTRUE;
+}
+
+//_____________________________________________
+Double_t FT2::HitDensity(double r2, double tgl) const
+{
+  // hit density from single collision at radius r
+  double den = fdNdY/(2.*TMath::Pi()*r2);
+  return den/TMath::Sqrt(1 + tgl*tgl);
+}
+
+//_____________________________________________
+Bool_t FT2::BiasAsFake(double yz[2], const double* extyz, const double *cov) const
+{
+  // assign instead of the "true" yz hit position a fake position such that it
+  // will have better chi2 wrt the extrapolation point
+  Double_t r00=cov[0], r01=cov[1], r11=cov[2];
+  Double_t det=r00*r11 - r01*r01;
+  if (TMath::Abs(det) < 1e-16) return kFALSE;
+  Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
+  double dy0 = yz[0]-extyz[0];
+  double dz0 = yz[1]-extyz[1];
+  double d2max = dy0*dy0*r00+dz0*dz0*r11+2.*dy0*dz0*r01;
+  dy0 = TMath::Sqrt(d2max*cov[0]);
+  dz0 = TMath::Sqrt(d2max*cov[2]);
+  double d2 = 1e9;
+  //  printf("d2max: %f dy:%f dz:%f\n",d2max,dy0,dz0);
+  double dy(dy0),dz(dz0);
+  while(d2>=d2max) {
+    dy = (gRandom->Rndm()-0.5)*2*dy0;
+    dz = (gRandom->Rndm()-0.5)*2*dz0;
+    d2 = dy*dy*r00+dz*dz*r11+2.*dy*dz*r01;
+  }
+  //
+#if DEBUG>1
+  printf("AddFake DY:%f DZ:%f Chi2Max:%f Chi2F:%f\n",dy,dz, d2max,d2);
+#endif
+  yz[0] = extyz[0]+dy;
+  yz[1] = extyz[1]+dz;
+  return kTRUE;
+}
diff --git a/ITS/UPGRADE/FT2/FT2.h b/ITS/UPGRADE/FT2/FT2.h
new file mode 100644 (file)
index 0000000..853d6c7
--- /dev/null
@@ -0,0 +1,172 @@
+#ifndef FT2_H
+#define FT2_H
+#include <TObject.h>
+#include "AliExternalTrackParam.h"
+#include "AliITSUAux.h"
+#include "AliESDpid.h"
+#include "AliPIDResponse.h"
+
+//#define DEBUG 11
+//#define DEBUG 1
+
+class AliITSUReconstructor;
+class AliITSURecoDet;
+class AliITSURecoLayer;
+class AliITSURecoSens;
+class TParticle;
+class AliVertex;
+
+const float kRidiculous = 999999;
+const float kMaxZTPC = 250;
+const double kTrackToler = 1.e-6; // tracking tolerance
+
+class FTProbe : public AliExternalTrackParam
+{
+ public:
+  FTProbe(){};
+  virtual ~FTProbe(){};
+  virtual Double_t GetTPCmomentum() const {return fTPCmomentum;}
+  virtual Double_t GetTPCsignal() const {return fTPCSignal;}
+  virtual UShort_t GetTPCsignalN() const {return fTPCSignalN;}
+  virtual Int_t GetTPCTrackingPID() const {return fAbsPdgCodeForTracking;}
+       
+ protected:
+  Double_t fProbeMass; // true mass
+  Double_t fTPCmomentum;       // momentum after TPC reconstruction
+  Double_t fTPCSignal; // TPC signal
+  UShort_t fTPCSignalN;
+  Double_t fAbsPdgCode;        // pdg code of particle
+  Int_t fAbsPdgCodeForTracking;
+  ClassDef(FTProbe,1)
+    };
+
+
+class FT2 : public TObject
+{
+ public:
+  enum {kMaxITSLr=7, kMaxHitPerLr=2};
+  struct FT2TPCLayer {
+  FT2TPCLayer(float xr=0,float x2x=0,float resRPhi=0,float resZ=0,float effL=0) :
+    x(xr),x2x0(x2x),rphiRes(resRPhi),zRes(resZ),eff(effL),isDead(resRPhi>kRidiculous||resZ>kRidiculous),
+      hitY(0),hitZ(0),hitSect(-1) {if (isDead) effL=-1;}
+    Float_t  x;
+    Float_t  x2x0;
+    Float_t  rphiRes;
+    Float_t  zRes;
+    Float_t  eff;
+    Bool_t   isDead;
+    //
+    Float_t  hitY;
+    Float_t  hitZ;
+    Float_t  hitSect;
+  };
+  typedef struct FT2TPCLayer FT2TPCLayer_t;
+       
+       
+  FT2();
+  virtual ~FT2();
+       
+  void InitEnvLocal();
+  void InitTPCSignalFiles(const char *TPCsignalElectrons,const char *TPCsignalMuons,const char *TPCsignalPions,const char *TPCsignalKaons,const char *TPCsignalProtons);
+  void InitTPCPIDResponse();
+  void InitDetector(Bool_t addTPC=kTRUE, Float_t sigYTPC=0.1, Float_t sigZTPC=0.1,
+                   Float_t effTPC=0.99, Float_t scEdge=2.6);
+  //
+  void PrintLayout();
+  //
+  Bool_t ProcessTrack(TParticle* trc, AliVertex* vtx);
+  void   SetSimMat(Bool_t v=kTRUE) {fSimMat = v;}
+  void   SetMinTPCHits(int v=60)  {fMinTPCHits=v;}
+  void   SetMinITSLrHit(int v=4)  {fMinITSLrHit=v;}
+  void   SetUsePIDForTracking(Bool_t usePID) {fUsePIDForTracking=usePID;}
+  //
+  Int_t    GetNClITSFakes()  const {return fNClITSFakes;}
+  Int_t    GetNClITS()  const {return fNClITS;}
+  Int_t    GetNClTPC()  const {return fNClTPC;}
+  Double_t GetChi2ITS() const {return fChi2ITS;}
+  Double_t GetChi2TPC() const {return fChi2TPC;}
+  FTProbe& GetProbe() const {return (FTProbe&)fProbe;}
+  AliExternalTrackParam& GetKalmanOut(int i) {return (AliExternalTrackParam&)fKalmanOutward[i];}
+  void   SetUseKalmanOut(Bool_t v=kTRUE)  {fUseKalmanOut = v;}
+  Bool_t GetUseKalmanOut()   const {return fUseKalmanOut;}
+  //
+  const Double_t* GetDCA()    const {return &fDCA[0];}
+  const Double_t* GetDCACov() const {return &fDCACov[0];}
+  static float GetMaxStepTGeo() {return fgMaxStepTGeo;}
+  static void  SetMaxStepTGeo(float stp=1.) {fgMaxStepTGeo = stp;}
+  //
+  void     SetdNdY(double v=-1) {fdNdY = v;}
+  Double_t GetdNdY() const {return fdNdY;}
+  Double_t HitDensity(double r2, double tgl) const;
+  Bool_t   BiasAsFake(double yz[2], const double* extyz, const double *cov) const;
+  Bool_t DiagonalizeErrors(const double *cov, double &sy2d, double &sz2d);
+  Int_t  GetITSPattern() const {return fITSPattern;}
+  Int_t  GetITSPatternFakes() const {return fITSPatternFake;}
+  //
+ protected:
+  void AddTPC(Float_t sigY=0.1, Float_t sigZ=0.1, Float_t eff=0.99, Float_t scEdge=2.6);
+  void AddTPCLayer(Float_t x, Float_t x2x0,Float_t sigY, Float_t sigZ, Float_t eff);
+  Bool_t InitProbe(TParticle* trc);
+  Bool_t PrepareProbe();
+  Bool_t ApplyMSEloss(double x2X0, double xrho);
+  Bool_t PropagateToX(double xTgt, int dir,Bool_t propErr,Bool_t simMat,Bool_t useTGeo);
+  Bool_t PropagateToR(double xTgt, int dir,Bool_t propErr,Bool_t simMat,Bool_t useTGeo);
+  Bool_t IsZero(double val, double tol=1e-9) const {return TMath::Abs(val)<tol;}
+  Bool_t PassActiveITSLayer(AliITSURecoLayer* lr);
+  Bool_t GetRoadWidth(AliITSURecoLayer* lrA,double *pr);
+  void   ResetCovMat(AliExternalTrackParam* trc);
+  Double_t UpdateKalman(AliExternalTrackParam* trc, double y,double z,double sigY,double sigZ,Bool_t randomize=kTRUE,Bool_t fake=kFALSE);
+  Bool_t ReconstructProbe();
+  AliPIDResponse::EDetPidStatus GetComputeTPCProbability (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
+  //
+ protected:
+  AliITSUReconstructor* fITSRec; // interface for reconstructor
+  AliITSURecoDet* fITS;          // interface for ITS
+  Bool_t fUsePIDForTracking;
+  Bool_t fIsTPC;                 // TPC added
+  AliPIDResponse* fPIDResponse;
+  TFile *fTPCSignalElectron;
+  TFile *fTPCSignalMuon;
+  TFile *fTPCSignalPion;
+  TFile *fTPCSignalKaon;
+  TFile *fTPCSignalProton;
+  Float_t fTPCSectorEdge;        // cut in cm on sector edge
+  Double_t fMaxSnpTPC;           // stop particle if snp>fMaxSnpTPC
+  std::vector<FT2TPCLayer_t> fTPCLayers;
+  std::vector<int>           fTPCHitLr;
+  Int_t                      fNTPCHits; //!
+  Int_t                      fMinTPCHits; // require min amount of TPC hits
+  Int_t                      fMinITSLrHit; // require min amount of ITS lr hit
+  Double_t fDCA[2],fDCACov[3];   //! dca to vertex and its covariance
+  //
+  FTProbe fProbe;  // track
+  AliExternalTrackParam* fKalmanOutward; //! parameters of outward kalman 
+  Bool_t  fUseKalmanOut;                 //! use KalmanOut estimate for fakes
+  //Double_t              fProbeMass; // probe mass
+  Double_t              fBz;        // bz
+  Bool_t                fSimMat;    // simulate material effects in probe preparation
+  //
+  Int_t                 fCurrITSLr; //! current ITS layer under tracking
+  Int_t                 fNClTPC;    //! N used TPC clusters
+  Int_t                 fNClITS;    //! N used ITS clusters
+  Int_t                 fNClITSFakes;    //! N used ITS Fake clusters
+  Int_t                 fITSPatternFake; //! fakes pattern for ITS
+  Int_t                 fITSPattern;     //! pattern for ITS clusters
+  Double_t              fChi2TPC;   //! total chi2 in TPC
+  Double_t              fChi2ITS;   //! total chi2 in ITS
+  //
+  // hit info in the ITS
+  Double_t fSigYITS,fSigZITS;       // nominal ITS later resolution
+  Int_t fNITSHits[kMaxITSLr]; //! n hits per ITS layer
+  Int_t fNITSLrHit;           //! n of ITS layers whith hit
+  AliITSURecoSens* fITSSensHit[kMaxITSLr][2]; //! hit sensors
+  Double_t fITSHitYZ[kMaxITSLr][kMaxHitPerLr][2]; //! tracking Y,Z of each hit
+  //
+  Double_t fdNdY;             //! if positive, use it for fakes simulation
+  //
+  static float fgMaxStepTGeo; // max step for tracking accounting for TGeo materials
+       
+  ClassDef(FT2,1)
+};
+
+#endif
diff --git a/ITS/UPGRADE/FT2/README b/ITS/UPGRADE/FT2/README
new file mode 100644 (file)
index 0000000..1c77dac
--- /dev/null
@@ -0,0 +1,17 @@
+Don Dez 11 20:22:25 CET 2014
+
+FT2: semi-analitical fast MC
+
+See testF.C macro for example of running.
+With dndy>=0 fakes will be generated (if useKalmanOut==kTRUE - 
+with quasi-smoothing procedure)
+
+Before running the test, one should have in the working directory
+the files from fullMC:
+geometry.root
+itsSegmentations.root
+ITS // with ITSU specific stuff
+GRP // grp created by MC
+
+
+To be completed ...
diff --git a/ITS/UPGRADE/FT2/test.C b/ITS/UPGRADE/FT2/test.C
new file mode 100644 (file)
index 0000000..2f85c52
--- /dev/null
@@ -0,0 +1,48 @@
+TH2* h=0;
+TH1F *hdcar=0,*hdcaz=0;
+
+void test()
+{
+  gSystem->Load("libITSUpgradeBase.so");
+  gSystem->Load("libITSUpgradeSim.so");
+  gSystem->Load("libITSUpgradeRec.so");
+  gROOT->ProcessLine(".L FT2.cxx+");
+  //
+  TParticle prt;
+  //
+  FT2* det = new FT2();
+  det->InitEnvLocal();
+  det->InitDetector();
+  det->SetSimMat(kTRUE); // simulate mat.effects in probe preparation
+  //
+  det->SetMinTPCHits(60);
+  det->SetMinITSLrHit(4);
+  //
+  AliESDVertex* vtx = new AliESDVertex();
+
+  h = new TH2F("h","h",100,0,TMath::Pi()*2,100,1,-1);
+  hdcar = new TH1F("hdcar","dcaXY",100,-100e-4,100e-4);
+  hdcaz = new TH1F("hdcaz","dcaZ",100,-100e-4,100e-4);
+
+  double pt = 1;//0.2;
+  for (int i=0;i<50000;i++) {
+    double phi = gRandom->Rndm()*TMath::Pi()*2;
+    double eta = gRandom->Rndm()*0.8;
+    double theta = 2*TMath::ATan(TMath::Exp(-eta));
+    double pz = pt/TMath::Tan(theta);
+    double px = pt*TMath::Cos(phi);
+    double py = pt*TMath::Sin(phi);
+    //    double pxyz[3]={0.2,0.7,0.4};
+    double en = TMath::Sqrt(pt*pt+pz*pz+0.14*0.14);
+    prt.SetPdgCode(211);
+    prt.SetMomentum(px,py,pz,en);
+    if (det->ProcessTrack(&prt,vtx)) {
+      AliExternalTrackParam& prb = det->GetProbe();
+      h->Fill(prb.Phi(),prb.GetSigmaY2());
+      const double* dca = det->GetDCA();
+      hdcar->Fill(dca[0]);
+      hdcaz->Fill(dca[1]);
+    }
+  }
+
+}
diff --git a/ITS/UPGRADE/FT2/testF.C b/ITS/UPGRADE/FT2/testF.C
new file mode 100644 (file)
index 0000000..12ee7ac
--- /dev/null
@@ -0,0 +1,108 @@
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TH1F.h>
+#include <TSystem.h>
+#include <TMath.h>
+#include <TParticle.h>
+#include <TRandom.h>
+#include <TROOT.h>
+#include "AliESDVertex.h"
+#include "AliExternalTrackParam.h"
+#include "FT2.h"
+
+FT2* det=0;
+#endif
+
+TH1F* hdca=0,*hdcaN=0;
+TH1F* hdcaZ=0,*hdcaZN=0;
+TH1F* hFdca=0,*hFdcaN=0;
+TH1F* hFdcaZ=0,*hFdcaZN=0;
+TH1F* hPatternITS=0;
+TH1F* hFPatternITS=0;
+
+void testF(int ntrials=10000, double dndy=2100., Bool_t useKalmanOut=kTRUE)
+{
+#if defined(__CINT__) && !defined(__MAKECINT__)
+  gSystem->Load("libITSUpgradeBase.so");
+  gSystem->Load("libITSUpgradeSim.so");
+  gSystem->Load("libITSUpgradeRec.so");
+  gROOT->ProcessLine(".L FT2.cxx+");
+  FT2* det=0;
+#endif
+  //
+  det = new FT2();
+  det->InitEnvLocal();
+  det->InitDetector();
+  det->SetSimMat(kTRUE);
+  det->SetMaxStepTGeo(1.);
+  det->SetdNdY(dndy);
+  det->SetUseKalmanOut(useKalmanOut);
+  //
+  hdca  = new TH1F("hdca","dca",  100,-0.1,.1);
+  hdcaN = new TH1F("hdcaN","dcaN",100,-10.,10.);
+  hdcaZ  = new TH1F("hdcaZ","dcaZ",  100,-0.1,.1);
+  hdcaZN = new TH1F("hdcaZN","dcaZN",100,-10.,10.);
+  //
+  hFdca  = new TH1F("hFdca","dca fake",  100,-0.1,.1);
+  hFdcaN = new TH1F("hFdcaN","dcaN fake",100,-10.,10.);
+  hFdcaZ  = new TH1F("hFdcaZ","dcaZ fake",  100,-0.1,.1);
+  hFdcaZN = new TH1F("hFdcaZN","dcaZN fake",100,-10.,10.);
+  //
+  hPatternITS  = new TH1F("itsPattern","ITS hits pattern"      ,7,-0.5,6.5);
+  hFPatternITS = new TH1F("itsFPattern","ITS fake hits pattern",7,-0.5,6.5);
+  //
+  AliESDVertex *vtx = new AliESDVertex();
+  double vcov[6] = {1e-4, 0, 1e-4, 0, 0, 1e-4};
+  vtx->SetCovarianceMatrix(vcov);
+  vtx->Print();
+  //
+  TParticle prt;
+  double p = 0.25;
+  for (int ntr=0;ntr<ntrials;ntr++) {
+    
+    vtx->SetXv(gRandom->Gaus(-0.1, 50e-4));
+    vtx->SetYv(gRandom->Gaus(0.2,  50e-4));
+    vtx->SetZv(gRandom->Gaus(0.5, 5.0));
+    //
+    double phi = gRandom->Rndm()*TMath::TwoPi();
+    double eta =  2*(gRandom->Rndm()-0.5)*0.8;
+    double theta = 2*TMath::ATan(TMath::Exp(-eta));
+    double pz = p*TMath::Cos(theta);
+    double pt = p*TMath::Sin(theta);
+    double pxyz[3]={pt*TMath::Cos(phi),pt*TMath::Sin(phi),pz};
+    double en = TMath::Sqrt(p*p+0.14*0.14);
+    prt.SetPdgCode(gRandom->Rndm()>0.5 ? 211 : -211);
+    //prt.SetPdgCode(-211);
+    prt.SetMomentum(pxyz[0],pxyz[1],pxyz[2],en);
+    prt.SetProductionVertex(vtx->GetX(),vtx->GetY(),vtx->GetZ(),0);
+    //
+    if (det->ProcessTrack(&prt,vtx)) {
+      //const AliExternalTrackParam& prob = det->GetProbe();
+      //      printf("%d %d\n",det->GetNClITS(),det->GetNClTPC());
+      //    prob.Print();
+      //      /*
+      const double* dca = det->GetDCA();
+      const double* cov = det->GetDCACov();
+      hdca->Fill(dca[0]);
+      hdcaN->Fill(dca[0]/TMath::Sqrt(cov[0]));
+      hdcaZ->Fill(dca[1]);
+      hdcaZN->Fill(dca[1]/TMath::Sqrt(cov[2]));
+      //
+      if (det->GetNClITSFakes()) {
+       hFdca->Fill(dca[0]);
+       hFdcaN->Fill(dca[0]/TMath::Sqrt(cov[0]));
+       hFdcaZ->Fill(dca[1]);
+       hFdcaZN->Fill(dca[1]/TMath::Sqrt(cov[2]));
+      }
+      int hits  = det->GetITSPattern();
+      int hitsF = det->GetITSPatternFakes();
+      for (int j=7;j--;) {
+       if (hits&(0x1<<j))  hPatternITS->Fill(j);
+       if (hitsF&(0x1<<j)) hFPatternITS->Fill(j);      
+      }
+    }
+    else {
+      printf("Failed on track %d\n",ntr);
+    }
+  }
+  //
+}