[u/mrichter/AliRoot.git] / TPC / AliTPCcalibCalib.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


///////////////////////////////////////////////////////////////////////////////
//                                                                           //
//   Component for redoing the reconstruction from the clusters and tracks
//      
//   The new calibration data used 
//
//   In reality it overwrites the content of the ESD 
//    

/*

  gSystem->Load("libANALYSIS");
  gSystem->Load("libTPCcalib");
  //
  gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
  gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
  AliXRDPROOFtoolkit tool; 
  TChain * chainCl = tool.MakeChain("calib.txt","Clusters",0,1);
  chainCl->Lookup();
  TChain * chainTr = tool.MakeChain("calib.txt","Tracks",0,1);
  chainTr->Lookup();


*/


//  marian.ivanov@cern.ch
// 
#include "AliTPCcalibCalib.h"
#include "TSystem.h"
#include "TFile.h"
#include "TTreeStream.h"
#include "AliLog.h"
#include "TTimeStamp.h"
#include "AliESDEvent.h"
#include "AliESDfriend.h"
#include "AliESDtrack.h"
#include "AliTracker.h"
#include "AliTPCClusterParam.h"

#include "AliTPCcalibDB.h"
#include "AliTPCTransform.h"
#include "AliTPCRecoParam.h"
#include "AliTPCclusterMI.h"
#include "AliTPCseed.h"
#include "AliTPCPointCorrection.h"

ClassImp(AliTPCcalibCalib)

AliTPCcalibCalib::AliTPCcalibCalib():
AliTPCcalibBase(),
  fApplyExBCorrection(1),      // apply ExB correction
  fApplyTOFCorrection(1),      // apply TOF correction
  fApplyPositionCorrection(1), // apply position correction
  fApplySectorAlignment(1),    // apply sector alignment
  fApplyRPhiCorrection(1),     // apply R-Phi correction
  fApplyRCorrection(1)         // apply Radial correction

{
  //
  // Constructor
  //
}


AliTPCcalibCalib::AliTPCcalibCalib(const Text_t *name, const Text_t *title) 
  :AliTPCcalibBase(),
  fApplyExBCorrection(1),      // apply ExB correction
  fApplyTOFCorrection(1),      // apply TOF correction
  fApplyPositionCorrection(1), // apply position correction
  fApplySectorAlignment(1),    // apply sector alignment
  fApplyRPhiCorrection(1),     // apply R-Phi correction
  fApplyRCorrection(1)         // apply Radial correction
{  
  SetName(name);
  SetTitle(title);
}


AliTPCcalibCalib::AliTPCcalibCalib(const AliTPCcalibCalib&calib):
  AliTPCcalibBase(calib),
  fApplyExBCorrection(calib.GetApplyExBCorrection()),
  fApplyTOFCorrection(calib.GetApplyTOFCorrection()),
  fApplyPositionCorrection(calib.GetApplyPositionCorrection()),
  fApplySectorAlignment(calib.GetApplySectorAlignment()),
  fApplyRPhiCorrection(calib.GetApplyRPhiCorrection()),
  fApplyRCorrection(calib.GetApplyRCorrection())

{
  //
  // copy constructor
  //
}

AliTPCcalibCalib &AliTPCcalibCalib::operator=(const AliTPCcalibCalib&calib){
  //
  //
  //
  ((AliTPCcalibBase *)this)->operator=(calib);
  return *this;
}


AliTPCcalibCalib::~AliTPCcalibCalib() {
  //
  // destructor
  //
}


