1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 ///////////////////////////////////////////////////////////////////////////////
19 // Component for redoing the reconstruction from the clusters and tracks
21 // The new calibration data used
23 // In reality it overwrites the content of the ESD
28 gSystem->Load("libANALYSIS");
29 gSystem->Load("libTPCcalib");
31 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
32 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
33 AliXRDPROOFtoolkit tool;
34 TChain * chainCl = tool.MakeChain("calib.txt","Clusters",0,1);
36 TChain * chainTr = tool.MakeChain("calib.txt","Tracks",0,1);
45 // marian.ivanov@cern.ch
47 #include "AliTPCcalibCalib.h"
50 #include "TTreeStream.h"
52 #include "TTimeStamp.h"
53 #include "AliESDEvent.h"
54 #include "AliESDfriend.h"
55 #include "AliESDtrack.h"
56 #include "AliTracker.h"
57 #include "AliTPCClusterParam.h"
58 #include "AliTPCParam.h"
60 #include "AliTPCcalibDB.h"
61 #include "AliTPCTransform.h"
62 #include "AliTPCRecoParam.h"
63 #include "AliTPCclusterMI.h"
64 #include "AliTPCseed.h"
65 #include "AliTPCPointCorrection.h"
66 #include <TGeoManager.h>
67 #include <TGeoPhysicalNode.h>
68 #include "TDatabasePDG.h"
69 ClassImp(AliTPCcalibCalib)
71 AliTPCcalibCalib::AliTPCcalibCalib():
73 fApplyExBCorrection(1), // apply ExB correction
74 fApplyTOFCorrection(1), // apply TOF correction
75 fApplyPositionCorrection(1), // apply position correction
76 fApplySectorAlignment(1), // apply sector alignment
77 fApplyRPhiCorrection(1), // apply R-Phi correction
78 fApplyRCorrection(1) // apply Radial correction
87 AliTPCcalibCalib::AliTPCcalibCalib(const Text_t *name, const Text_t *title)
89 fApplyExBCorrection(1), // apply ExB correction
90 fApplyTOFCorrection(1), // apply TOF correction
91 fApplyPositionCorrection(1), // apply position correction
92 fApplySectorAlignment(1), // apply sector alignment
93 fApplyRPhiCorrection(1), // apply R-Phi correction
94 fApplyRCorrection(1) // apply Radial correction
101 AliTPCcalibCalib::AliTPCcalibCalib(const AliTPCcalibCalib&calib):
102 AliTPCcalibBase(calib),
103 fApplyExBCorrection(calib.GetApplyExBCorrection()),
104 fApplyTOFCorrection(calib.GetApplyTOFCorrection()),
105 fApplyPositionCorrection(calib.GetApplyPositionCorrection()),
106 fApplySectorAlignment(calib.GetApplySectorAlignment()),
107 fApplyRPhiCorrection(calib.GetApplyRPhiCorrection()),
108 fApplyRCorrection(calib.GetApplyRCorrection())
116 AliTPCcalibCalib &AliTPCcalibCalib::operator=(const AliTPCcalibCalib&calib){
120 ((AliTPCcalibBase *)this)->operator=(calib);
125 AliTPCcalibCalib::~AliTPCcalibCalib() {
132 void AliTPCcalibCalib::Process(AliESDEvent *event){
139 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
143 if (ESDfriend->TestSkipBit()) return;
144 if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
145 Int_t ntracks=ESDfriend->GetNumberOfTracks();
146 //AliTPCcalibDB::Instance()->SetExBField(fMagF);
152 for (Int_t i=0;i<ntracks;++i) {
153 AliESDtrack *track = event->GetTrack(i);
154 AliESDfriendTrack *friendTrack = (AliESDfriendTrack*) ESDfriend->GetTrack(i);
155 if (!friendTrack) continue;
156 //track->SetFriendTrack(friendTrack);
157 fCurrentFriendTrack=friendTrack;
158 const AliExternalTrackParam * trackIn = track->GetInnerParam();
159 const AliExternalTrackParam * trackOut = track->GetOuterParam();
160 AliExternalTrackParam * tpcOut = (AliExternalTrackParam *)friendTrack->GetTPCOut();
161 if (!trackIn) continue;
162 if (!trackOut) continue;
163 if (!tpcOut) continue;
164 TObject *calibObject;
165 AliTPCseed *seed = 0;
166 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
167 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
170 RefitTrack(track, seed, event->GetMagneticField());
171 (*tpcOut)=*(track->GetOuterParam());
176 Bool_t AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Float_t magesd){
179 // if magesd==0 forget the curvature
182 // 0 - Setup transform object
184 const Double_t kxIFC = 83.; // position of IFC
185 const Double_t kxOFC = 250.; // position of OFC
186 const Double_t kaFC = 1.; // amplitude
187 const Double_t ktFC = 5.0; // slope of error
188 //cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
189 //cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
191 static Int_t streamCounter=0;
193 AliESDfriendTrack *friendTrack = fCurrentFriendTrack;
195 AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
196 AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters();
197 transform->SetCurrentRun(fRun);
198 transform->SetCurrentTimeStamp((UInt_t)fTime);
199 if(!fApplyExBCorrection) { // disable ExB correction in transform
200 if(transform->GetCurrentRecoParam())
201 transform->GetCurrentRecoParamNonConst()->SetUseExBCorrection(0);
203 if(!fApplyTOFCorrection) { // disable TOF correction in transform
204 if(transform->GetCurrentRecoParam())
205 transform->GetCurrentRecoParamNonConst()->SetUseTOFCorrection(kFALSE);
209 // First apply calibration
211 // AliTPCPointCorrection * corr = AliTPCPointCorrection::Instance();
212 for (Int_t irow=0;irow<159;irow++) {
213 AliTPCclusterMI *cluster=seed->GetClusterPointer(irow);
214 if (!cluster) continue;
215 AliTPCclusterMI cl0(*cluster);
216 Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
217 Int_t i[1]={cluster->GetDetector()};
219 transform->Transform(x,i,0,1);
221 // get position correction
224 if (cluster->GetDetector()>35) ipad=1;
225 Float_t dy =0;//AliTPCClusterParam::SPosCorrection(0,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
226 Float_t dz =0;//AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
235 if (transform->GetCurrentRecoParam()->GetUseSectorAlignment()){
236 if (!param->IsGeoRead()) param->ReadGeoMatrices();
237 TGeoHMatrix *mat = param->GetClusterMatrix(cluster->GetDetector());
239 Double_t pos[3]= {cluster->GetX(),cluster->GetY(),cluster->GetZ()};
240 Double_t posC[3]={cluster->GetX(),cluster->GetY(),cluster->GetZ()};
241 if (mat) mat->LocalToMaster(pos,posC);
243 // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
245 cluster->SetX(posC[0]);
246 cluster->SetY(posC[1]);
247 cluster->SetZ(posC[2]);
252 if (fStreamLevel>2 && streamCounter<20*fStreamLevel ){
253 // dump debug info if required
254 TTreeSRedirector *cstream = GetDebugStreamer();
256 (*cstream)<<"Clusters"<<
257 "run="<<fRun<< // run number
258 "event="<<fEvent<< // event number
259 "time="<<fTime<< // time stamp of event
260 "trigger="<<fTrigger<< // trigger
261 "triggerClass="<<&fTriggerClass<< // trigger
262 "mag="<<fMagF<< // magnetic field
274 Int_t ncl = seed->GetNumberOfClusters();
275 const Double_t kResetCov=4.;
276 const Double_t kSigma=5.;
278 for (Int_t i=0;i<15;i++) covar[i]=0;
279 covar[0]=kSigma*kSigma;
280 covar[2]=kSigma*kSigma;
281 covar[5]=kSigma*kSigma/Float_t(ncl*ncl);
282 covar[9]=kSigma*kSigma/Float_t(ncl*ncl);
284 if (TMath::Abs(magesd)<0.05) {
285 covar[14]=0.025*0.025;
290 AliExternalTrackParam * trackInOld = (AliExternalTrackParam*)track->GetInnerParam();
291 AliExternalTrackParam * trackOuter = (AliExternalTrackParam*)track->GetOuterParam();
292 AliExternalTrackParam * trackOutOld = (AliExternalTrackParam *)friendTrack->GetTPCOut();
293 Double_t mass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
297 AliExternalTrackParam trackIn = *trackOutOld;
298 trackIn.ResetCovariance(kResetCov);
299 trackIn.AddCovariance(covar);
300 if (TMath::Abs(magesd)<0.05) {
301 ((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
302 ((Double_t&)(trackIn.GetCovariance()[14]))=covar[14]; // fix the line
306 Int_t nclIn=0,nclOut=0;
310 for (Int_t irow=159; irow>0; irow--){
311 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
313 if (cl->GetX()<80) continue;
314 Int_t sector = cl->GetDetector();
315 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
316 if (TMath::Abs(dalpha)>0.01){
317 if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
319 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
320 Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
321 AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
324 cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
325 cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
327 // Double_t bz = AliTracker::GetBz(xyz);
329 // if (!trackIn.PropagateTo(r[0],bz)) continue;
330 if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
332 if (RejectCluster(cl,&trackIn)) continue;
334 trackIn.Update(&r[1],cov);
337 AliExternalTrackParam trackOut = trackIn;
338 trackOut.ResetCovariance(kResetCov);
339 trackOut.AddCovariance(covar);
340 if (TMath::Abs(magesd)<0.05) {
341 ((Double_t&)(trackOut.GetParameter()[4]))=0.000000001;
342 ((Double_t&)(trackOut.GetCovariance()[14]))=covar[14]; // fix the line
348 //Bool_t lastEdge=kFALSE;
349 for (Int_t irow=0; irow<160; irow++){
350 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
352 if (cl->GetX()<80) continue;
353 Int_t sector = cl->GetDetector();
354 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
356 if (TMath::Abs(dalpha)>0.01){
357 if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
359 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
361 Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
362 AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
365 cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
366 cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
367 trackOut.GetXYZ(xyz);
368 //Double_t bz = AliTracker::GetBz(xyz);
369 // if (!trackOut.PropagateTo(r[0],bz)) continue;
370 if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;
372 if (RejectCluster(cl,&trackOut)) continue;
374 trackOut.Update(&r[1],cov);
375 //if (cl->GetType()<0) lastEdge=kTRUE;
376 //if (cl->GetType()>=0) lastEdge=kFALSE;
383 trackIn.ResetCovariance(kResetCov);
384 if (TMath::Abs(magesd)<0.05) {
385 ((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
386 ((Double_t&)(trackIn.GetCovariance()[14]))=covar[14]; // fix the line
389 // Refit in one more time
391 for (Int_t irow=159; irow>0; irow--){
392 AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
394 if (cl->GetX()<80) continue;
395 Int_t sector = cl->GetDetector();
396 Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackIn.GetAlpha();
397 if (TMath::Abs(dalpha)>0.01){
398 if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
400 Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
401 Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
402 AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
405 cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
406 cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
409 //Double_t bz = AliTracker::GetBz(xyz);
411 // if (!trackIn.PropagateTo(r[0],bz)) continue;
412 if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1,kFALSE)) continue;
414 if (RejectCluster(cl,&trackIn)) continue;
416 trackIn.Update(&r[1],cov);
420 trackIn.Rotate(trackInOld->GetAlpha());
421 trackOut.Rotate(trackOutOld->GetAlpha());
423 trackInOld->GetXYZ(xyz);
424 Double_t bz = AliTracker::GetBz(xyz);
425 trackIn.PropagateTo(trackInOld->GetX(),bz);
427 trackOutOld->GetXYZ(xyz);
428 bz = AliTracker::GetBz(xyz);
429 trackOut.PropagateTo(trackOutOld->GetX(),bz);
432 if (fStreamLevel>0 && streamCounter<100*fStreamLevel){
433 TTreeSRedirector *cstream = GetDebugStreamer();
435 (*cstream)<<"Tracks"<<
436 "run="<<fRun<< // run number
437 "event="<<fEvent<< // event number
438 "time="<<fTime<< // time stamp of event
439 "trigger="<<fTrigger<< // trigger
440 "triggerClass="<<&fTriggerClass<< // trigger
441 "mag="<<fMagF<< // magnetic field
445 "TrIn0.="<<trackInOld<<
446 "TrOut0.="<<trackOutOld<<
447 "TrIn1.="<<&trackIn<<
448 "TrOut1.="<<&trackOut<<
453 // And now rewrite ESDtrack and TPC seed
456 (*trackInOld) = trackIn;
457 (*trackOutOld) = trackOut;
458 (*trackOuter) = trackOut;
459 AliExternalTrackParam *t = &trackIn;
460 //track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
461 seed->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
462 seed->SetNumberOfClusters((nclIn+nclOut)/2);
468 Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
470 // check the acceptance of cluster
471 // Cut on edge effects
473 if (!param) return kTRUE;
474 Float_t kEdgeCut=2.5;
477 Bool_t isReject = kFALSE;
478 Float_t edgeY = cl->GetX()*TMath::Tan(TMath::Pi()/18);
479 Float_t dist = edgeY - TMath::Abs(cl->GetY());
480 dist = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
481 if (dist<kEdgeCut) isReject=kTRUE;
483 Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
484 AliTPCseed::GetError(cl, param,cov[0],cov[2]);
485 if (param->GetSigmaY2()<0 || param->GetSigmaZ2()<0){
486 AliError("Wrong parameters");
489 Double_t py = (cl->GetY()-param->GetY())/TMath::Sqrt(cov[0]*cov[0]+param->GetSigmaY2());
490 Double_t pz = (cl->GetZ()-param->GetZ())/TMath::Sqrt(cov[2]*cov[2]+param->GetSigmaZ2());
492 if ((py*py+pz*pz)>kSigmaCut*kSigmaCut) isReject=kTRUE;