update TRD code for PbPb 2011
[u/mrichter/AliRoot.git] / TRD / Cal / AliTRDCalTrkAttach.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 // Container for calibration parameters for AliTRDseedV1::AttachClusters()
19 // .... Longer description of content ...
20 //
21 // For calibration procedure check AliTRDtrackleOflHelper::CalibrateAttach()
22 // .... some reference to the calibration procedure .........
23 //                                                                      //
24 // Authors:                                                             //
25 //   Alex Bercuci <a.bercuci@gsi.de>                                    //
26 //                                                                      //
27 //////////////////////////////////////////////////////////////////////////
28
29 #include <TFile.h>
30 #include <TROOT.h>
31 #include <TMath.h>
32 #include "TGraphErrors.h"
33 #include "TObjArray.h"
34
35 #include "AliLog.h"
36
37 #include "AliTRDcalibDB.h"
38 #include "AliTRDgeometry.h"
39 #include "AliTRDCalTrkAttach.h"
40
41 ClassImp(AliTRDCalTrkAttach)
42
43 //______________________________________________________________________________
44 AliTRDCalTrkAttach::AliTRDCalTrkAttach()
45   :TNamed("AliTRDCalTrkAttach", "Calibration of AliTRDseedV1::AttachClusters")
46   ,fRClikeLimit(0.65)
47   ,fScaleCov(2.)
48   ,fLike(NULL)
49 {
50 // Default constructor
51
52   fNsgmDy[0] = 5; fNsgmDy[1] = 7;
53   fLikeMinRelDecrease[0] = .2; fLikeMinRelDecrease[1] = .3;
54 }
55
56 //______________________________________________________________________________
57 AliTRDCalTrkAttach::~AliTRDCalTrkAttach()
58 {
59 // Destructor
60   if(fLike) delete fLike;
61 }
62
63
64 //______________________________________________________________________________
65 Double_t AliTRDCalTrkAttach::CookLikelihood(Bool_t chg, Int_t ly, Float_t pt, Float_t phiTrk, Int_t n, Double_t dyr, Double_t dphi, Double_t sr) const
66 {
67 // Calculate likelihood for a segment to belong to a tracklet
68 // Based on calibrated values
69
70   if(n<4){
71     AliDebug(4, Form("Failed basic cut[s] : n[%d] ...", n));
72     return 0.;
73   }
74   //check likelihood array  
75   if (!fLike || fLike->GetEntries() != AliTRDgeometry::kNlayer*kNcharge*kNcalib) {
76     AliError("No usable AttachClusters calib object.");
77     return 0.;
78   }
79   Int_t offset(ly*kNcalib*kNcharge); // offset for layer
80   TGraphErrors *g(NULL);
81   if(!(g = (TGraphErrors*)fLike->At(offset+Int_t(kResPos)*kNcharge+Int_t(chg)))){
82     AliError("Failed retrieving AttachClusters graph.");
83     return 0.;
84   }
85   // Interpolate p_t
86   Int_t npts(g->GetN()), ip(0), jp(-1);
87   Double_t x0, y0, x1, y1, dd(0.), invdx(0.), f[4]={0., 0., 0., 0.};
88   for(Int_t kp(0); kp<npts; kp++){
89     g->GetPoint(kp, x1, y1);
90     if(x1>=pt){jp=kp; break;}
91   }
92   Bool_t boundary(kFALSE);
93   if(jp<0){
94     jp = npts-1; 
95     g->GetPoint(jp, x1, y1);
96     ip = npts-1;
97     boundary = kTRUE;
98   }else if(jp==0){ 
99     ip = 0;
100     boundary = kTRUE;
101   }else{ 
102     ip = jp-1;
103     g->GetPoint(ip, x0, y0);
104     invdx = 1./(x0-x1);
105   }
106   // process pt dependences
107   // +++ process dy
108   Double_t dym = boundary?y1:((pt*(y0-y1) + (x0*y1-x1*y0))*invdx),
109            sym = 0.5*(g->GetErrorY(ip)+g->GetErrorY(jp));
110   dd      = (dyr - dym)/sym; dd*=dd;
111   f[0] = TMath::Exp(-0.5*dd);
112   // +++ process dphi
113   if(!(g = (TGraphErrors*)fLike->At(offset+Int_t(kResAng)*kNcharge+Int_t(chg)))){
114     AliError("Failed retrieving AttachClusters graph.");
115     return 0.;
116   }
117   g->GetPoint(ip, x0, y0);g->GetPoint(jp, x1, y1);
118   Double_t dpm = boundary?y1:((pt*(y0-y1) + (x0*y1-x1*y0))*invdx),
119            spm = 0.5*(g->GetErrorY(ip)+g->GetErrorY(jp));
120   dd      = (dphi - dpm)/spm; dd*=dd;
121   f[1] = TMath::Exp(-0.5*dd);
122   // +++ process no of clusters
123   if(!(g = (TGraphErrors*)fLike->At(offset+Int_t(kNclMean)*kNcharge+Int_t(chg)))){
124     AliError("Failed retrieving AttachClusters graph.");
125     return 0.;
126   }
127   g->GetPoint(ip, x0, y0);g->GetPoint(jp, x1, y1);
128   Double_t nm = boundary?y1:((pt*(y0-y1) + (x0*y1-x1*y0))*invdx);
129   f[2] = (nm-TMath::Abs(n-nm))/nm;
130  
131   // process phi dependences
132   // +++ process <s>/s
133   if(!(g = (TGraphErrors*)fLike->At(offset+Int_t(kSigma)*kNcharge+Int_t(chg)))){
134     AliError("Failed retrieving AttachClusters graph.");
135     return 0.;
136   }
137   // Interpolate phi [deg]
138   npts=g->GetN(); jp=-1;
139   for(Int_t kp(0); kp<npts; kp++){
140     g->GetPoint(kp, x1, y1);
141     if(x1>=phiTrk){jp=kp; break;}
142   }
143   if(jp<0){
144     jp = npts-1; 
145     g->GetPoint(jp, x1, y1);
146     ip = jp;
147     boundary = kTRUE;
148   }else if(jp==0){ 
149     ip = jp;
150     boundary = kTRUE;
151   }else{ 
152     ip = jp-1;
153     g->GetPoint(ip, x0, y0);
154     invdx = 1./(x0-x1);
155     boundary = kFALSE;
156   }
157   Double_t sm = boundary?y1:((phiTrk*(y0-y1) + (x0*y1-x1*y0))*invdx),
158            ssm = 0.5*(g->GetErrorY(ip)+g->GetErrorY(jp));
159   dd      = (sr - sm)/ssm; dd*=dd;
160   f[3] = 1.;//TMath::Exp(-0.5*dd);
161
162   // Calculate likelihood
163   Double_t length = f[0]*f[0]+f[1]*f[1]+f[2]*f[2]+f[3]*f[3];
164   length = TMath::Sqrt(length);
165   Double_t cosTht = f[0]+f[1]+f[2]+f[3];
166   cosTht /= (4.*length);
167   AliDebug(2, Form("Like[%5.3f] ThtLike[%6.2f](deg)\n"
168     "    F_dy (%+6.2f)=%4.2f\n" 
169     "    F_phi(%+6.2f)=%4.2f\n"
170     "    F_ncl(%+6d)=%4.2f\n"
171     "    F_<s>(%+6.2f)=%4.2f", 
172     length, TMath::ACos(cosTht)*TMath::RadToDeg(),
173     dyr, f[0], dphi, f[1], n, f[2], sr, f[3]));
174
175   return length;
176 }
177
178 //______________________________________________________________________________
179 Bool_t AliTRDCalTrkAttach::LoadReferences(const Char_t *file)
180 {
181 // Load calibration data from file
182
183   if(!file || !TFile::Open(file)){
184     AliError("Parametrization file missing or unreadable.");
185     return kFALSE;
186   }
187   TGraphErrors *g(NULL);
188   Char_t co[kNcalib] = {'y', 'p', 's', 'n'},
189          cs[2] = {'n', 'p'};
190
191   if(fLike) fLike->Clear();
192   else fLike = new TObjArray(AliTRDgeometry::kNlayer*kNcharge*kNcalib);
193   for(Int_t ily(0); ily<AliTRDgeometry::kNlayer; ily++){
194     for(Int_t icalib(0); icalib<kNcalib; icalib++){
195       for(Int_t isgn(0); isgn<kNcharge; isgn++){
196         if(!(g = (TGraphErrors*)gFile->Get(Form("%c%c%d", co[icalib], cs[isgn], ily)))) return kFALSE;
197         fLike->AddAt(g, kNcharge*(ily*kNcalib+icalib)+isgn);
198       }
199     }
200   }
201   return kTRUE;
202 }