void     AliTPCcalibCalib::Process(AliESDEvent *event){
  //
  // 
  //
  if (!event) {
    return;
  }  
  AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
  if (!ESDfriend) {
   return;
  }
  
  if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
  Int_t ntracks=event->GetNumberOfTracks();   
  //AliTPCcalibDB::Instance()->SetExBField(fMagF);

  //
  //
  //

  for (Int_t i=0;i<ntracks;++i) {
    AliESDtrack *track = event->GetTrack(i);     
    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
 
    const AliExternalTrackParam * trackIn  = track->GetInnerParam();
    const AliExternalTrackParam * trackOut = track->GetOuterParam();
    AliExternalTrackParam * tpcOut   = (AliExternalTrackParam *)friendTrack->GetTPCOut();
    if (!trackIn) continue;
    if (!trackOut) continue;
    if (!tpcOut) continue;   
    TObject *calibObject;
    AliTPCseed *seed = 0;
    for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
      if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
    }
    if (!seed) continue;
    RefitTrack(track, seed,event->GetMagneticField());
    (*tpcOut)=*(track->GetOuterParam());  
  }
  return;
}

Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Float_t magesd){
  //
  // Refit track
  // if magesd==0 forget the curvature

  //
  // 0 - Setup transform object
  //
  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
  transform->SetCurrentRun(fRun);
  transform->SetCurrentTimeStamp((UInt_t)fTime);
  if(!fApplyExBCorrection) { // disable ExB correction in transform
    if(transform->GetCurrentRecoParam())
      transform->GetCurrentRecoParamNonConst()->SetUseExBCorrection(0);
  }
  if(!fApplyTOFCorrection) { // disable TOF correction in transform
    if(transform->GetCurrentRecoParam())
      transform->GetCurrentRecoParamNonConst()->SetUseTOFCorrection(kFALSE);
  }

  //
  // First apply calibration
  //
  AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
  for (Int_t irow=0;irow<159;irow++) {
    AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
    if (!cluster) continue; 
    AliTPCclusterMI cl0(*cluster);
    Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
    Int_t i[1]={cluster->GetDetector()};

    transform->Transform(x,i,0,1);
    //
    // get position correction
    //
    Int_t ipad=0;
    if (cluster->GetDetector()>35) ipad=1;
    Float_t dy =AliTPCClusterParam::SPosCorrection(0,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
    Float_t dz =AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
    // if(fApplyPositionCorrection) {
    //x[1]-=dy;
    //x[2]-=dz;
    // }
    //
    // Apply sector alignment
    //
    Double_t dxq = AliTPCPointCorrection::SGetCorrectionSector(0,cluster->GetDetector()%36,cluster->GetX(),
							       cluster->GetY(),cluster->GetZ());
    Double_t dyq = AliTPCPointCorrection::SGetCorrectionSector(1,cluster->GetDetector()%36,cluster->GetX(),
							       cluster->GetY(),cluster->GetZ());
    Double_t dzq = AliTPCPointCorrection::SGetCorrectionSector(2,cluster->GetDetector()%36,cluster->GetX(),
							       cluster->GetY(),cluster->GetZ());
    if (fApplySectorAlignment){
      x[0]-=dxq;
      x[1]-=dyq;
      x[2]-=dzq;
    }
    //
    // Apply r-phi correction  - To be done on track level- knowing the track angle !!!
    //
    Double_t corrclY =  
      corr->RPhiCOGCorrection(cluster->GetDetector(),cluster->GetRow(), cluster->GetPad(),
				  cluster->GetY(),cluster->GetY(), cluster->GetZ(), 0., cluster->GetMax(),2.5);
    // R correction
    Double_t corrR   = corr->CorrectionOutR0(kFALSE,kFALSE,cluster->GetX(),cluster->GetY(),cluster->GetZ(),cluster->GetDetector());

    if (fApplyRPhiCorrection){
      if (cluster->GetY()>0) x[1]+=corrclY;  // rphi correction
      if (cluster->GetY()<0) x[1]-=corrclY;  // rphi correction
    }
    if (fApplyRCorrection){      
      x[0]+=corrR;                           // radial correction
    }

    //
    //
    //
    cluster->SetX(x[0]);
    cluster->SetY(x[1]);
    cluster->SetZ(x[2]);
    if (fStreamLevel>2){
      TTreeSRedirector *cstream = GetDebugStreamer();
      if (cstream){
	(*cstream)<<"Clusters"<<
	  "run="<<fRun<<              //  run number
	  "event="<<fEvent<<          //  event number
	  "time="<<fTime<<            //  time stamp of event
	  "trigger="<<fTrigger<<      //  trigger
	  "triggerClass="<<&fTriggerClass<<      //  trigger	  
	  "mag="<<fMagF<<             //  magnetic field
	  "cl0.="<<&cl0<<
	  "cl.="<<cluster<<
	  "cy="<<dy<<
	  "cz="<<dz<<
	  "cY="<<corrclY<<
	  "cR="<<corrR<<
	  "dxq="<<dxq<<
	  "dyq="<<dyq<<
	  "dzq="<<dzq<<
	  "\n";
      }
    }
  }
  //
  //
  //
  Int_t ncl = seed->GetNumberOfClusters();
  Double_t covar[15];
  for (Int_t i=0;i<15;i++) covar[i]=0;
  covar[0]=10.*10.;
  covar[2]=10.*10.;
  covar[5]=10.*10./(64.*64.);
  covar[9]=10.*10./(64.*64.);
  covar[14]=0.3*0.3;
  if (TMath::Abs(magesd)<0.05) {
     covar[14]=0.025*0.025;
  }
  // 
  // And now do refit
  //
  AliExternalTrackParam * trackInOld  = (AliExternalTrackParam*)track->GetInnerParam();
  AliExternalTrackParam * trackOutOld = (AliExternalTrackParam*)track->GetOuterParam();


  AliExternalTrackParam trackIn  = *trackOutOld;
  if (TMath::Abs(magesd)<0.05) {
    ((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
  }  
  trackIn.ResetCovariance(20.);
  trackIn.AddCovariance(covar);
  Double_t xyz[3];
  Int_t nclIn=0,nclOut=0;
  //
  // Refit in
  //

  for (Int_t irow=159; irow>0; irow--){
    AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
    if (!cl) continue;
    if (cl->GetX()<80) continue;
    Int_t sector = cl->GetDetector();
    Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
    if (TMath::Abs(dalpha)>0.01){
      if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
    }
    Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
    Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
    AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
    cov[0]*=cov[0];
    cov[2]*=cov[2];
    trackIn.GetXYZ(xyz);
    Double_t bz = AliTracker::GetBz(xyz);

    if (!trackIn.PropagateTo(r[0],bz)) continue;
    if (RejectCluster(cl,&trackIn)) continue;
    nclIn++;
    trackIn.Update(&r[1],cov);    
  }
  //
  AliExternalTrackParam trackOut = trackIn;
  if (TMath::Abs(magesd)<0.05) {
    ((Double_t&)(trackOut.GetParameter()[4]))=0.000000001;
  }
  trackOut.ResetCovariance(20.);
  trackOut.AddCovariance(covar);
  //
  // Refit out
  //
  //Bool_t lastEdge=kFALSE;
  for (Int_t irow=0; irow<160; irow++){
    AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
    if (!cl) continue;
    if (cl->GetX()<80) continue;
    Int_t sector = cl->GetDetector();
    Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();

    if (TMath::Abs(dalpha)>0.01){
      if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
    }
    Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};

    Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation    
    AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
    cov[0]*=cov[0];
    cov[2]*=cov[2];
    trackOut.GetXYZ(xyz);
    Double_t bz = AliTracker::GetBz(xyz);
    if (!trackOut.PropagateTo(r[0],bz)) continue;
    if (RejectCluster(cl,&trackOut)) continue;
    nclOut++;
    trackOut.Update(&r[1],cov);    	
    //if (cl->GetType()<0) lastEdge=kTRUE;
    //if (cl->GetType()>=0) lastEdge=kFALSE;    
  }
  //
  //
  //
  nclIn=0;
  trackIn  = trackOut;
  trackIn.ResetCovariance(10.);
  //
  // Refit in one more time
  //
  for (Int_t irow=159; irow>0; irow--){
    AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
    if (!cl) continue;
    if (cl->GetX()<80) continue;
    Int_t sector = cl->GetDetector();
    Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
    if (TMath::Abs(dalpha)>0.01){
      if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
    }
    Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
    Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
    AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
    cov[0]*=cov[0];
    cov[2]*=cov[2];
    trackIn.GetXYZ(xyz);
    Double_t bz = AliTracker::GetBz(xyz);

    if (!trackIn.PropagateTo(r[0],bz)) continue;
    if (RejectCluster(cl,&trackIn)) continue;
    nclIn++;
    trackIn.Update(&r[1],cov);    
  }


  trackIn.Rotate(trackInOld->GetAlpha());
  trackOut.Rotate(trackOutOld->GetAlpha());
  //
  trackInOld->GetXYZ(xyz);
  Double_t bz = AliTracker::GetBz(xyz);
  trackIn.PropagateTo(trackInOld->GetX(),bz);
  //
  trackOutOld->GetXYZ(xyz);
  bz = AliTracker::GetBz(xyz);  
  trackOut.PropagateTo(trackOutOld->GetX(),bz);
  

  if (fStreamLevel>0){
    TTreeSRedirector *cstream = GetDebugStreamer();
    if (cstream){
      (*cstream)<<"Tracks"<<
	"run="<<fRun<<              //  run number
	"event="<<fEvent<<          //  event number
	"time="<<fTime<<            //  time stamp of event
	"trigger="<<fTrigger<<      //  trigger
	"triggerClass="<<&fTriggerClass<<      //  trigger
	"mag="<<fMagF<<             //  magnetic field
	"nclIn="<<nclIn<<
	"nclOut="<<nclOut<<
	"ncl="<<ncl<<
	"TrIn0.="<<trackInOld<<
	"TrOut0.="<<trackOutOld<<
	"TrIn1.="<<&trackIn<<
	"TrOut1.="<<&trackOut<<
	"\n";
    }
  }
  //
  // And now rewrite ESDtrack and TPC seed
  //
 
  (*trackInOld)  = trackIn;
  (*trackOutOld) = trackOut;
  AliExternalTrackParam *t = &trackIn;
  track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
  seed->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
  seed->SetNumberOfClusters((nclIn+nclOut)/2);
  return kTRUE;
}


Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
  //
  // check the acceptance of cluster
  // Cut on edge effects
  //
  Float_t kEdgeCut=2.5;
  Float_t kSigmaCut=6;

  Bool_t isReject = kFALSE;
  Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
  Float_t dist  = edgeY - TMath::Abs(cl->GetY());
  if (param)  dist  = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
  if (dist<kEdgeCut) isReject=kTRUE;

  Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation    
  AliTPCseed::GetError(cl, param,cov[0],cov[2]);
  if (param->GetSigmaY2()<0 || param->GetSigmaZ2()<0){
    AliError("Wrong parameters");
    return kFALSE;
  }
  Double_t py = (cl->GetY()-param->GetY())/TMath::Sqrt(cov[0]*cov[0]+param->GetSigmaY2());
  Double_t pz = (cl->GetZ()-param->GetZ())/TMath::Sqrt(cov[2]*cov[2]+param->GetSigmaZ2());
  //
  if ((py*py+pz*pz)>kSigmaCut*kSigmaCut) isReject=kTRUE;
  
  return isReject;
}
Commit	Line	Data
9dcfce73	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16
	17	///////////////////////////////////////////////////////////////////////////////
	18	// //
	19	// Component for redoing the reconstruction from the clusters and tracks
	20	//
	21	// The new calibration data used
	22	//
	23	// In reality it overwrites the content of the ESD
	24	//
	25
3bf5c7a6	26	/*
	27
	28	gSystem->Load("libANALYSIS");
	29	gSystem->Load("libTPCcalib");
	30	//
	31	gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
	32	gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
	33	AliXRDPROOFtoolkit tool;
11c7a854	34	TChain * chainCl = tool.MakeChain("calib.txt","Clusters",0,1);
	35	chainCl->Lookup();
	36	TChain * chainTr = tool.MakeChain("calib.txt","Tracks",0,1);
	37	chainTr->Lookup();
3bf5c7a6	38
	39
	40
	41	*/
	42
	43
	44
9dcfce73	45	// marian.ivanov@cern.ch
	46	//
	47	#include "AliTPCcalibCalib.h"
	48	#include "TSystem.h"
	49	#include "TFile.h"
	50	#include "TTreeStream.h"
	51	#include "AliLog.h"
	52	#include "TTimeStamp.h"
	53	#include "AliESDEvent.h"
	54	#include "AliESDfriend.h"
	55	#include "AliESDtrack.h"
	56	#include "AliTracker.h"
11c7a854	57	#include "AliTPCClusterParam.h"
9dcfce73	58
	59	#include "AliTPCcalibDB.h"
	60	#include "AliTPCTransform.h"
7af539c6	61	#include "AliTPCRecoParam.h"
9dcfce73	62	#include "AliTPCclusterMI.h"
9dcfce73	63	#include "AliTPCseed.h"
15e48021	64	#include "AliTPCPointCorrection.h"
9dcfce73	65
	66	ClassImp(AliTPCcalibCalib)
	67
	68	AliTPCcalibCalib::AliTPCcalibCalib():
7af539c6	69	AliTPCcalibBase(),
	70	fApplyExBCorrection(1), // apply ExB correction
	71	fApplyTOFCorrection(1), // apply TOF correction
	72	fApplyPositionCorrection(1), // apply position correction
	73	fApplySectorAlignment(1), // apply sector alignment
	74	fApplyRPhiCorrection(1), // apply R-Phi correction
	75	fApplyRCorrection(1) // apply Radial correction
	76
9dcfce73	77	{
	78	//
	79	// Constructor
	80	//
	81	}
	82
	83
	84	AliTPCcalibCalib::AliTPCcalibCalib(const Text_t name, const Text_t title)
7af539c6	85	:AliTPCcalibBase(),
	86	fApplyExBCorrection(1), // apply ExB correction
	87	fApplyTOFCorrection(1), // apply TOF correction
	88	fApplyPositionCorrection(1), // apply position correction
	89	fApplySectorAlignment(1), // apply sector alignment
	90	fApplyRPhiCorrection(1), // apply R-Phi correction
	91	fApplyRCorrection(1) // apply Radial correction
9dcfce73	92	{
	93	SetName(name);
	94	SetTitle(title);
	95	}
	96
	97
	98	AliTPCcalibCalib::AliTPCcalibCalib(const AliTPCcalibCalib&calib):
7af539c6	99	AliTPCcalibBase(calib),
	100	fApplyExBCorrection(calib.GetApplyExBCorrection()),
	101	fApplyTOFCorrection(calib.GetApplyTOFCorrection()),
	102	fApplyPositionCorrection(calib.GetApplyPositionCorrection()),
	103	fApplySectorAlignment(calib.GetApplySectorAlignment()),
	104	fApplyRPhiCorrection(calib.GetApplyRPhiCorrection()),
	105	fApplyRCorrection(calib.GetApplyRCorrection())
	106
9dcfce73	107	{
	108	//
	109	// copy constructor
	110	//
	111	}
	112
	113	AliTPCcalibCalib &AliTPCcalibCalib::operator=(const AliTPCcalibCalib&calib){
	114	//
	115	//
	116	//
	117	((AliTPCcalibBase *)this)->operator=(calib);
	118	return *this;
	119	}
	120
	121
	122	AliTPCcalibCalib::~AliTPCcalibCalib() {
	123	//
	124	// destructor
	125	//
	126	}
	127
	128
	129	void AliTPCcalibCalib::Process(AliESDEvent *event){
	130	//
	131	//
	132	//
	133	if (!event) {
	134	return;
	135	}
	136	AliESDfriend ESDfriend=static_cast<AliESDfriend>(event->FindListObject("AliESDfriend"));
	137	if (!ESDfriend) {
	138	return;
	139	}
	140
	141	if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
	142	Int_t ntracks=event->GetNumberOfTracks();
91350f5b	143	//AliTPCcalibDB::Instance()->SetExBField(fMagF);
5c3e0d17	144
9dcfce73	145	//
	146	//
	147	//
	148
	149	for (Int_t i=0;i<ntracks;++i) {
3e55050f	150	AliESDtrack *track = event->GetTrack(i);
	151	AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
	152
	153	const AliExternalTrackParam * trackIn = track->GetInnerParam();
9dcfce73	154	const AliExternalTrackParam * trackOut = track->GetOuterParam();
3e55050f	155	AliExternalTrackParam * tpcOut = (AliExternalTrackParam *)friendTrack->GetTPCOut();
9dcfce73	156	if (!trackIn) continue;
9dcfce73	157	if (!trackOut) continue;
3e55050f	158	if (!tpcOut) continue;
9dcfce73	159	TObject *calibObject;
	160	AliTPCseed *seed = 0;
	161	for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
	162	if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
	163	}
	164	if (!seed) continue;
38b1a1ca	165	RefitTrack(track, seed,event->GetMagneticField());
3e55050f	166	(tpcOut)=(track->GetOuterParam());
9dcfce73	167	}
	168	return;
	169	}
	170
