[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 * chain = tool.MakeChain("esd.txt","esdTree",0,5000);
  chain->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 "AliTPCcalibDB.h"
#include "AliTPCTransform.h"
#include "AliTPCclusterMI.h"
#include "AliTPCseed.h"

ClassImp(AliTPCcalibCalib)

AliTPCcalibCalib::AliTPCcalibCalib():
    AliTPCcalibBase()
{
  //
  // Constructor
  //
}


AliTPCcalibCalib::AliTPCcalibCalib(const Text_t *name, const Text_t *title) 
  :AliTPCcalibBase()
{  
  SetName(name);
  SetTitle(title);
}


AliTPCcalibCalib::AliTPCcalibCalib(const AliTPCcalibCalib&calib):
  AliTPCcalibBase(calib)
{
  //
  // 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();   
  //
  //
  //

  for (Int_t i=0;i<ntracks;++i) {
    AliESDtrack *track = event->GetTrack(i);  
    const AliExternalTrackParam * trackIn = track->GetInnerParam();
    const AliExternalTrackParam * trackOut = track->GetOuterParam();
    if (!trackIn) continue;
    if (!trackOut) continue;
   
    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
    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);
  }
  return;
}

Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
  //
  // Refit track
  //

  //
  // First apply calibration
  //

  AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform();
  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);
    //
    cluster->SetX(x[0]);
    cluster->SetY(x[1]);
    cluster->SetZ(x[2]);
    if (fStreamLevel>2){
      TTreeSRedirector *cstream = GetDebugStreamer();
      if (cstream){
	(*cstream)<<"Clusters"<<
	  "cl0.="<<&cl0<<
	  "cl.="<<cluster<<
	  "\n";
      }
    }
  }
  Int_t ncl = seed->GetNumberOfClusters();  
  // 
  // And now do refit
  //
  AliExternalTrackParam * trackInOld  = (AliExternalTrackParam*)track->GetTPCInnerParam();
  AliExternalTrackParam * trackOutOld = (AliExternalTrackParam*)track->GetOuterParam();

  AliExternalTrackParam trackIn  = *trackOutOld;
  AliExternalTrackParam trackOut = *trackInOld;
  trackIn.ResetCovariance(10.);
  trackOut.ResetCovariance(10.);
  Double_t xyz[3];
  Int_t nclIn=0,nclOut=0;
  //
  // Refit out
  //
  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)
      trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.));
    if (RejectCluster(cl,&trackOut)) continue;
    Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
    Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
    trackOut.GetXYZ(xyz);
    Double_t bz = AliTracker::GetBz(xyz);
    if (trackOut.PropagateTo(r[0],bz)) nclOut++;
    trackOut.Update(&r[1],cov);    
    
  }
  //
  // 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)
      trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.));
    if (RejectCluster(cl,&trackIn)) continue;
    Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
    Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
    trackOut.GetXYZ(xyz);
    Double_t bz = AliTracker::GetBz(xyz);

    if (trackIn.PropagateTo(r[0],bz)) 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"<<
	"nclIn="<<nclIn<<
	"nclOut="<<nclOut<<
	"ncl="<<ncl<<
	"TrIn0.="<<trackInOld<<
	"TrOut0.="<<trackOutOld<<
	"TrIn1.="<<&trackIn<<
	"TrOut1.="<<&trackOut<<
	"\n";
    }
  }
  //
  // And now rewrite ESDtrack
  //
 
  (*trackInOld)  = trackIn;
  (*trackOutOld) = trackOut;
  AliExternalTrackParam *t = &trackIn;
  track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
  return kTRUE;
}


Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
  //
  // check the acceptance of cluster
  // Cut on edge effects
  //
  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<3) isReject=kTRUE;
  if (cl->GetType()<0) 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;
	34	TChain * chain = tool.MakeChain("esd.txt","esdTree",0,5000);
	35	chain->Lookup();
	36
	37
	38
	39	*/
	40
	41
	42
