Corrected UInt_t <-> Int_t conversion
[u/mrichter/AliRoot.git] / TPC / 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 "AliTPCPointCorrection.h"
66 #include <TGeoManager.h>
67 #include <TGeoPhysicalNode.h>
68 #include "TDatabasePDG.h"
69 ClassImp(AliTPCcalibCalib)
70
71 AliTPCcalibCalib::AliTPCcalibCalib():
72 AliTPCcalibBase(),
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
79
80 {
81   //
82   // Constructor
83   //
84 }
85
86
87 AliTPCcalibCalib::AliTPCcalibCalib(const Text_t *name, const Text_t *title) 
88   :AliTPCcalibBase(),
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
95 {  
96   SetName(name);
97   SetTitle(title);
98 }
99
100
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())
109
110 {
111   //
112   // copy constructor
113   //
114 }
115
116 AliTPCcalibCalib &AliTPCcalibCalib::operator=(const AliTPCcalibCalib&calib){
117   //
118   //
119   //
120   ((AliTPCcalibBase *)this)->operator=(calib);
121   return *this;
122 }
123
124
125 AliTPCcalibCalib::~AliTPCcalibCalib() {
126   //
127   // destructor
128   //
129 }
130
131
132 void     AliTPCcalibCalib::Process(AliESDEvent *event){
133   //
134   // 
135   //
136   if (!event) {
137     return;
138   }  
139   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
140   if (!ESDfriend) {
141    return;
142   }
143   
144   if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
145   Int_t ntracks=event->GetNumberOfTracks();   
146   //AliTPCcalibDB::Instance()->SetExBField(fMagF);
147
148   //
149   //
150   //
151
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;
168     }
169     if (!seed) continue;
170     RefitTrack(track, seed, event->GetMagneticField());
171     (*tpcOut)=*(track->GetOuterParam());  
172   }
173   return;
174 }
175
176 Bool_t  AliTPCcalibCalib::RefitTrack(AliESDtrack * track, AliTPCseed *seed, Float_t magesd){
177   //
178   // Refit track
179   // if magesd==0 forget the curvature
180
181   //
182   // 0 - Setup transform object
183   //
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));
190
191   static Int_t streamCounter=0;
192   streamCounter++;
193   AliESDfriendTrack *friendTrack = fCurrentFriendTrack;
194
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);
202   }
203   if(!fApplyTOFCorrection) { // disable TOF correction in transform
204     if(transform->GetCurrentRecoParam())
205       transform->GetCurrentRecoParamNonConst()->SetUseTOFCorrection(kFALSE);
206   }
207
208   //
209   // First apply calibration
210   //
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()};
218
219     transform->Transform(x,i,0,1);
220     //
221     // get position correction
222     //
223     Int_t ipad=0;
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());
227     //
228     cluster->SetX(x[0]);
229     cluster->SetY(x[1]);
230     cluster->SetZ(x[2]);
231     
232     //
233     // Apply alignemnt
234     //
235     if (transform->GetCurrentRecoParam()->GetUseSectorAlignment()){
236       if (!param->IsGeoRead()) param->ReadGeoMatrices();
237       TGeoHMatrix  *mat = param->GetClusterMatrix(cluster->GetDetector());
238       //TGeoHMatrix  mat;
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);
242       else{
243         // chack Loading of Geo matrices from GeoManager - TEMPORARY FIX
244       }
245       cluster->SetX(posC[0]);
246       cluster->SetY(posC[1]);
247       cluster->SetZ(posC[2]);
248     }
249
250
251
252     if (fStreamLevel>2 && streamCounter<20*fStreamLevel ){
253       // dump debug info if required
254       TTreeSRedirector *cstream = GetDebugStreamer();
255       if (cstream){
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
263           "cl0.="<<&cl0<<
264           "cl.="<<cluster<<
265           "cy="<<dy<<
266           "cz="<<dz<<
267           "\n";
268       }
269     }
270   }
271   //
272   //
273   //
274   Int_t ncl = seed->GetNumberOfClusters();
275   const Double_t kResetCov=4.;
276   const Double_t kSigma=5.;
277   Double_t covar[15];
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);
283   covar[14]=0.2*0.2;
284   if (TMath::Abs(magesd)<0.05) {
285      covar[14]=0.025*0.025;
286   }
287   // 
288   // And now do refit
289   //
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();
294
295     
296
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
303   }  
304
305   Double_t xyz[3];
306   Int_t nclIn=0,nclOut=0;
307   //
308   // Refit in
309   //
310   for (Int_t irow=159; irow>0; irow--){
311     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
312     if (!cl) continue;
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;
318     }
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]);
322     cov[0]*=cov[0];
323     cov[2]*=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));
326     trackIn.GetXYZ(xyz);
327     //    Double_t bz = AliTracker::GetBz(xyz);
328
329     //    if (!trackIn.PropagateTo(r[0],bz)) continue;
330     if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1.,kFALSE)) continue;
331
332     if (RejectCluster(cl,&trackIn)) continue;
333     nclIn++;
334     trackIn.Update(&r[1],cov);    
335   }
336   //
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
343   }
344
345   //
346   // Refit out
347   //
348   //Bool_t lastEdge=kFALSE;
349   for (Int_t irow=0; irow<160; irow++){
350     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
351     if (!cl) continue;
352     if (cl->GetX()<80) continue;
353     Int_t sector = cl->GetDetector();
354     Float_t dalpha = TMath::DegToRad()*(sector%18*20.+10.)-trackOut.GetAlpha();
355
356     if (TMath::Abs(dalpha)>0.01){
357       if (!trackOut.Rotate(TMath::DegToRad()*(sector%18*20.+10.))) break;
358     }
359     Double_t r[3]={cl->GetX(),cl->GetY(),cl->GetZ()};
360
361     Double_t cov[3]={0.01,0.,0.01}; //TODO: correct error parametrisation    
362     AliTPCseed::GetError(cl, &trackOut,cov[0],cov[2]);
363     cov[0]*=cov[0];
364     cov[2]*=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;
371
372     if (RejectCluster(cl,&trackOut)) continue;
373     nclOut++;
374     trackOut.Update(&r[1],cov);         
375     //if (cl->GetType()<0) lastEdge=kTRUE;
376     //if (cl->GetType()>=0) lastEdge=kFALSE;    
377   }
378   //
379   //
380   //
381   nclIn=0;
382   trackIn  = trackOut;
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
387   }
388   //
389   // Refit in one more time
390   //
391   for (Int_t irow=159; irow>0; irow--){
392     AliTPCclusterMI *cl=seed->GetClusterPointer(irow);
393     if (!cl) continue;
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;
399     }
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]);
403     cov[0]*=cov[0];
404     cov[2]*=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));
407
408     trackIn.GetXYZ(xyz);
409     //Double_t bz = AliTracker::GetBz(xyz);
410
411     //    if (!trackIn.PropagateTo(r[0],bz)) continue;
412     if (!AliTracker::PropagateTrackToBxByBz(&trackIn, r[0],mass,1,kFALSE)) continue;
413
414     if (RejectCluster(cl,&trackIn)) continue;
415     nclIn++;
416     trackIn.Update(&r[1],cov);    
417   }
418
419
420   trackIn.Rotate(trackInOld->GetAlpha());
421   trackOut.Rotate(trackOutOld->GetAlpha());
422   //
423   trackInOld->GetXYZ(xyz);
424   Double_t bz = AliTracker::GetBz(xyz);
425   trackIn.PropagateTo(trackInOld->GetX(),bz);
426   //
427   trackOutOld->GetXYZ(xyz);
428   bz = AliTracker::GetBz(xyz);  
429   trackOut.PropagateTo(trackOutOld->GetX(),bz);
430   
431
432   if (fStreamLevel>0 && streamCounter<100*fStreamLevel){
433     TTreeSRedirector *cstream = GetDebugStreamer();
434     if (cstream){
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
442         "nclIn="<<nclIn<<
443         "nclOut="<<nclOut<<
444         "ncl="<<ncl<<
445         "TrIn0.="<<trackInOld<<
446         "TrOut0.="<<trackOutOld<<
447         "TrIn1.="<<&trackIn<<
448         "TrOut1.="<<&trackOut<<
449         "\n";
450     }
451   }
452   //
453   // And now rewrite ESDtrack and TPC seed
454   //
455  
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);
463   return kTRUE;
464 }
465
466
467
468 Bool_t AliTPCcalibCalib::RejectCluster(AliTPCclusterMI* cl, AliExternalTrackParam * param){
469   //
470   // check the acceptance of cluster
471   // Cut on edge effects
472   //
473   if (!param) return kTRUE;
474   Float_t kEdgeCut=2.5;
475   Float_t kSigmaCut=6;
476
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;
482
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");
487     return kFALSE;
488   }
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());
491   //
492   if ((py*py+pz*pz)>kSigmaCut*kSigmaCut) isReject=kTRUE;
493   
494   return isReject;
495 }
496
497
498