38b1a1ca	171	Bool_t AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Float_t magesd){
9dcfce73	172	//
9dcfce73	173	// Refit track
38b1a1ca	174	// if magesd==0 forget the curvature
9dcfce73	175
15e48021	176	//
	177	// 0 - Setup transform object
	178	//
	179	AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
	180	transform->SetCurrentRun(fRun);
	181	transform->SetCurrentTimeStamp((UInt_t)fTime);
7af539c6	182	if(!fApplyExBCorrection) { // disable ExB correction in transform
	183	if(transform->GetCurrentRecoParam())
	184	transform->GetCurrentRecoParamNonConst()->SetUseExBCorrection(0);
	185	}
	186	if(!fApplyTOFCorrection) { // disable TOF correction in transform
	187	if(transform->GetCurrentRecoParam())
	188	transform->GetCurrentRecoParamNonConst()->SetUseTOFCorrection(kFALSE);
	189	}
	190
9dcfce73	191	//
	192	// First apply calibration
	193	//
15e48021	194	AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
9dcfce73	195	for (Int_t irow=0;irow<159;irow++) {
	196	AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
	197	if (!cluster) continue;
	198	AliTPCclusterMI cl0(*cluster);
	199	Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
	200	Int_t i[1]={cluster->GetDetector()};
7af539c6	201
9dcfce73	202	transform->Transform(x,i,0,1);
11c7a854	203	//
	204	// get position correction
	205	//
	206	Int_t ipad=0;
	207	if (cluster->GetDetector()>35) ipad=1;
	208	Float_t dy =AliTPCClusterParam::SPosCorrection(0,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
	209	Float_t dz =AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
7af539c6	210	// if(fApplyPositionCorrection) {
11c7a854	211	//x[1]-=dy;
11c7a854	212	//x[2]-=dz;
7af539c6	213	// }
15e48021	214	//
	215	// Apply sector alignment
	216	//
	217	Double_t dxq = AliTPCPointCorrection::SGetCorrectionSector(0,cluster->GetDetector()%36,cluster->GetX(),
	218	cluster->GetY(),cluster->GetZ());
	219	Double_t dyq = AliTPCPointCorrection::SGetCorrectionSector(1,cluster->GetDetector()%36,cluster->GetX(),
	220	cluster->GetY(),cluster->GetZ());
	221	Double_t dzq = AliTPCPointCorrection::SGetCorrectionSector(2,cluster->GetDetector()%36,cluster->GetX(),
	222	cluster->GetY(),cluster->GetZ());
7af539c6	223	if (fApplySectorAlignment){
15e48021	224	x[0]-=dxq;
	225	x[1]-=dyq;
	226	x[2]-=dzq;
	227	}
	228	//
	229	// Apply r-phi correction - To be done on track level- knowing the track angle !!!
	230	//
	231	Double_t corrclY =
	232	corr->RPhiCOGCorrection(cluster->GetDetector(),cluster->GetRow(), cluster->GetPad(),
	233	cluster->GetY(),cluster->GetY(), cluster->GetZ(), 0., cluster->GetMax(),2.5);
	234	// R correction
	235	Double_t corrR = corr->CorrectionOutR0(kFALSE,kFALSE,cluster->GetX(),cluster->GetY(),cluster->GetZ(),cluster->GetDetector());
	236
7af539c6	237	if (fApplyRPhiCorrection){
15e48021	238	if (cluster->GetY()>0) x[1]+=corrclY; // rphi correction
15e48021	239	if (cluster->GetY()<0) x[1]-=corrclY; // rphi correction
7af539c6	240	}
7af539c6	241	if (fApplyRCorrection){
15e48021	242	x[0]+=corrR; // radial correction
15e48021	243	}
11c7a854	244
15e48021	245	//
15e48021	246	//
9dcfce73	247	//
	248	cluster->SetX(x[0]);
	249	cluster->SetY(x[1]);
	250	cluster->SetZ(x[2]);
	251	if (fStreamLevel>2){
	252	TTreeSRedirector *cstream = GetDebugStreamer();
	253	if (cstream){
	254	(*cstream)<<"Clusters"<<
108953e9	255	"run="<<fRun<< // run number
	256	"event="<<fEvent<< // event number
	257	"time="<<fTime<< // time stamp of event
	258	"trigger="<<fTrigger<< // trigger
15e48021	259	"triggerClass="<<&fTriggerClass<< // trigger
108953e9	260	"mag="<<fMagF<< // magnetic field
9dcfce73	261	"cl0.="<<&cl0<<
9dcfce73	262	"cl.="<<cluster<<
11c7a854	263	"cy="<<dy<<
11c7a854	264	"cz="<<dz<<
fe8454d9	265	"cY="<<corrclY<<
15e48021	266	"cR="<<corrR<<
	267	"dxq="<<dxq<<
	268	"dyq="<<dyq<<
	269	"dzq="<<dzq<<
9dcfce73	270	"\n";
	271	}
	272	}
	273	}
15e48021	274	//
	275	//
	276	//
11c7a854	277	Int_t ncl = seed->GetNumberOfClusters();
	278	Double_t covar[15];
	279	for (Int_t i=0;i<15;i++) covar[i]=0;
	280	covar[0]=10.*10.;
	281	covar[2]=10.*10.;
	282	covar[5]=10.10./(64.64.);
	283	covar[9]=10.10./(64.64.);
38b1a1ca	284	covar[14]=0.3*0.3;
	285	if (TMath::Abs(magesd)<0.05) {
	286	covar[14]=0.025*0.025;
	287	}
9dcfce73	288	//
	289	// And now do refit
	290	//
11c7a854	291	AliExternalTrackParam * trackInOld = (AliExternalTrackParam*)track->GetInnerParam();
9dcfce73	292	AliExternalTrackParam * trackOutOld = (AliExternalTrackParam*)track->GetOuterParam();
9dcfce73	293
3e55050f	294
9dcfce73	295	AliExternalTrackParam trackIn = *trackOutOld;
38b1a1ca	296	if (TMath::Abs(magesd)<0.05) {
	297	((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
	298	}
	299	trackIn.ResetCovariance(20.);
11c7a854	300	trackIn.AddCovariance(covar);
9dcfce73	301	Double_t xyz[3];
	302	Int_t nclIn=0,nclOut=0;
	303	//
15e48021	304	// Refit in
	305	//
	306
	307	for (Int_t irow=159; irow>0; irow--){
	308	AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
	309	if (!cl) continue;
	310	if (cl->GetX()<80) continue;
	311	Int_t sector = cl->GetDetector();
	312	Float_t dalpha = TMath::DegToRad()(sector%1820.+10.)-trackIn.GetAlpha();
91f8fa1d	313	if (TMath::Abs(dalpha)>0.01){
	314	if (!trackIn.Rotate(TMath::DegToRad()(sector%1820.+10.))) break;
	315	}
15e48021	316	Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
	317	Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
	318	AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
	319	cov[0]*=cov[0];
	320	cov[2]*=cov[2];
	321	trackIn.GetXYZ(xyz);
	322	Double_t bz = AliTracker::GetBz(xyz);
	323
	324	if (!trackIn.PropagateTo(r[0],bz)) continue;
	325	if (RejectCluster(cl,&trackIn)) continue;
	326	nclIn++;
	327	trackIn.Update(&r[1],cov);
	328	}
	329	//
	330	AliExternalTrackParam trackOut = trackIn;
38b1a1ca	331	if (TMath::Abs(magesd)<0.05) {
	332	((Double_t&)(trackOut.GetParameter()[4]))=0.000000001;
	333	}
	334	trackOut.ResetCovariance(20.);
15e48021	335	trackOut.AddCovariance(covar);
15e48021	336	//
9dcfce73	337	// Refit out
9dcfce73	338	//
15e48021	339	//Bool_t lastEdge=kFALSE;
9dcfce73	340	for (Int_t irow=0; irow<160; irow++){
	341	AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
	342	if (!cl) continue;
	343	if (cl->GetX()<80) continue;
	344	Int_t sector = cl->GetDetector();
	345	Float_t dalpha = TMath::DegToRad()(sector%1820.+10.)-trackOut.GetAlpha();
	346
91f8fa1d	347	if (TMath::Abs(dalpha)>0.01){
	348	if (!trackOut.Rotate(TMath::DegToRad()(sector%1820.+10.))) break;
	349	}
9dcfce73	350	Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
91350f5b	351
	352	Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
	353	AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
	354	cov[0]*=cov[0];
	355	cov[2]*=cov[2];
9dcfce73	356	trackOut.GetXYZ(xyz);
9dcfce73	357	Double_t bz = AliTracker::GetBz(xyz);
15e48021	358	if (!trackOut.PropagateTo(r[0],bz)) continue;
42b40d07	359	if (RejectCluster(cl,&trackOut)) continue;
15e48021	360	nclOut++;
	361	trackOut.Update(&r[1],cov);
	362	//if (cl->GetType()<0) lastEdge=kTRUE;
	363	//if (cl->GetType()>=0) lastEdge=kFALSE;
9dcfce73	364	}
9dcfce73	365	//
9dcfce73	366	//
15e48021	367	//
	368	nclIn=0;
	369	trackIn = trackOut;
	370	trackIn.ResetCovariance(10.);
	371	//
	372	// Refit in one more time
	373	//
9dcfce73	374	for (Int_t irow=159; irow>0; irow--){
	375	AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
	376	if (!cl) continue;
	377	if (cl->GetX()<80) continue;
	378	Int_t sector = cl->GetDetector();
	379	Float_t dalpha = TMath::DegToRad()(sector%1820.+10.)-trackIn.GetAlpha();
91f8fa1d	380	if (TMath::Abs(dalpha)>0.01){
	381	if (!trackIn.Rotate(TMath::DegToRad()(sector%1820.+10.))) break;
	382	}
9dcfce73	383	Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
9dcfce73	384	Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
91350f5b	385	AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
	386	cov[0]*=cov[0];
	387	cov[2]*=cov[2];
	388	trackIn.GetXYZ(xyz);
9dcfce73	389	Double_t bz = AliTracker::GetBz(xyz);
9dcfce73	390
15e48021	391	if (!trackIn.PropagateTo(r[0],bz)) continue;
42b40d07	392	if (RejectCluster(cl,&trackIn)) continue;
15e48021	393	nclIn++;
9dcfce73	394	trackIn.Update(&r[1],cov);
9dcfce73	395	}
15e48021	396
15e48021	397
3bf5c7a6	398	trackIn.Rotate(trackInOld->GetAlpha());
	399	trackOut.Rotate(trackOutOld->GetAlpha());
	400	//
	401	trackInOld->GetXYZ(xyz);
	402	Double_t bz = AliTracker::GetBz(xyz);
	403	trackIn.PropagateTo(trackInOld->GetX(),bz);
	404	//
	405	trackOutOld->GetXYZ(xyz);
	406	bz = AliTracker::GetBz(xyz);
	407	trackOut.PropagateTo(trackOutOld->GetX(),bz);
	408
9dcfce73	409
	410	if (fStreamLevel>0){
	411	TTreeSRedirector *cstream = GetDebugStreamer();
	412	if (cstream){
	413	(*cstream)<<"Tracks"<<
108953e9	414	"run="<<fRun<< // run number
	415	"event="<<fEvent<< // event number
	416	"time="<<fTime<< // time stamp of event
	417	"trigger="<<fTrigger<< // trigger
15e48021	418	"triggerClass="<<&fTriggerClass<< // trigger
108953e9	419	"mag="<<fMagF<< // magnetic field
9dcfce73	420	"nclIn="<<nclIn<<
	421	"nclOut="<<nclOut<<
	422	"ncl="<<ncl<<
	423	"TrIn0.="<<trackInOld<<
	424	"TrOut0.="<<trackOutOld<<
	425	"TrIn1.="<<&trackIn<<
	426	"TrOut1.="<<&trackOut<<
	427	"\n";
	428	}
	429	}
3bf5c7a6	430	//
981e9de5	431	// And now rewrite ESDtrack and TPC seed
3bf5c7a6	432	//
	433
	434	(*trackInOld) = trackIn;
	435	(*trackOutOld) = trackOut;
	436	AliExternalTrackParam *t = &trackIn;
	437	track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
981e9de5	438	seed->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
fe8454d9	439	seed->SetNumberOfClusters((nclIn+nclOut)/2);
2cfb8d90	440	return kTRUE;
9dcfce73	441	}
	442
	443
	444
	445	Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
	446	//
	447	// check the acceptance of cluster
	448	// Cut on edge effects
	449	//
15e48021	450	Float_t kEdgeCut=2.5;
	451	Float_t kSigmaCut=6;
	452
9dcfce73	453	Bool_t isReject = kFALSE;
	454	Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
	455	Float_t dist = edgeY - TMath::Abs(cl->GetY());
	456	if (param) dist = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
15e48021	457	if (dist<kEdgeCut) isReject=kTRUE;
	458
	459	Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
	460	AliTPCseed::GetError(cl, param,cov[0],cov[2]);
91f8fa1d	461	if (param->GetSigmaY2()<0 \|\| param->GetSigmaZ2()<0){
d3ce44cb	462	AliError("Wrong parameters");
	463	return kFALSE;
	464	}
15e48021	465	Double_t py = (cl->GetY()-param->GetY())/TMath::Sqrt(cov[0]*cov[0]+param->GetSigmaY2());
	466	Double_t pz = (cl->GetZ()-param->GetZ())/TMath::Sqrt(cov[2]*cov[2]+param->GetSigmaZ2());
	467	//
	468	if ((pypy+pzpz)>kSigmaCut*kSigmaCut) isReject=kTRUE;
	469
9dcfce73	470	return isReject;
	471	}
	472
	473
91350f5b	474