9dcfce73	43	// marian.ivanov@cern.ch
	44	//
	45	#include "AliTPCcalibCalib.h"
	46	#include "TSystem.h"
	47	#include "TFile.h"
	48	#include "TTreeStream.h"
	49	#include "AliLog.h"
	50	#include "TTimeStamp.h"
	51	#include "AliESDEvent.h"
	52	#include "AliESDfriend.h"
	53	#include "AliESDtrack.h"
	54	#include "AliTracker.h"
	55
	56	#include "AliTPCcalibDB.h"
	57	#include "AliTPCTransform.h"
	58	#include "AliTPCclusterMI.h"
	59	#include "AliTPCseed.h"
	60
	61	ClassImp(AliTPCcalibCalib)
	62
	63	AliTPCcalibCalib::AliTPCcalibCalib():
	64	AliTPCcalibBase()
	65	{
	66	//
	67	// Constructor
	68	//
	69	}
	70
	71
	72	AliTPCcalibCalib::AliTPCcalibCalib(const Text_t name, const Text_t title)
	73	:AliTPCcalibBase()
	74	{
	75	SetName(name);
	76	SetTitle(title);
	77	}
	78
	79
	80	AliTPCcalibCalib::AliTPCcalibCalib(const AliTPCcalibCalib&calib):
	81	AliTPCcalibBase(calib)
	82	{
	83	//
	84	// copy constructor
	85	//
	86	}
	87
	88	AliTPCcalibCalib &AliTPCcalibCalib::operator=(const AliTPCcalibCalib&calib){
	89	//
	90	//
	91	//
	92	((AliTPCcalibBase *)this)->operator=(calib);
	93	return *this;
	94	}
	95
	96
	97	AliTPCcalibCalib::~AliTPCcalibCalib() {
	98	//
	99	// destructor
	100	//
	101	}
	102
	103
	104	void AliTPCcalibCalib::Process(AliESDEvent *event){
	105	//
	106	//
107	//
108	if (!event) {
109	return;
110	}
111	AliESDfriend ESDfriend=static_cast<AliESDfriend>(event->FindListObject("AliESDfriend"));
112	if (!ESDfriend) {
113	return;
114	}
115
116	if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
117	Int_t ntracks=event->GetNumberOfTracks();
118	//
119	//
120	//
121
122	for (Int_t i=0;i<ntracks;++i) {
123	AliESDtrack *track = event->GetTrack(i);
124	const AliExternalTrackParam * trackIn = track->GetInnerParam();
125	const AliExternalTrackParam * trackOut = track->GetOuterParam();
126	if (!trackIn) continue;
127	if (!trackOut) continue;
128
129	AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
130	TObject *calibObject;
131	AliTPCseed *seed = 0;
132	for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
133	if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
134	}
135	if (!seed) continue;
136	RefitTrack(track, seed);
137	}
138	return;
139	}
140
141	Bool_t AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed){
142	//
143	// Refit track
144	//
145
146	//
147	// First apply calibration
148	//
149
150	AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform();
151	for (Int_t irow=0;irow<159;irow++) {
152	AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
153	if (!cluster) continue;
154	AliTPCclusterMI cl0(*cluster);
155	Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
156	Int_t i[1]={cluster->GetDetector()};
157	transform->Transform(x,i,0,1);
9dcfce73	158	//
	159	cluster->SetX(x[0]);
	160	cluster->SetY(x[1]);
	161	cluster->SetZ(x[2]);
	162	if (fStreamLevel>2){
	163	TTreeSRedirector *cstream = GetDebugStreamer();
	164	if (cstream){
	165	(*cstream)<<"Clusters"<<
	166	"cl0.="<<&cl0<<
	167	"cl.="<<cluster<<
	168	"\n";
	169	}
	170	}
	171	}
	172	Int_t ncl = seed->GetNumberOfClusters();
	173	//
	174	// And now do refit
	175	//
3bf5c7a6	176	AliExternalTrackParam * trackInOld = (AliExternalTrackParam*)track->GetTPCInnerParam();
9dcfce73	177	AliExternalTrackParam * trackOutOld = (AliExternalTrackParam*)track->GetOuterParam();
	178
	179	AliExternalTrackParam trackIn = *trackOutOld;
	180	AliExternalTrackParam trackOut = *trackInOld;
	181	trackIn.ResetCovariance(10.);
	182	trackOut.ResetCovariance(10.);
	183	Double_t xyz[3];
	184	Int_t nclIn=0,nclOut=0;
	185	//
	186	// Refit out
	187	//
	188	for (Int_t irow=0; irow<160; irow++){
	189	AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
	190	if (!cl) continue;
	191	if (cl->GetX()<80) continue;
	192	Int_t sector = cl->GetDetector();
	193	Float_t dalpha = TMath::DegToRad()(sector%1820.+10.)-trackOut.GetAlpha();
	194
	195	if (TMath::Abs(dalpha)>0.01)
	196	trackOut.Rotate(TMath::DegToRad()(sector%1820.+10.));
3bf5c7a6	197	if (RejectCluster(cl,&trackOut)) continue;
9dcfce73	198	Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
	199	Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
	200	trackOut.GetXYZ(xyz);
	201	Double_t bz = AliTracker::GetBz(xyz);
	202	if (trackOut.PropagateTo(r[0],bz)) nclOut++;
	203	trackOut.Update(&r[1],cov);
	204
	205	}
	206	//
	207	// Refit in
	208	//
	209
	210	for (Int_t irow=159; irow>0; irow--){
	211	AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
	212	if (!cl) continue;
	213	if (cl->GetX()<80) continue;
	214	Int_t sector = cl->GetDetector();
	215	Float_t dalpha = TMath::DegToRad()(sector%1820.+10.)-trackIn.GetAlpha();
	216	if (TMath::Abs(dalpha)>0.01)
	217	trackIn.Rotate(TMath::DegToRad()(sector%1820.+10.));
3bf5c7a6	218	if (RejectCluster(cl,&trackIn)) continue;
9dcfce73	219	Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
	220	Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
	221	trackOut.GetXYZ(xyz);
	222	Double_t bz = AliTracker::GetBz(xyz);
	223
	224	if (trackIn.PropagateTo(r[0],bz)) nclIn++;
	225	trackIn.Update(&r[1],cov);
	226	}
3bf5c7a6	227	trackIn.Rotate(trackInOld->GetAlpha());
	228	trackOut.Rotate(trackOutOld->GetAlpha());
	229	//
	230	trackInOld->GetXYZ(xyz);
	231	Double_t bz = AliTracker::GetBz(xyz);
	232	trackIn.PropagateTo(trackInOld->GetX(),bz);
	233	//
	234	trackOutOld->GetXYZ(xyz);
	235	bz = AliTracker::GetBz(xyz);
	236	trackOut.PropagateTo(trackOutOld->GetX(),bz);
	237
9dcfce73	238
	239	if (fStreamLevel>0){
	240	TTreeSRedirector *cstream = GetDebugStreamer();
	241	if (cstream){
	242	(*cstream)<<"Tracks"<<
	243	"nclIn="<<nclIn<<
	244	"nclOut="<<nclOut<<
	245	"ncl="<<ncl<<
	246	"TrIn0.="<<trackInOld<<
	247	"TrOut0.="<<trackOutOld<<
	248	"TrIn1.="<<&trackIn<<
	249	"TrOut1.="<<&trackOut<<
	250	"\n";
	251	}
	252	}
3bf5c7a6	253	//
	254	// And now rewrite ESDtrack
	255	//
	256
	257	(*trackInOld) = trackIn;
	258	(*trackOutOld) = trackOut;
	259	AliExternalTrackParam *t = &trackIn;
	260	track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
2cfb8d90	261	return kTRUE;
9dcfce73	262	}
	263
	264
	265
	266	Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
	267	//
	268	// check the acceptance of cluster
	269	// Cut on edge effects
	270	//
	271	Bool_t isReject = kFALSE;
	272	Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
	273	Float_t dist = edgeY - TMath::Abs(cl->GetY());
	274	if (param) dist = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
	275	if (dist<3) isReject=kTRUE;
	276	if (cl->GetType()<0) isReject=kTRUE;
	277	return isReject;
	278	}
	279
	280