[u/mrichter/AliRoot.git] / TPC / Calib / 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 "AliTPCParam.h"

#include "AliTPCcalibDB.h"
#include "AliTPCTransform.h"
#include "AliTPCRecoParam.h"
#include "AliTPCclusterMI.h"
#include "AliTPCseed.h"
#include <TGeoManager.h>
#include <TGeoPhysicalNode.h>
#include "TDatabasePDG.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 (ESDfriend->TestSkipBit()) return;
  if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
  Int_t ntracks=ESDfriend->GetNumberOfTracks();   
  //AliTPCcalibDB::Instance()->SetExBField(fMagF);

  //
  //
  //

  for (Int_t i=0;i<ntracks;++i) {
    AliESDtrack *track = event->GetTrack(i);     
    AliESDfriendTrack *friendTrack = (AliESDfriendTrack*) ESDfriend->GetTrack(i);
    if (!friendTrack) continue;
    //track->SetFriendTrack(friendTrack);
    fCurrentFriendTrack=friendTrack;
    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
  //
  const Double_t kxIFC = 83.;   // position of IFC
  const Double_t kxOFC = 250.;  // position of OFC
  const Double_t kaFC = 1.;    // amplitude 
  const Double_t ktFC = 5.0;    // slope of error 
  //cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
  //cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));

  static Int_t streamCounter=0;
  streamCounter++;
  AliESDfriendTrack *friendTrack = fCurrentFriendTrack;

  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
  AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
  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();
  TVectorD vec(5, seed->GetParameter());
  for (Int_t irow=0;irow<159;irow++) {
    AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
    if (!cluster) continue; 
    AliTPCclusterMI cl0(*cluster);
    Double_t x[3]={static_cast<Double_t>(cluster->GetRow()),cluster->GetPad(),cluster->GetTimeBin()};
    Int_t i[1]={cluster->GetDetector()};
    AliTPCTrackerPoint * point = seed->GetTrackPoint(irow);
    Double_t ty=0,tz=0;
     
    if (point){
      ty = TMath::Abs(point->GetAngleY());
      tz = TMath::Abs(point->GetAngleZ()*TMath::Sqrt(1+ty*ty));      
    }
    transform->Transform(x,i,0,1);
    //
    // get position correction
    //
    Int_t ipad=0;
    if (cluster->GetDetector()>35) ipad=1;
    Float_t dy =0;//AliTPCClusterParam::SPosCorrection(0,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
    Float_t dz =0;//AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
    //
    cluster->SetX(x[0]);
    cluster->SetY(x[1]);
    cluster->SetZ(x[2]);
    
    //
    // Apply alignemnt
    //
    if (transform->GetCurrentRecoParam()->GetUseSectorAlignment()){
      if (!param->IsGeoRead()) param->ReadGeoMatrices();
      TGeoHMatrix  *mat = param->GetClusterMatrix(cluster->GetDetector());
      //TGeoHMatrix  mat;
      Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
      Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
      if (mat) mat->LocalToMaster(pos,posC);
      else{
	// chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
      }
      cluster->SetX(posC[0]);
      cluster->SetY(posC[1]);
      cluster->SetZ(posC[2]);
    }


    if (fStreamLevel>2 && gRandom->Rndm()<0.1 ){
      // dump debug info if required
      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<<
	  "ty="<<ty<<
	  "tz="<<tz<<
	  "vec.="<<&vec<<  //track parameters
	  "\n";
      }
    }
  }
  //
  //
  //
  Int_t ncl = seed->GetNumberOfClusters();
  const Double_t kResetCov=4.;
  const Double_t kSigma=5.;
  Double_t covar[15];
  for (Int_t i=0;i<15;i++) covar[i]=0;
  covar[0]=kSigma*kSigma;
  covar[2]=kSigma*kSigma;
  covar[5]=kSigma*kSigma/Float_t(ncl*ncl);
  covar[9]=kSigma*kSigma/Float_t(ncl*ncl);
  covar[14]=0.2*0.2;
  if (TMath::Abs(magesd)<0.05) {
     covar[14]=0.025*0.025;
  }
  // 
  // And now do refit
  //
  AliExternalTrackParam * trackInOld  = (AliExternalTrackParam*)track->GetInnerParam();
  AliExternalTrackParam * trackOuter = (AliExternalTrackParam*)track->GetOuterParam();
  AliExternalTrackParam * trackOutOld   = (AliExternalTrackParam *)friendTrack->GetTPCOut();
  Double_t mass =    TDatabasePDG::Instance()->GetParticle("pi+")->Mass();

    
  AliExternalTrackParam trackIn  = *trackOutOld;
  trackIn.ResetCovariance(kResetCov);
  trackIn.AddCovariance(covar);
  if (TMath::Abs(magesd)<0.05) {
    ((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
    ((Double_t&)(trackIn.GetCovariance()[14]))=covar[14];  // fix the line
  }  

  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];
    cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
    cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
    trackIn.GetXYZ(xyz);
    //    Double_t bz = AliTracker::GetBz(xyz);

    //    if (!trackIn.PropagateTo(r[0],bz)) continue;
    if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;

    if (RejectCluster(cl,&trackIn)) continue;
    nclIn++;
    trackIn.Update(&r[1],cov);    
  }
  //
  AliExternalTrackParam trackOut = trackIn;
  trackOut.ResetCovariance(kResetCov);
  trackOut.AddCovariance(covar);
  if (TMath::Abs(magesd)<0.05) {
    ((Double_t&)(trackOut.GetParameter()[4]))=0.000000001;
    ((Double_t&)(trackOut.GetCovariance()[14]))=covar[14];  // fix the line
  }

  //
  // 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];
    cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
    cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
    trackOut.GetXYZ(xyz);
    //Double_t bz = AliTracker::GetBz(xyz);
    //    if (!trackOut.PropagateTo(r[0],bz)) continue;
    if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) 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(kResetCov);
  if (TMath::Abs(magesd)<0.05) {
    ((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
    ((Double_t&)(trackIn.GetCovariance()[14]))=covar[14];  // fix the line
  }
  //
  // 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];
    cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
    cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));

    trackIn.GetXYZ(xyz);
    //Double_t bz = AliTracker::GetBz(xyz);

    //    if (!trackIn.PropagateTo(r[0],bz)) continue;
    if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1,kFALSE)) 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 && streamCounter<100*fStreamLevel){
    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<<
	"seed.="<<seed<<
	"track.="<<track<<
	"TrIn0.="<<trackInOld<<
	"TrOut0.="<<trackOutOld<<
	"TrIn1.="<<&trackIn<<
	"TrOut1.="<<&trackOut<<
	"\n";
    }
  }
  //
  // And now rewrite ESDtrack and TPC seed
  //
 
  (*trackInOld)  = trackIn;
  (*trackOutOld) = trackOut;
  (*trackOuter) = 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
  //
  if (!param) return kTRUE;
  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());
  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"
