.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / TPC / TPCcalib / AliTPCcalibCalib.cxx
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
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 * chainCl = tool.MakeChain("calib.txt","Clusters",0,1);
35   chainCl->Lookup();
36   TChain * chainTr = tool.MakeChain("calib.txt","Tracks",0,1);
37   chainTr->Lookup();
38
39
40
41 */
42
43
44
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"
57 #include "AliTPCClusterParam.h"
58 #include "AliTPCParam.h"
59
60 #include "AliTPCcalibDB.h"
61 #include "AliTPCTransform.h"
62 #include "AliTPCRecoParam.h"
63 #include "AliTPCclusterMI.h"
64 #include "AliTPCseed.h"
65 #include <TGeoManager.h>
66 #include <TGeoPhysicalNode.h>
67 #include "TDatabasePDG.h"
68 ClassImp(AliTPCcalibCalib)
69
70 AliTPCcalibCalib::AliTPCcalibCalib():
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
79 {
80   //
81   // Constructor
82   //
83 }
84
85
86 AliTPCcalibCalib::AliTPCcalibCalib(const Text_t *name, const Text_t *title) 
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
94 {  
95   SetName(name);
96   SetTitle(title);
97 }
98
99
100 AliTPCcalibCalib::AliTPCcalibCalib(const AliTPCcalibCalib&calib):
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
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   }
142   if (ESDfriend->TestSkipBit()) return;
143   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
144   Int_t ntracks=ESDfriend->GetNumberOfTracks();   
145   //AliTPCcalibDB::Instance()->SetExBField(fMagF);
146
147   //
148   //
149   //
150
151   for (Int_t i=0;i<ntracks;++i) {
152     AliESDtrack *track = event->GetTrack(i);     
153     AliESDfriendTrack *friendTrack = (AliESDfriendTrack*) ESDfriend->GetTrack(i);
154     if (!friendTrack) continue;
155     //track->SetFriendTrack(friendTrack);
156     fCurrentFriendTrack=friendTrack;
157     const AliExternalTrackParam * trackIn  = track->GetInnerParam();
158     const AliExternalTrackParam * trackOut = track->GetOuterParam();
159     AliExternalTrackParam * tpcOut   = (AliExternalTrackParam *)friendTrack->GetTPCOut();
160     if (!trackIn) continue;
161     if (!trackOut) continue;
162     if (!tpcOut) continue;   
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;
169     RefitTrack(track, seed, event->GetMagneticField());
170     (*tpcOut)=*(track->GetOuterParam());  
171   }
172   return;
173 }
174
175 Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Float_t magesd){
176   //
177   // Refit track
178   // if magesd==0 forget the curvature
179
180   //
181   // 0 - Setup transform object
182   //
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
190   static Int_t streamCounter=0;
191   streamCounter++;
192   AliESDfriendTrack *friendTrack = fCurrentFriendTrack;
193
194   AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ;
195   AliTPCParam     *param     = AliTPCcalibDB::Instance()->GetParameters();
196   transform->SetCurrentRun(fRun);
197   transform->SetCurrentTimeStamp((UInt_t)fTime);
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
207   //
208   // First apply calibration
209   //
210   //  AliTPCPointCorrection * corr =  AliTPCPointCorrection::Instance();
211   TVectorD vec(5, seed->GetParameter());
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]={static_cast<Double_t>(cluster->GetRow()),cluster->GetPad(),cluster->GetTimeBin()};
217     Int_t i[1]={cluster->GetDetector()};
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+ty*ty));      
224     }
225     transform->Transform(x,i,0,1);
226     //
227     // get position correction
228     //
229     Int_t ipad=0;
230     if (cluster->GetDetector()>35) ipad=1;
231     Float_t dy =0;//AliTPCClusterParam::SPosCorrection(0,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
232     Float_t dz =0;//AliTPCClusterParam::SPosCorrection(1,ipad,cluster->GetPad(),cluster->GetTimeBin(),cluster->GetZ(),cluster->GetSigmaY2(),cluster->GetSigmaZ2(),cluster->GetMax());
233     //
234     cluster->SetX(x[0]);
235     cluster->SetY(x[1]);
236     cluster->SetZ(x[2]);
237     
238     //
239     // Apply alignemnt
240     //
241     if (transform->GetCurrentRecoParam()->GetUseSectorAlignment()){
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     }
255
256
257
258     if (fStreamLevel>2 && gRandom->Rndm()<0.1 ){
259       // dump debug info if required
260       TTreeSRedirector *cstream = GetDebugStreamer();
261       if (cstream){
262         (*cstream)<<"Clusters"<<
263           "run="<<fRun<<              //  run number
264           "event="<<fEvent<<          //  event number
265           "time="<<fTime<<            //  time stamp of event
266           "trigger="<<fTrigger<<      //  trigger
267           "triggerClass="<<&fTriggerClass<<      //  trigger      
268           "mag="<<fMagF<<             //  magnetic field
269           "cl0.="<<&cl0<<
270           "cl.="<<cluster<<
271           "cy="<<dy<<
272           "cz="<<dz<<
273           "ty="<<ty<<
274           "tz="<<tz<<
275           "vec.="<<&vec<<  //track parameters
276           "\n";
277       }
278     }
279   }
280   //
281   //
282   //
283   Int_t ncl = seed->GetNumberOfClusters();
284   const Double_t kResetCov=4.;
285   const Double_t kSigma=5.;
286   Double_t covar[15];
287   for (Int_t i=0;i<15;i++) covar[i]=0;
288   covar[0]=kSigma*kSigma;
289   covar[2]=kSigma*kSigma;
290   covar[5]=kSigma*kSigma/Float_t(ncl*ncl);
291   covar[9]=kSigma*kSigma/Float_t(ncl*ncl);
292   covar[14]=0.2*0.2;
293   if (TMath::Abs(magesd)<0.05) {
294      covar[14]=0.025*0.025;
295   }
296   // 
297   // And now do refit
298   //
299   AliExternalTrackParam * trackInOld  = (AliExternalTrackParam*)track->GetInnerParam();
300   AliExternalTrackParam * trackOuter = (AliExternalTrackParam*)track->GetOuterParam();
301   AliExternalTrackParam * trackOutOld   = (AliExternalTrackParam *)friendTrack->GetTPCOut();
302   Double_t mass =    TDatabasePDG::Instance()->GetParticle("pi+")->Mass();
303
304     
305
306   AliExternalTrackParam trackIn  = *trackOutOld;
307   trackIn.ResetCovariance(kResetCov);
308   trackIn.AddCovariance(covar);
309   if (TMath::Abs(magesd)<0.05) {
310     ((Double_t&)(trackIn.GetParameter()[4]))=0.000000001;
311     ((Double_t&)(trackIn.GetCovariance()[14]))=covar[14];  // fix the line
312   }  
313
314   Double_t xyz[3];
315   Int_t nclIn=0,nclOut=0;
316   //
317   // Refit in
318   //
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%18*20.+10.)-trackIn.GetAlpha();
325     if (TMath::Abs(dalpha)>0.01){
326       if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
327     }
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];
333     cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
334     cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
335     trackIn.GetXYZ(xyz);
336     //    Double_t bz = AliTracker::GetBz(xyz);
337
338     //    if (!trackIn.PropagateTo(r[0],bz)) continue;
339     if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
340
341     if (RejectCluster(cl,&trackIn)) continue;
342     nclIn++;
343     trackIn.Update(&r[1],cov);    
344   }
345   //
346   AliExternalTrackParam trackOut = trackIn;
347   trackOut.ResetCovariance(kResetCov);
348   trackOut.AddCovariance(covar);
349   if (TMath::Abs(magesd)<0.05) {
350     ((Double_t&)(trackOut.GetParameter()[4]))=0.000000001;
351     ((Double_t&)(trackOut.GetCovariance()[14]))=covar[14];  // fix the line
352   }
353
354   //
355   // Refit out
356   //
357   //Bool_t lastEdge=kFALSE;
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%18*20.+10.)-trackOut.GetAlpha();
364
365     if (TMath::Abs(dalpha)>0.01){
366       if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
367     }
368     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
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];
374     cov[0]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
375     cov[2]+= kaFC*(TMath::Exp(-TMath::Abs(cl->GetX()-kxIFC)/ktFC)+TMath::Exp(-TMath::Abs(cl->GetX()-kxOFC)/ktFC));
376     trackOut.GetXYZ(xyz);
377     //Double_t bz = AliTracker::GetBz(xyz);
378     //    if (!trackOut.PropagateTo(r[0],bz)) continue;
379     if (!AliTracker::PropagateTrackToBxByBz(&trackOut, r[0],mass,1.,kFALSE)) continue;
380
381     if (RejectCluster(cl,&trackOut)) continue;
382     nclOut++;
383     trackOut.Update(&r[1],cov);         
384     //if (cl->GetType()<0) lastEdge=kTRUE;
385     //if (cl->GetType()>=0) lastEdge=kFALSE;    
386   }
387   //
388   //
389   //
390   nclIn=0;
391   trackIn  = trackOut;
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   }
397   //
398   // Refit in one more time
399   //
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%18*20.+10.)-trackIn.GetAlpha();
406     if (TMath::Abs(dalpha)>0.01){
407       if (!trackIn.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
408     }
409     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
410     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation
411     AliTPCseed::GetError(cl, &trackIn,cov[0],cov[2]);
412     cov[0]*=cov[0];
413     cov[2]*=cov[2];
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
417     trackIn.GetXYZ(xyz);
418     //Double_t bz = AliTracker::GetBz(xyz);
419
420     //    if (!trackIn.PropagateTo(r[0],bz)) continue;
421     if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1,kFALSE)) continue;
422
423     if (RejectCluster(cl,&trackIn)) continue;
424     nclIn++;
425     trackIn.Update(&r[1],cov);    
426   }
427
428
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   
440
441   if (fStreamLevel>0 && streamCounter<100*fStreamLevel){
442     TTreeSRedirector *cstream = GetDebugStreamer();
443     if (cstream){
444       (*cstream)<<"Tracks"<<
445         "run="<<fRun<<              //  run number
446         "event="<<fEvent<<          //  event number
447         "time="<<fTime<<            //  time stamp of event
448         "trigger="<<fTrigger<<      //  trigger
449         "triggerClass="<<&fTriggerClass<<      //  trigger
450         "mag="<<fMagF<<             //  magnetic field
451         "nclIn="<<nclIn<<
452         "nclOut="<<nclOut<<
453         "ncl="<<ncl<<
454         "seed.="<<seed<<
455         "track.="<<track<<
456         "TrIn0.="<<trackInOld<<
457         "TrOut0.="<<trackOutOld<<
458         "TrIn1.="<<&trackIn<<
459         "TrOut1.="<<&trackOut<<
460         "\n";
461     }
462   }
463   //
464   // And now rewrite ESDtrack and TPC seed
465   //
466  
467   (*trackInOld)  = trackIn;
468   (*trackOutOld) = trackOut;
469   (*trackOuter) = trackOut;
470   AliExternalTrackParam *t = &trackIn;
471   //track->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
472   seed->Set(t->GetX(),t->GetAlpha(),t->GetParameter(),t->GetCovariance());
473   seed->SetNumberOfClusters((nclIn+nclOut)/2);
474   return kTRUE;
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   //
484   if (!param) return kTRUE;
485   Float_t kEdgeCut=2.5;
486   Float_t kSigmaCut=6;
487
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());
491   dist  = TMath::Abs(edgeY - TMath::Abs(param->GetY()));
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]);
496   if (param->GetSigmaY2()<0 || param->GetSigmaZ2()<0){
497     AliError("Wrong parameters");
498     return kFALSE;
499   }
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 ((py*py+pz*pz)>kSigmaCut*kSigmaCut) isReject=kTRUE;
504   
505   return isReject;
506 }
507
508
509