b5738316	58	#include "AliTPCParam.h"
9dcfce73	59
	60	#include "AliTPCcalibDB.h"
	61	#include "AliTPCTransform.h"
7af539c6	62	#include "AliTPCRecoParam.h"
9dcfce73	63	#include "AliTPCclusterMI.h"
9dcfce73	64	#include "AliTPCseed.h"
b5738316	65	#include <TGeoManager.h>
b5738316	66	#include <TGeoPhysicalNode.h>
be67055b	67	#include "TDatabasePDG.h"
9dcfce73	68	ClassImp(AliTPCcalibCalib)
	69
	70	AliTPCcalibCalib::AliTPCcalibCalib():
7af539c6	71	AliTPCcalibBase(),
	72	fApplyExBCorrection(1), // apply ExB correction
	73	fApplyTOFCorrection(1), // apply TOF correction
	74	fApplyPositionCorrection(1), // apply position correction
	75	fApplySectorAlignment(1), // apply sector alignment
	76	fApplyRPhiCorrection(1), // apply R-Phi correction
	77	fApplyRCorrection(1) // apply Radial correction
	78
9dcfce73	79	{
	80	//
	81	// Constructor
	82	//
	83	}
	84
	85
	86	AliTPCcalibCalib::AliTPCcalibCalib(const Text_t name, const Text_t title)
7af539c6	87	:AliTPCcalibBase(),
	88	fApplyExBCorrection(1), // apply ExB correction
	89	fApplyTOFCorrection(1), // apply TOF correction
	90	fApplyPositionCorrection(1), // apply position correction
	91	fApplySectorAlignment(1), // apply sector alignment
	92	fApplyRPhiCorrection(1), // apply R-Phi correction
	93	fApplyRCorrection(1) // apply Radial correction
9dcfce73	94	{
	95	SetName(name);
	96	SetTitle(title);
	97	}
	98
	99
	100	AliTPCcalibCalib::AliTPCcalibCalib(const AliTPCcalibCalib&calib):
7af539c6	101	AliTPCcalibBase(calib),
	102	fApplyExBCorrection(calib.GetApplyExBCorrection()),
	103	fApplyTOFCorrection(calib.GetApplyTOFCorrection()),
	104	fApplyPositionCorrection(calib.GetApplyPositionCorrection()),
	105	fApplySectorAlignment(calib.GetApplySectorAlignment()),
	106	fApplyRPhiCorrection(calib.GetApplyRPhiCorrection()),
	107	fApplyRCorrection(calib.GetApplyRCorrection())
	108
9dcfce73	109	{
	110	//
	111	// copy constructor
	112	//
	113	}
	114
	115	AliTPCcalibCalib &AliTPCcalibCalib::operator=(const AliTPCcalibCalib&calib){
	116	//
	117	//
	118	//
	119	((AliTPCcalibBase *)this)->operator=(calib);
	120	return *this;
	121	}
	122
	123
	124	AliTPCcalibCalib::~AliTPCcalibCalib() {
	125	//
	126	// destructor
	127	//
	128	}
	129
	130
	131	void AliTPCcalibCalib::Process(AliESDEvent *event){
	132	//
	133	//
	134	//
	135	if (!event) {
	136	return;
	137	}
	138	AliESDfriend ESDfriend=static_cast<AliESDfriend>(event->FindListObject("AliESDfriend"));
	139	if (!ESDfriend) {
	140	return;
	141	}
d0c09227	142	if (ESDfriend->TestSkipBit()) return;
9dcfce73	143	if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
19d6a12a	144	Int_t ntracks=ESDfriend->GetNumberOfTracks();
91350f5b	145	//AliTPCcalibDB::Instance()->SetExBField(fMagF);
5c3e0d17	146
9dcfce73	147	//
	148	//
	149	//
	150
	151	for (Int_t i=0;i<ntracks;++i) {
3e55050f	152	AliESDtrack *track = event->GetTrack(i);
e527a1b9	153	AliESDfriendTrack friendTrack = (AliESDfriendTrack) ESDfriend->GetTrack(i);
e527a1b9	154	if (!friendTrack) continue;
4486a91f	155	//track->SetFriendTrack(friendTrack);
4486a91f	156	fCurrentFriendTrack=friendTrack;
3e55050f	157	const AliExternalTrackParam * trackIn = track->GetInnerParam();
9dcfce73	158	const AliExternalTrackParam * trackOut = track->GetOuterParam();
3e55050f	159	AliExternalTrackParam * tpcOut = (AliExternalTrackParam *)friendTrack->GetTPCOut();
9dcfce73	160	if (!trackIn) continue;
9dcfce73	161	if (!trackOut) continue;
3e55050f	162	if (!tpcOut) continue;
9dcfce73	163	TObject *calibObject;
	164	AliTPCseed *seed = 0;
	165	for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
	166	if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
	167	}
	168	if (!seed) continue;
e527a1b9	169	RefitTrack(track, seed, event->GetMagneticField());
3e55050f	170	(tpcOut)=(track->GetOuterParam());
9dcfce73	171	}
	172	return;
	173	}
	174
38b1a1ca	175	Bool_t AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Float_t magesd){
9dcfce73	176	//
9dcfce73	177	// Refit track
38b1a1ca	178	// if magesd==0 forget the curvature
9dcfce73	179
15e48021	180	//
	181	// 0 - Setup transform object
	182	//
76c58ee2	183	const Double_t kxIFC = 83.; // position of IFC
	184	const Double_t kxOFC = 250.; // position of OFC
	185	const Double_t kaFC = 1.; // amplitude
	186	const Double_t ktFC = 5.0; // slope of error
	187	//cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
	188	//cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
	189
4486a91f	190	static Int_t streamCounter=0;
	191	streamCounter++;
	192	AliESDfriendTrack *friendTrack = fCurrentFriendTrack;
e527a1b9	193
15e48021	194	AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
b96c3aef	195	AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
15e48021	196	transform->SetCurrentRun(fRun);
15e48021	197	transform->SetCurrentTimeStamp((UInt_t)fTime);
7af539c6	198	if(!fApplyExBCorrection) { // disable ExB correction in transform
	199	if(transform->GetCurrentRecoParam())
	200	transform->GetCurrentRecoParamNonConst()->SetUseExBCorrection(0);
	201	}
	202	if(!fApplyTOFCorrection) { // disable TOF correction in transform
	203	if(transform->GetCurrentRecoParam())
	204	transform->GetCurrentRecoParamNonConst()->SetUseTOFCorrection(kFALSE);
	205	}
	206
9dcfce73	207	//
	208	// First apply calibration
	209	//
064244d8	210	// AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
8fae6121	211	TVectorD vec(5, seed->GetParameter());
9dcfce73	212	for (Int_t irow=0;irow<159;irow++) {
	213	AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
	214	if (!cluster) continue;
	215	AliTPCclusterMI cl0(*cluster);
2942f542	216	Double_t x[3]={static_cast<Double_t>(cluster->GetRow()),cluster->GetPad(),cluster->GetTimeBin()};
9dcfce73	217	Int_t i[1]={cluster->GetDetector()};
8fae6121 MK	218	AliTPCTrackerPoint * point = seed->GetTrackPoint(irow);
	219	Double_t ty=0,tz=0;
	220
	221	if (point){
	222	ty = TMath::Abs(point->GetAngleY());
	223	tz = TMath::Abs(point->GetAngleZ()TMath::Sqrt(1+tyty));
	224	}
9dcfce73	225	transform->Transform(x,i,0,1);
11c7a854	226	//
	227	// get position correction
	228	//
	229	Int_t ipad=0;
	230	if (cluster->GetDetector()>35) ipad=1;
f26def9e	231	Float_t dy =0;//AliTPCClusterParam::SPosCorrection(0,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
f26def9e	232	Float_t dz =0;//AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
9dcfce73	233	//
	234	cluster->SetX(x[0]);
	235	cluster->SetY(x[1]);
	236	cluster->SetZ(x[2]);
b96c3aef	237
	238	//
	239	// Apply alignemnt
	240	//
b322e06a	241	if (transform->GetCurrentRecoParam()->GetUseSectorAlignment()){
b5738316	242	if (!param->IsGeoRead()) param->ReadGeoMatrices();
	243	TGeoHMatrix *mat = param->GetClusterMatrix(cluster->GetDetector());
	244	//TGeoHMatrix mat;
	245	Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
	246	Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
	247	if (mat) mat->LocalToMaster(pos,posC);
	248	else{
	249	// chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
	250	}
	251	cluster->SetX(posC[0]);
	252	cluster->SetY(posC[1]);
	253	cluster->SetZ(posC[2]);
	254	}
b96c3aef	255
	256
	257
8fae6121	258	if (fStreamLevel>2 && gRandom->Rndm()<0.1 ){
4486a91f	259	// dump debug info if required
9dcfce73	260	TTreeSRedirector *cstream = GetDebugStreamer();
	261	if (cstream){
	262	(*cstream)<<"Clusters"<<
108953e9	263	"run="<<fRun<< // run number
	264	"event="<<fEvent<< // event number
	265	"time="<<fTime<< // time stamp of event
	266	"trigger="<<fTrigger<< // trigger
15e48021	267	"triggerClass="<<&fTriggerClass<< // trigger
108953e9	268	"mag="<<fMagF<< // magnetic field
9dcfce73	269	"cl0.="<<&cl0<<
9dcfce73	270	"cl.="<<cluster<<
11c7a854	271	"cy="<<dy<<
11c7a854	272	"cz="<<dz<<
8fae6121 MK	273	"ty="<<ty<<
	274	"tz="<<tz<<
	275	"vec.="<<&vec<< //track parameters
9dcfce73	276	"\n";
	277	}
	278	}
	279	}
15e48021	280	//
	281	//
	282	//
11c7a854	283	Int_t ncl = seed->GetNumberOfClusters();
b5738316	284	const Double_t kResetCov=4.;
b5738316	285	const Double_t kSigma=5.;
11c7a854	286	Double_t covar[15];
11c7a854	287	for (Int_t i=0;i<15;i++) covar[i]=0;
b5738316	288	covar[0]=kSigma*kSigma;
	289	covar[2]=kSigma*kSigma;
	290	covar[5]=kSigmakSigma/Float_t(nclncl);
	291	covar[9]=kSigmakSigma/Float_t(nclncl);
	292	covar[14]=0.2*0.2;
38b1a1ca	293	if (TMath::Abs(magesd)<0.05) {
	294	covar[14]=0.025*0.025;
	295	}
9dcfce73	296	//
	297	// And now do refit
	298	//
11c7a854	299	AliExternalTrackParam * trackInOld = (AliExternalTrackParam*)track->GetInnerParam();
e527a1b9	300	AliExternalTrackParam * trackOuter = (AliExternalTrackParam*)track->GetOuterParam();
e527a1b9	301	AliExternalTrackParam * trackOutOld = (AliExternalTrackParam *)friendTrack->GetTPCOut();
be67055b	302	Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
be67055b	303
e527a1b9	304
3e55050f	305
9dcfce73	306	AliExternalTrackParam trackIn = *trackOutOld;
b5738316	307	trackIn.ResetCovariance(kResetCov);
b5738316	308	trackIn.AddCovariance(covar);
38b1a1ca	309	if (TMath::Abs(magesd)<0.05) {
38b1a1ca	310	((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
b5738316	311	((Double_t&)(trackIn.GetCovariance()[14]))=covar[14]; // fix the line
38b1a1ca	312	}
b5738316	313
9dcfce73	314	Double_t xyz[3];
	315	Int_t nclIn=0,nclOut=0;
	316	//
15e48021	317	// Refit in
15e48021	318	//
15e48021	319	for (Int_t irow=159; irow>0; irow--){
	320	AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
	321	if (!cl) continue;
	322	if (cl->GetX()<80) continue;
	323	Int_t sector = cl->GetDetector();
	324	Float_t dalpha = TMath::DegToRad()(sector%1820.+10.)-trackIn.GetAlpha();
91f8fa1d	325	if (TMath::Abs(dalpha)>0.01){
	326	if (!trackIn.Rotate(TMath::DegToRad()(sector%1820.+10.))) break;
	327	}
15e48021	328	Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
	329	Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
	330	AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
	331	cov[0]*=cov[0];
	332	cov[2]*=cov[2];
76c58ee2	333	cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
76c58ee2	334	cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
15e48021	335	trackIn.GetXYZ(xyz);
be67055b	336	// Double_t bz = AliTracker::GetBz(xyz);
	337
	338	// if (!trackIn.PropagateTo(r[0],bz)) continue;
b322e06a	339	if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
15e48021	340
15e48021	341	if (RejectCluster(cl,&trackIn)) continue;
	342	nclIn++;
	343	trackIn.Update(&r[1],cov);
	344	}
	345	//
	346	AliExternalTrackParam trackOut = trackIn;
b5738316	347	trackOut.ResetCovariance(kResetCov);
b5738316	348	trackOut.AddCovariance(covar);
38b1a1ca	349	if (TMath::Abs(magesd)<0.05) {
38b1a1ca	350	((Double_t&)(trackOut.GetParameter()[4]))=0.000000001;
b5738316	351	((Double_t&)(trackOut.GetCovariance()[14]))=covar[14]; // fix the line
38b1a1ca	352	}
b5738316	353
15e48021	354	//
9dcfce73	355	// Refit out
9dcfce73	356	//
15e48021	357	//Bool_t lastEdge=kFALSE;
9dcfce73	358	for (Int_t irow=0; irow<160; irow++){
	359	AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
	360	if (!cl) continue;
	361	if (cl->GetX()<80) continue;
	362	Int_t sector = cl->GetDetector();
	363	Float_t dalpha = TMath::DegToRad()(sector%1820.+10.)-trackOut.GetAlpha();
	364
91f8fa1d	365	if (TMath::Abs(dalpha)>0.01){
	366	if (!trackOut.Rotate(TMath::DegToRad()(sector%1820.+10.))) break;
	367	}
9dcfce73	368	Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
91350f5b	369
	370	Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
	371	AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
	372	cov[0]*=cov[0];
	373	cov[2]*=cov[2];
76c58ee2	374	cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
76c58ee2	375	cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
9dcfce73	376	trackOut.GetXYZ(xyz);
be67055b	377	//Double_t bz = AliTracker::GetBz(xyz);
be67055b	378	// if (!trackOut.PropagateTo(r[0],bz)) continue;
b322e06a	379	if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;
be67055b	380
42b40d07	381	if (RejectCluster(cl,&trackOut)) continue;
15e48021	382	nclOut++;
	383	trackOut.Update(&r[1],cov);
	384	//if (cl->GetType()<0) lastEdge=kTRUE;
	385	//if (cl->GetType()>=0) lastEdge=kFALSE;
9dcfce73	386	}
9dcfce73	387	//
9dcfce73	388	//
15e48021	389	//
	390	nclIn=0;
	391	trackIn = trackOut;
b5738316	392	trackIn.ResetCovariance(kResetCov);
	393	if (TMath::Abs(magesd)<0.05) {
	394	((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
	395	((Double_t&)(trackIn.GetCovariance()[14]))=covar[14]; // fix the line
	396	}
15e48021	397	//
	398	// Refit in one more time
	399	//
9dcfce73	400	for (Int_t irow=159; irow>0; irow--){
	401	AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
	402	if (!cl) continue;
	403	if (cl->GetX()<80) continue;
	404	Int_t sector = cl->GetDetector();
	405	Float_t dalpha = TMath::DegToRad()(sector%1820.+10.)-trackIn.GetAlpha();
91f8fa1d	406	if (TMath::Abs(dalpha)>0.01){
	407	if (!trackIn.Rotate(TMath::DegToRad()(sector%1820.+10.))) break;
	408	}
9dcfce73	409	Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
9dcfce73	410	Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
91350f5b	411	AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
	412	cov[0]*=cov[0];
	413	cov[2]*=cov[2];
76c58ee2	414	cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
	415	cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
	416
91350f5b	417	trackIn.GetXYZ(xyz);
be67055b	418	//Double_t bz = AliTracker::GetBz(xyz);
	419
	420	// if (!trackIn.PropagateTo(r[0],bz)) continue;
b322e06a	421	if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1,kFALSE)) continue;
9dcfce73	422
42b40d07	423	if (RejectCluster(cl,&trackIn)) continue;
15e48021	424	nclIn++;
9dcfce73	425	trackIn.Update(&r[1],cov);
9dcfce73	426	}
15e48021	427
15e48021	428
3bf5c7a6	429	trackIn.Rotate(trackInOld->GetAlpha());
	430	trackOut.Rotate(trackOutOld->GetAlpha());
	431	//
	432	trackInOld->GetXYZ(xyz);
	433	Double_t bz = AliTracker::GetBz(xyz);
	434	trackIn.PropagateTo(trackInOld->GetX(),bz);
	435	//
	436	trackOutOld->GetXYZ(xyz);
	437	bz = AliTracker::GetBz(xyz);
	438	trackOut.PropagateTo(trackOutOld->GetX(),bz);
	439
9dcfce73	440
4486a91f	441	if (fStreamLevel>0 && streamCounter<100*fStreamLevel){
9dcfce73	442	TTreeSRedirector *cstream = GetDebugStreamer();
	443	if (cstream){
	444	(*cstream)<<"Tracks"<<
108953e9	445	"run="<<fRun<< // run number
	446	"event="<<fEvent<< // event number
	447	"time="<<fTime<< // time stamp of event
	448	"trigger="<<fTrigger<< // trigger
15e48021	449	"triggerClass="<<&fTriggerClass<< // trigger
108953e9	450	"mag="<<fMagF<< // magnetic field
9dcfce73	451	"nclIn="<<nclIn<<
	452	"nclOut="<<nclOut<<
	453	"ncl="<<ncl<<
8fae6121 MK	454	"seed.="<<seed<<
8fae6121 MK	455	"track.="<<track<<
9dcfce73	456	"TrIn0.="<<trackInOld<<
	457	"TrOut0.="<<trackOutOld<<
	458	"TrIn1.="<<&trackIn<<
	459	"TrOut1.="<<&trackOut<<
	460	"\n";
	461	}
	462	}
3bf5c7a6	463	//
981e9de5	464	// And now rewrite ESDtrack and TPC seed
3bf5c7a6	465	//
	466
	467	(*trackInOld) = trackIn;
	468	(*trackOutOld) = trackOut;
e527a1b9	469	(*trackOuter) = trackOut;
3bf5c7a6	470	AliExternalTrackParam *t = &trackIn;
b322e06a	471	//track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
981e9de5	472	seed->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
fe8454d9	473	seed->SetNumberOfClusters((nclIn+nclOut)/2);
2cfb8d90	474	return kTRUE;
9dcfce73	475	}
	476
	477
	478
	479	Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
	480	//
	481	// check the acceptance of cluster
	482	// Cut on edge effects
	483	//
768e46f1	484	if (!param) return kTRUE;
15e48021	485	Float_t kEdgeCut=2.5;
	486	Float_t kSigmaCut=6;
	487
9dcfce73	488	Bool_t isReject = kFALSE;
	489	Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
	490	Float_t dist = edgeY - TMath::Abs(cl->GetY());
768e46f1	491	dist = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
15e48021	492	if (dist<kEdgeCut) isReject=kTRUE;
	493
	494	Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
	495	AliTPCseed::GetError(cl, param,cov[0],cov[2]);
91f8fa1d	496	if (param->GetSigmaY2()<0 \|\| param->GetSigmaZ2()<0){
d3ce44cb	497	AliError("Wrong parameters");
	498	return kFALSE;
	499	}
15e48021	500	Double_t py = (cl->GetY()-param->GetY())/TMath::Sqrt(cov[0]*cov[0]+param->GetSigmaY2());
	501	Double_t pz = (cl->GetZ()-param->GetZ())/TMath::Sqrt(cov[2]*cov[2]+param->GetSigmaZ2());
	502	//
	503	if ((pypy+pzpz)>kSigmaCut*kSigmaCut) isReject=kTRUE;
	504
9dcfce73	505	return isReject;
	506	}
	507
	508
91350f5b	509