| 12ca5da1 |
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 | // // |
| dcf3a564 |
19 | // TPC cluster error, shape and charge parameterization as function |
| 20 | // of drift length, and inclination angle // |
| 21 | // |
| 22 | // Following notation is used in following |
| 23 | // Int_t dim 0 - y direction |
| 24 | // 1 - z direction |
| 25 | // |
| 26 | // Int_t type 0 - short pads |
| 27 | // 1 - medium pads |
| 28 | // 2 - long pads |
| 29 | // Float_t z - drift length |
| 30 | // |
| d028aade |
31 | // Float_t angle - tangent of inclination angle at given dimension |
| dcf3a564 |
32 | // |
| 33 | // Implemented parameterization |
| 34 | // |
| 35 | // |
| 36 | // 1. Resolution as function of drift length and inclination angle |
| 37 | // 1.a) GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle) |
| 38 | // Simple error parameterization as derived from analytical formula |
| 39 | // only linear term in drift length and angle^2 |
| 40 | // The formula is valid only with precission +-5% |
| 41 | // Separate parameterization for differnt pad geometry |
| 42 | // 1.b) GetError0Par |
| 43 | // Parabolic term correction - better precision |
| 44 | // |
| 45 | // 1.c) GetError1 - JUST FOR Study |
| 46 | // Similar to GetError1 |
| 47 | // The angular and diffusion effect is scaling with pad length |
| 48 | // common parameterization for different pad length |
| 49 | // |
| 96305e49 |
50 | // 2. Error parameterization using charge |
| 51 | // 2.a) GetErrorQ |
| 52 | // GetError0+ |
| 53 | // adding 1/Q component to diffusion and angluar part |
| 54 | // 2.b) GetErrorQPar |
| 55 | // GetError0Par+ |
| 56 | // adding 1/Q component to diffusion and angluar part |
| 57 | // 2.c) GetErrorQParScaled - Just for study |
| 58 | // One parameterization for all pad shapes |
| 59 | // Smaller precission as previous one |
| dcf3a564 |
60 | // |
| 96305e49 |
61 | // |
| d028aade |
62 | // Example how to retrieve the paramterization: |
| 63 | /* |
| d028aade |
64 | AliCDBManager::Instance()->SetDefaultStorage("local://$ALICE_ROOT"); |
| 236a0d03 |
65 | AliCDBManager::Instance()->SetRun(0) |
| d028aade |
66 | AliTPCClusterParam * param = AliTPCcalibDB::Instance()->GetClusterParam(); |
| 67 | |
| 68 | // |
| 69 | // |
| 70 | AliTPCClusterParam::SetInstance(param); |
| 71 | TF1 f1("f1","AliTPCClusterParam::SGetError0Par(1,0,x,0)",0,250); |
| d028aade |
72 | */ |
| 236a0d03 |
73 | |
| 74 | // EXAMPLE hot to create parameterization |
| 75 | /* |
| 76 | // Note resol is the resolution tree created by AliTPCcalibTracks |
| 77 | // |
| 78 | AliTPCClusterParam *param = new AliTPCClusterParam; |
| 79 | param->FitData(Resol); |
| 80 | AliTPCClusterParam::SetInstance(param); |
| 81 | |
| 82 | */ |
| 83 | |
| 96305e49 |
84 | // |
| d028aade |
85 | // // |
| 12ca5da1 |
86 | /////////////////////////////////////////////////////////////////////////////// |
| 87 | #include "AliTPCClusterParam.h" |
| 88 | #include "TMath.h" |
| 89 | #include "TFile.h" |
| 90 | #include "TTree.h" |
| 91 | #include <TVectorF.h> |
| 92 | #include <TLinearFitter.h> |
| 93 | #include <TH1F.h> |
| 94 | #include <TProfile2D.h> |
| 0a65832b |
95 | #include <TVectorD.h> |
| 96 | #include <TObjArray.h> |
| db2fdcfb |
97 | #include "AliTPCcalibDB.h" |
| 12ca5da1 |
98 | |
| 99 | ClassImp(AliTPCClusterParam) |
| 100 | |
| 101 | |
| 102 | AliTPCClusterParam* AliTPCClusterParam::fgInstance = 0; |
| 103 | |
| 104 | |
| 105 | /* |
| 106 | Example usage fitting parameterization: |
| 107 | TFile fres("resol.root"); //tree with resolution and shape |
| 108 | TTree * treeRes =(TTree*)fres.Get("Resol"); |
| 109 | |
| 110 | AliTPCClusterParam param; |
| 111 | param.SetInstance(¶m); |
| 112 | param.FitResol(treeRes); |
| 113 | param.FitRMS(treeRes); |
| 114 | TFile fparam("TPCClusterParam.root","recreate"); |
| 115 | param.Write("Param"); |
| 116 | // |
| 117 | // |
| 118 | TFile fparam("TPCClusterParam.root"); |
| 119 | AliTPCClusterParam *param2 = (AliTPCClusterParam *) fparam.Get("Param"); |
| 120 | param2->SetInstance(param2); |
| 121 | param2->Test(treeRes); |
| 122 | |
| 123 | |
| 124 | treeRes->Draw("(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol","Dim==0&&QMean<0") |
| 125 | |
| 126 | */ |
| 127 | |
| 128 | |
| 129 | |
| 130 | |
| 131 | //_ singleton implementation __________________________________________________ |
| 132 | AliTPCClusterParam* AliTPCClusterParam::Instance() |
| 133 | { |
| 134 | // |
| 135 | // Singleton implementation |
| 136 | // Returns an instance of this class, it is created if neccessary |
| 137 | // |
| 138 | if (fgInstance == 0){ |
| 139 | fgInstance = new AliTPCClusterParam(); |
| 140 | } |
| 141 | return fgInstance; |
| 142 | } |
| 143 | |
| 144 | |
| f1c2a4a3 |
145 | AliTPCClusterParam::AliTPCClusterParam(): |
| 146 | TObject(), |
| 38caa778 |
147 | fRatio(0), |
| b17540e4 |
148 | fQNorm(0), |
| 149 | fQpadTnorm(0), // q pad normalization - Total charge |
| 150 | fQpadMnorm(0) // q pad normalization - Max charge |
| 151 | // |
| f1c2a4a3 |
152 | { |
| 153 | // |
| 154 | // Default constructor |
| 155 | // |
| b17540e4 |
156 | fPosQTnorm[0] = 0; fPosQTnorm[1] = 0; fPosQTnorm[2] = 0; |
| 157 | fPosQMnorm[0] = 0; fPosQMnorm[1] = 0; fPosQMnorm[2] = 0; |
| 2e5bcb67 |
158 | // |
| 159 | fPosYcor[0] = 0; fPosYcor[1] = 0; fPosYcor[2] = 0; |
| 160 | fPosZcor[0] = 0; fPosZcor[1] = 0; fPosZcor[2] = 0; |
| f1c2a4a3 |
161 | } |
| 38caa778 |
162 | |
| 163 | AliTPCClusterParam::AliTPCClusterParam(const AliTPCClusterParam& param): |
| 164 | TObject(param), |
| 165 | fRatio(0), |
| b17540e4 |
166 | fQNorm(0), |
| 167 | fQpadTnorm(new TVectorD(*(param.fQpadTnorm))), // q pad normalization - Total charge |
| 168 | fQpadMnorm(new TVectorD(*(param.fQpadMnorm))) // q pad normalization - Max charge |
| 169 | |
| 38caa778 |
170 | { |
| 171 | // |
| 172 | // copy constructor |
| 173 | // |
| 174 | memcpy(this, ¶m,sizeof(AliTPCClusterParam)); |
| 175 | if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone(); |
| b17540e4 |
176 | // |
| 177 | if (param.fPosQTnorm[0]){ |
| 178 | fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0])); |
| 179 | fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1])); |
| 180 | fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2])); |
| 181 | // |
| 182 | fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0])); |
| 183 | fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1])); |
| 184 | fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2])); |
| 185 | } |
| 2e5bcb67 |
186 | if (param.fPosYcor[0]){ |
| 187 | fPosYcor[0] = new TVectorD(*(param.fPosYcor[0])); |
| 188 | fPosYcor[1] = new TVectorD(*(param.fPosYcor[1])); |
| 189 | fPosYcor[2] = new TVectorD(*(param.fPosYcor[2])); |
| 190 | // |
| 191 | fPosZcor[0] = new TVectorD(*(param.fPosZcor[0])); |
| 192 | fPosZcor[1] = new TVectorD(*(param.fPosZcor[1])); |
| 193 | fPosZcor[2] = new TVectorD(*(param.fPosZcor[2])); |
| 194 | } |
| 195 | |
| 38caa778 |
196 | } |
| 197 | |
| b17540e4 |
198 | |
| 38caa778 |
199 | AliTPCClusterParam & AliTPCClusterParam::operator=(const AliTPCClusterParam& param){ |
| 200 | // |
| 201 | // Assignment operator |
| 202 | // |
| 203 | if (this != ¶m) { |
| 204 | memcpy(this, ¶m,sizeof(AliTPCClusterParam)); |
| 205 | if (param.fQNorm) fQNorm = (TObjArray*) param.fQNorm->Clone(); |
| b17540e4 |
206 | if (param.fPosQTnorm[0]){ |
| 207 | fPosQTnorm[0] = new TVectorD(*(param.fPosQTnorm[0])); |
| 208 | fPosQTnorm[1] = new TVectorD(*(param.fPosQTnorm[1])); |
| 209 | fPosQTnorm[2] = new TVectorD(*(param.fPosQTnorm[2])); |
| 210 | // |
| 211 | fPosQMnorm[0] = new TVectorD(*(param.fPosQMnorm[0])); |
| 212 | fPosQMnorm[1] = new TVectorD(*(param.fPosQMnorm[1])); |
| 213 | fPosQMnorm[2] = new TVectorD(*(param.fPosQMnorm[2])); |
| 214 | } |
| 2e5bcb67 |
215 | if (param.fPosYcor[0]){ |
| 216 | fPosYcor[0] = new TVectorD(*(param.fPosYcor[0])); |
| 217 | fPosYcor[1] = new TVectorD(*(param.fPosYcor[1])); |
| 218 | fPosYcor[2] = new TVectorD(*(param.fPosYcor[2])); |
| 219 | // |
| 220 | fPosZcor[0] = new TVectorD(*(param.fPosZcor[0])); |
| 221 | fPosZcor[1] = new TVectorD(*(param.fPosZcor[1])); |
| 222 | fPosZcor[2] = new TVectorD(*(param.fPosZcor[2])); |
| 223 | } |
| 38caa778 |
224 | } |
| 225 | return *this; |
| 226 | } |
| 227 | |
| 228 | |
| f1c2a4a3 |
229 | AliTPCClusterParam::~AliTPCClusterParam(){ |
| 230 | // |
| 231 | // destructor |
| 232 | // |
| 233 | if (fQNorm) fQNorm->Delete(); |
| 234 | delete fQNorm; |
| b17540e4 |
235 | if (fPosQTnorm[0]){ |
| 236 | delete fPosQTnorm[0]; |
| 237 | delete fPosQTnorm[1]; |
| 238 | delete fPosQTnorm[2]; |
| 239 | // |
| 240 | delete fPosQMnorm[0]; |
| 241 | delete fPosQMnorm[1]; |
| 242 | delete fPosQMnorm[2]; |
| 243 | } |
| 2e5bcb67 |
244 | if (fPosYcor[0]){ |
| 245 | delete fPosYcor[0]; |
| 246 | delete fPosYcor[1]; |
| 247 | delete fPosYcor[2]; |
| 248 | // |
| 249 | delete fPosZcor[0]; |
| 250 | delete fPosZcor[1]; |
| 251 | delete fPosZcor[2]; |
| 252 | } |
| f1c2a4a3 |
253 | } |
| 12ca5da1 |
254 | |
| 255 | |
| 256 | void AliTPCClusterParam::FitResol0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){ |
| 257 | // |
| 258 | // Fit z - angular dependence of resolution |
| 259 | // |
| 260 | // Int_t dim=0, type=0; |
| 261 | char varVal[100]; |
| 262 | sprintf(varVal,"Resol:AngleM:Zm"); |
| 263 | char varErr[100]; |
| 264 | sprintf(varErr,"Sigma:AngleS:Zs"); |
| 265 | char varCut[100]; |
| 266 | sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type); |
| 267 | // |
| 268 | Int_t entries = tree->Draw(varVal,varCut); |
| 269 | Float_t px[10000], py[10000], pz[10000]; |
| 270 | Float_t ex[10000], ey[10000], ez[10000]; |
| 271 | // |
| 272 | tree->Draw(varErr,varCut); |
| 273 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 274 | ex[ipoint]= tree->GetV3()[ipoint]; |
| 275 | ey[ipoint]= tree->GetV2()[ipoint]; |
| 276 | ez[ipoint]= tree->GetV1()[ipoint]; |
| 277 | } |
| 278 | tree->Draw(varVal,varCut); |
| 279 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 280 | px[ipoint]= tree->GetV3()[ipoint]; |
| 281 | py[ipoint]= tree->GetV2()[ipoint]; |
| 282 | pz[ipoint]= tree->GetV1()[ipoint]; |
| 283 | } |
| 284 | |
| 285 | // |
| 286 | TLinearFitter fitter(3,"hyp2"); |
| 287 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 288 | Float_t val = pz[ipoint]*pz[ipoint]; |
| 289 | Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]); |
| 290 | Double_t x[2]; |
| 291 | x[0] = px[ipoint]; |
| 292 | x[1] = py[ipoint]*py[ipoint]; |
| 293 | fitter.AddPoint(x,val,err); |
| 294 | } |
| 295 | fitter.Eval(); |
| 296 | TVectorD param(3); |
| 297 | fitter.GetParameters(param); |
| 298 | param0[0] = param[0]; |
| 299 | param0[1] = param[1]; |
| 300 | param0[2] = param[2]; |
| 301 | Float_t chi2 = fitter.GetChisquare()/entries; |
| 302 | param0[3] = chi2; |
| 303 | error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2)); |
| 304 | error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2)); |
| 305 | error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2)); |
| 306 | } |
| 307 | |
| 308 | |
| 309 | void AliTPCClusterParam::FitResol0Par(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){ |
| 310 | // |
| 311 | // Fit z - angular dependence of resolution |
| 312 | // |
| 313 | // Int_t dim=0, type=0; |
| 314 | char varVal[100]; |
| 315 | sprintf(varVal,"Resol:AngleM:Zm"); |
| 316 | char varErr[100]; |
| 317 | sprintf(varErr,"Sigma:AngleS:Zs"); |
| 318 | char varCut[100]; |
| 319 | sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type); |
| 320 | // |
| 321 | Int_t entries = tree->Draw(varVal,varCut); |
| 322 | Float_t px[10000], py[10000], pz[10000]; |
| 323 | Float_t ex[10000], ey[10000], ez[10000]; |
| 324 | // |
| 325 | tree->Draw(varErr,varCut); |
| 326 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 327 | ex[ipoint]= tree->GetV3()[ipoint]; |
| 328 | ey[ipoint]= tree->GetV2()[ipoint]; |
| 329 | ez[ipoint]= tree->GetV1()[ipoint]; |
| 330 | } |
| 331 | tree->Draw(varVal,varCut); |
| 332 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 333 | px[ipoint]= tree->GetV3()[ipoint]; |
| 334 | py[ipoint]= tree->GetV2()[ipoint]; |
| 335 | pz[ipoint]= tree->GetV1()[ipoint]; |
| 336 | } |
| 337 | |
| 338 | // |
| 339 | TLinearFitter fitter(6,"hyp5"); |
| 340 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 341 | Float_t val = pz[ipoint]*pz[ipoint]; |
| 342 | Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]); |
| 343 | Double_t x[6]; |
| 344 | x[0] = px[ipoint]; |
| 345 | x[1] = py[ipoint]*py[ipoint]; |
| 346 | x[2] = x[0]*x[0]; |
| 347 | x[3] = x[1]*x[1]; |
| 348 | x[4] = x[0]*x[1]; |
| 349 | fitter.AddPoint(x,val,err); |
| 350 | } |
| 351 | fitter.Eval(); |
| 352 | TVectorD param(6); |
| 353 | fitter.GetParameters(param); |
| 354 | param0[0] = param[0]; |
| 355 | param0[1] = param[1]; |
| 356 | param0[2] = param[2]; |
| 357 | param0[3] = param[3]; |
| 358 | param0[4] = param[4]; |
| 359 | param0[5] = param[5]; |
| 360 | Float_t chi2 = fitter.GetChisquare()/entries; |
| 361 | param0[6] = chi2; |
| 362 | error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2)); |
| 363 | error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2)); |
| 364 | error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2)); |
| 365 | error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2)); |
| 366 | error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2)); |
| 367 | error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2)); |
| 368 | } |
| 369 | |
| 370 | |
| 371 | |
| 372 | |
| 373 | |
| 374 | void AliTPCClusterParam::FitResol1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){ |
| 375 | // |
| 376 | // Fit z - angular dependence of resolution - pad length scaling |
| 377 | // |
| 378 | // Int_t dim=0, type=0; |
| 379 | char varVal[100]; |
| 380 | sprintf(varVal,"Resol:AngleM*sqrt(Length):Zm/Length"); |
| 381 | char varErr[100]; |
| 382 | sprintf(varErr,"Sigma:AngleS:Zs"); |
| 383 | char varCut[100]; |
| 384 | sprintf(varCut,"Dim==%d&&QMean<0",dim); |
| 385 | // |
| 386 | Int_t entries = tree->Draw(varVal,varCut); |
| 387 | Float_t px[10000], py[10000], pz[10000]; |
| 388 | Float_t ex[10000], ey[10000], ez[10000]; |
| 389 | // |
| 390 | tree->Draw(varErr,varCut); |
| 391 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 392 | ex[ipoint]= tree->GetV3()[ipoint]; |
| 393 | ey[ipoint]= tree->GetV2()[ipoint]; |
| 394 | ez[ipoint]= tree->GetV1()[ipoint]; |
| 395 | } |
| 396 | tree->Draw(varVal,varCut); |
| 397 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 398 | px[ipoint]= tree->GetV3()[ipoint]; |
| 399 | py[ipoint]= tree->GetV2()[ipoint]; |
| 400 | pz[ipoint]= tree->GetV1()[ipoint]; |
| 401 | } |
| 402 | |
| 403 | // |
| 404 | TLinearFitter fitter(3,"hyp2"); |
| 405 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 406 | Float_t val = pz[ipoint]*pz[ipoint]; |
| 407 | Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]); |
| 408 | Double_t x[2]; |
| 409 | x[0] = px[ipoint]; |
| 410 | x[1] = py[ipoint]*py[ipoint]; |
| 411 | fitter.AddPoint(x,val,err); |
| 412 | } |
| 413 | fitter.Eval(); |
| 414 | TVectorD param(3); |
| 415 | fitter.GetParameters(param); |
| 416 | param0[0] = param[0]; |
| 417 | param0[1] = param[1]; |
| 418 | param0[2] = param[2]; |
| 419 | Float_t chi2 = fitter.GetChisquare()/entries; |
| 420 | param0[3] = chi2; |
| 421 | error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2)); |
| 422 | error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2)); |
| 423 | error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2)); |
| 424 | } |
| 425 | |
| 426 | void AliTPCClusterParam::FitResolQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){ |
| 427 | // |
| 428 | // Fit z - angular dependence of resolution - Q scaling |
| 429 | // |
| 430 | // Int_t dim=0, type=0; |
| 431 | char varVal[100]; |
| 432 | sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean"); |
| 433 | char varVal0[100]; |
| 434 | sprintf(varVal0,"Resol:AngleM:Zm"); |
| 435 | // |
| 436 | char varErr[100]; |
| 437 | sprintf(varErr,"Sigma:AngleS:Zs"); |
| 438 | char varCut[100]; |
| 439 | sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type); |
| 440 | // |
| 441 | Int_t entries = tree->Draw(varVal,varCut); |
| 442 | Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000]; |
| 443 | Float_t ex[20000], ey[20000], ez[20000]; |
| 444 | // |
| 445 | tree->Draw(varErr,varCut); |
| 446 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 447 | ex[ipoint]= tree->GetV3()[ipoint]; |
| 448 | ey[ipoint]= tree->GetV2()[ipoint]; |
| 449 | ez[ipoint]= tree->GetV1()[ipoint]; |
| 450 | } |
| 451 | tree->Draw(varVal,varCut); |
| 452 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 453 | px[ipoint]= tree->GetV3()[ipoint]; |
| 454 | py[ipoint]= tree->GetV2()[ipoint]; |
| 455 | pz[ipoint]= tree->GetV1()[ipoint]; |
| 456 | } |
| 457 | tree->Draw(varVal0,varCut); |
| 458 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 459 | pu[ipoint]= tree->GetV3()[ipoint]; |
| 460 | pt[ipoint]= tree->GetV2()[ipoint]; |
| 461 | } |
| 462 | |
| 463 | // |
| 464 | TLinearFitter fitter(5,"hyp4"); |
| 465 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 466 | Float_t val = pz[ipoint]*pz[ipoint]; |
| 467 | Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]); |
| 468 | Double_t x[4]; |
| 469 | x[0] = pu[ipoint]; |
| 470 | x[1] = pt[ipoint]*pt[ipoint]; |
| 471 | x[2] = px[ipoint]; |
| 472 | x[3] = py[ipoint]*py[ipoint]; |
| 473 | fitter.AddPoint(x,val,err); |
| 474 | } |
| 475 | |
| 476 | fitter.Eval(); |
| 477 | TVectorD param(5); |
| 478 | fitter.GetParameters(param); |
| 479 | param0[0] = param[0]; |
| 480 | param0[1] = param[1]; |
| 481 | param0[2] = param[2]; |
| 482 | param0[3] = param[3]; |
| 483 | param0[4] = param[4]; |
| 484 | Float_t chi2 = fitter.GetChisquare()/entries; |
| 485 | param0[5] = chi2; |
| 486 | error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2)); |
| 487 | error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2)); |
| 488 | error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2)); |
| 489 | error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2)); |
| 490 | error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2)); |
| 491 | } |
| 492 | |
| 493 | void AliTPCClusterParam::FitResolQPar(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){ |
| 494 | // |
| 495 | // Fit z - angular dependence of resolution - Q scaling - parabolic correction |
| 496 | // |
| 497 | // Int_t dim=0, type=0; |
| 498 | char varVal[100]; |
| 499 | sprintf(varVal,"Resol:AngleM/sqrt(QMean):Zm/QMean"); |
| 500 | char varVal0[100]; |
| 501 | sprintf(varVal0,"Resol:AngleM:Zm"); |
| 502 | // |
| 503 | char varErr[100]; |
| 504 | sprintf(varErr,"Sigma:AngleS:Zs"); |
| 505 | char varCut[100]; |
| 506 | sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type); |
| 507 | // |
| 508 | Int_t entries = tree->Draw(varVal,varCut); |
| 509 | Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000]; |
| 510 | Float_t ex[20000], ey[20000], ez[20000]; |
| 511 | // |
| 512 | tree->Draw(varErr,varCut); |
| 513 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 514 | ex[ipoint]= tree->GetV3()[ipoint]; |
| 515 | ey[ipoint]= tree->GetV2()[ipoint]; |
| 516 | ez[ipoint]= tree->GetV1()[ipoint]; |
| 517 | } |
| 518 | tree->Draw(varVal,varCut); |
| 519 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 520 | px[ipoint]= tree->GetV3()[ipoint]; |
| 521 | py[ipoint]= tree->GetV2()[ipoint]; |
| 522 | pz[ipoint]= tree->GetV1()[ipoint]; |
| 523 | } |
| 524 | tree->Draw(varVal0,varCut); |
| 525 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 526 | pu[ipoint]= tree->GetV3()[ipoint]; |
| 527 | pt[ipoint]= tree->GetV2()[ipoint]; |
| 528 | } |
| 529 | |
| 530 | // |
| 531 | TLinearFitter fitter(8,"hyp7"); |
| 532 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 533 | Float_t val = pz[ipoint]*pz[ipoint]; |
| 534 | Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]); |
| 535 | Double_t x[7]; |
| 536 | x[0] = pu[ipoint]; |
| 537 | x[1] = pt[ipoint]*pt[ipoint]; |
| 538 | x[2] = x[0]*x[0]; |
| 539 | x[3] = x[1]*x[1]; |
| 540 | x[4] = x[0]*x[1]; |
| 541 | x[5] = px[ipoint]; |
| 542 | x[6] = py[ipoint]*py[ipoint]; |
| 543 | // |
| 544 | fitter.AddPoint(x,val,err); |
| 545 | } |
| 546 | |
| 547 | fitter.Eval(); |
| 548 | TVectorD param(8); |
| 549 | fitter.GetParameters(param); |
| 550 | param0[0] = param[0]; |
| 551 | param0[1] = param[1]; |
| 552 | param0[2] = param[2]; |
| 553 | param0[3] = param[3]; |
| 554 | param0[4] = param[4]; |
| 555 | param0[5] = param[5]; |
| 556 | param0[6] = param[6]; |
| 557 | param0[7] = param[7]; |
| 558 | |
| 559 | Float_t chi2 = fitter.GetChisquare()/entries; |
| 560 | param0[8] = chi2; |
| 561 | error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2)); |
| 562 | error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2)); |
| 563 | error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2)); |
| 564 | error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2)); |
| 565 | error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2)); |
| 566 | error[5] = (fitter.GetParError(5)*TMath::Sqrt(chi2)); |
| 567 | error[6] = (fitter.GetParError(6)*TMath::Sqrt(chi2)); |
| 568 | error[7] = (fitter.GetParError(7)*TMath::Sqrt(chi2)); |
| 569 | } |
| 570 | |
| 571 | |
| 572 | |
| 573 | void AliTPCClusterParam::FitRMS0(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){ |
| 574 | // |
| 575 | // Fit z - angular dependence of resolution |
| 576 | // |
| 577 | // Int_t dim=0, type=0; |
| 578 | char varVal[100]; |
| 579 | sprintf(varVal,"RMSm:AngleM:Zm"); |
| 580 | char varErr[100]; |
| 581 | sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs"); |
| 582 | char varCut[100]; |
| 583 | sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type); |
| 584 | // |
| 585 | Int_t entries = tree->Draw(varVal,varCut); |
| 586 | Float_t px[10000], py[10000], pz[10000]; |
| 587 | Float_t ex[10000], ey[10000], ez[10000]; |
| 588 | // |
| 589 | tree->Draw(varErr,varCut); |
| 590 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 591 | ex[ipoint]= tree->GetV3()[ipoint]; |
| 592 | ey[ipoint]= tree->GetV2()[ipoint]; |
| 593 | ez[ipoint]= tree->GetV1()[ipoint]; |
| 594 | } |
| 595 | tree->Draw(varVal,varCut); |
| 596 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 597 | px[ipoint]= tree->GetV3()[ipoint]; |
| 598 | py[ipoint]= tree->GetV2()[ipoint]; |
| 599 | pz[ipoint]= tree->GetV1()[ipoint]; |
| 600 | } |
| 601 | |
| 602 | // |
| 603 | TLinearFitter fitter(3,"hyp2"); |
| 604 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 605 | Float_t val = pz[ipoint]*pz[ipoint]; |
| 606 | Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]); |
| 607 | Double_t x[2]; |
| 608 | x[0] = px[ipoint]; |
| 609 | x[1] = py[ipoint]*py[ipoint]; |
| 610 | fitter.AddPoint(x,val,err); |
| 611 | } |
| 612 | fitter.Eval(); |
| 613 | TVectorD param(3); |
| 614 | fitter.GetParameters(param); |
| 615 | param0[0] = param[0]; |
| 616 | param0[1] = param[1]; |
| 617 | param0[2] = param[2]; |
| 618 | Float_t chi2 = fitter.GetChisquare()/entries; |
| 619 | param0[3] = chi2; |
| 620 | error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2)); |
| 621 | error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2)); |
| 622 | error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2)); |
| 623 | } |
| 624 | |
| 625 | void AliTPCClusterParam::FitRMS1(TTree * tree, Int_t dim, Float_t *param0, Float_t *error){ |
| 626 | // |
| 627 | // Fit z - angular dependence of resolution - pad length scaling |
| 628 | // |
| 629 | // Int_t dim=0, type=0; |
| 630 | char varVal[100]; |
| 631 | sprintf(varVal,"RMSm:AngleM*Length:Zm"); |
| 632 | char varErr[100]; |
| 633 | sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Pad"); |
| 634 | char varCut[100]; |
| 635 | sprintf(varCut,"Dim==%d&&QMean<0",dim); |
| 636 | // |
| 637 | Int_t entries = tree->Draw(varVal,varCut); |
| 638 | Float_t px[10000], py[10000], pz[10000]; |
| 639 | Float_t type[10000], ey[10000], ez[10000]; |
| 640 | // |
| 641 | tree->Draw(varErr,varCut); |
| 642 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 643 | type[ipoint] = tree->GetV3()[ipoint]; |
| 644 | ey[ipoint] = tree->GetV2()[ipoint]; |
| 645 | ez[ipoint] = tree->GetV1()[ipoint]; |
| 646 | } |
| 647 | tree->Draw(varVal,varCut); |
| 648 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 649 | px[ipoint]= tree->GetV3()[ipoint]; |
| 650 | py[ipoint]= tree->GetV2()[ipoint]; |
| 651 | pz[ipoint]= tree->GetV1()[ipoint]; |
| 652 | } |
| 653 | |
| 654 | // |
| 655 | TLinearFitter fitter(4,"hyp3"); |
| 656 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 657 | Float_t val = pz[ipoint]*pz[ipoint]; |
| 658 | Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]); |
| 659 | Double_t x[3]; |
| 660 | x[0] = (type[ipoint]<0.5)? 0.:1.; |
| 661 | x[1] = px[ipoint]; |
| 662 | x[2] = py[ipoint]*py[ipoint]; |
| 663 | fitter.AddPoint(x,val,err); |
| 664 | } |
| 665 | fitter.Eval(); |
| 666 | TVectorD param(4); |
| 667 | fitter.GetParameters(param); |
| 668 | param0[0] = param[0]; |
| 669 | param0[1] = param[0]+param[1]; |
| 670 | param0[2] = param[2]; |
| 671 | param0[3] = param[3]; |
| 672 | Float_t chi2 = fitter.GetChisquare()/entries; |
| 673 | param0[4] = chi2; |
| 674 | error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2)); |
| 675 | error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2)); |
| 676 | error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2)); |
| 677 | error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2)); |
| 678 | } |
| 679 | |
| 680 | void AliTPCClusterParam::FitRMSQ(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t *error){ |
| 681 | // |
| 682 | // Fit z - angular dependence of resolution - Q scaling |
| 683 | // |
| 684 | // Int_t dim=0, type=0; |
| 685 | char varVal[100]; |
| 686 | sprintf(varVal,"RMSm:AngleM/sqrt(QMean):Zm/QMean"); |
| 687 | char varVal0[100]; |
| 688 | sprintf(varVal0,"RMSm:AngleM:Zm"); |
| 689 | // |
| 690 | char varErr[100]; |
| 691 | sprintf(varErr,"sqrt((1./(100.*sqrt(12.))^2)+RMSe0^2):AngleS:Zs"); |
| 692 | char varCut[100]; |
| 693 | sprintf(varCut,"Dim==%d&&Pad==%d&&QMean>0",dim,type); |
| 694 | // |
| 695 | Int_t entries = tree->Draw(varVal,varCut); |
| 696 | Float_t px[20000], py[20000], pz[20000], pu[20000], pt[20000]; |
| 697 | Float_t ex[20000], ey[20000], ez[20000]; |
| 698 | // |
| 699 | tree->Draw(varErr,varCut); |
| 700 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 701 | ex[ipoint]= tree->GetV3()[ipoint]; |
| 702 | ey[ipoint]= tree->GetV2()[ipoint]; |
| 703 | ez[ipoint]= tree->GetV1()[ipoint]; |
| 704 | } |
| 705 | tree->Draw(varVal,varCut); |
| 706 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 707 | px[ipoint]= tree->GetV3()[ipoint]; |
| 708 | py[ipoint]= tree->GetV2()[ipoint]; |
| 709 | pz[ipoint]= tree->GetV1()[ipoint]; |
| 710 | } |
| 711 | tree->Draw(varVal0,varCut); |
| 712 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 713 | pu[ipoint]= tree->GetV3()[ipoint]; |
| 714 | pt[ipoint]= tree->GetV2()[ipoint]; |
| 715 | } |
| 716 | |
| 717 | // |
| 718 | TLinearFitter fitter(5,"hyp4"); |
| 719 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 720 | Float_t val = pz[ipoint]*pz[ipoint]; |
| 721 | Float_t err = 2*pz[ipoint]*TMath::Sqrt(ez[ipoint]*ez[ipoint]+fRatio*fRatio*pz[ipoint]*pz[ipoint]); |
| 722 | Double_t x[4]; |
| 723 | x[0] = pu[ipoint]; |
| 724 | x[1] = pt[ipoint]*pt[ipoint]; |
| 725 | x[2] = px[ipoint]; |
| 726 | x[3] = py[ipoint]*py[ipoint]; |
| 727 | fitter.AddPoint(x,val,err); |
| 728 | } |
| 729 | |
| 730 | fitter.Eval(); |
| 731 | TVectorD param(5); |
| 732 | fitter.GetParameters(param); |
| 733 | param0[0] = param[0]; |
| 734 | param0[1] = param[1]; |
| 735 | param0[2] = param[2]; |
| 736 | param0[3] = param[3]; |
| 737 | param0[4] = param[4]; |
| 738 | Float_t chi2 = fitter.GetChisquare()/entries; |
| 739 | param0[5] = chi2; |
| 740 | error[0] = (fitter.GetParError(0)*TMath::Sqrt(chi2)); |
| 741 | error[1] = (fitter.GetParError(1)*TMath::Sqrt(chi2)); |
| 742 | error[2] = (fitter.GetParError(2)*TMath::Sqrt(chi2)); |
| 743 | error[3] = (fitter.GetParError(3)*TMath::Sqrt(chi2)); |
| 744 | error[4] = (fitter.GetParError(4)*TMath::Sqrt(chi2)); |
| 745 | } |
| 746 | |
| 747 | |
| 748 | void AliTPCClusterParam::FitRMSSigma(TTree * tree, Int_t dim, Int_t type, Float_t *param0, Float_t */*error*/){ |
| 749 | // |
| 750 | // Fit z - angular dependence of resolution - Q scaling |
| 751 | // |
| 752 | // Int_t dim=0, type=0; |
| 753 | char varVal[100]; |
| 754 | sprintf(varVal,"RMSs:RMSm"); |
| 755 | // |
| 756 | char varCut[100]; |
| 757 | sprintf(varCut,"Dim==%d&&Pad==%d&&QMean<0",dim,type); |
| 758 | // |
| 759 | Int_t entries = tree->Draw(varVal,varCut); |
| 760 | Float_t px[20000], py[20000]; |
| 761 | // |
| 762 | tree->Draw(varVal,varCut); |
| 763 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 764 | px[ipoint]= tree->GetV2()[ipoint]; |
| 765 | py[ipoint]= tree->GetV1()[ipoint]; |
| 766 | } |
| 767 | TLinearFitter fitter(2,"pol1"); |
| 768 | for (Int_t ipoint=0; ipoint<entries; ipoint++){ |
| 769 | Float_t val = py[ipoint]; |
| 770 | Float_t err = fRatio*px[ipoint]; |
| 771 | Double_t x[4]; |
| 772 | x[0] = px[ipoint]; |
| 236a0d03 |
773 | if (err>0) fitter.AddPoint(x,val,err); |
| 12ca5da1 |
774 | } |
| 775 | fitter.Eval(); |
| 776 | param0[0]= fitter.GetParameter(0); |
| 777 | param0[1]= fitter.GetParameter(1); |
| 778 | } |
| 779 | |
| 780 | |
| 781 | |
| 782 | Float_t AliTPCClusterParam::GetError0(Int_t dim, Int_t type, Float_t z, Float_t angle){ |
| 783 | // |
| 784 | // |
| 785 | // |
| 786 | Float_t value=0; |
| 787 | value += fParamS0[dim][type][0]; |
| 788 | value += fParamS0[dim][type][1]*z; |
| 789 | value += fParamS0[dim][type][2]*angle*angle; |
| 790 | value = TMath::Sqrt(TMath::Abs(value)); |
| 791 | return value; |
| 792 | } |
| 793 | |
| 794 | |
| 795 | Float_t AliTPCClusterParam::GetError0Par(Int_t dim, Int_t type, Float_t z, Float_t angle){ |
| 796 | // |
| 797 | // |
| 798 | // |
| 799 | Float_t value=0; |
| 800 | value += fParamS0Par[dim][type][0]; |
| 801 | value += fParamS0Par[dim][type][1]*z; |
| 802 | value += fParamS0Par[dim][type][2]*angle*angle; |
| 803 | value += fParamS0Par[dim][type][3]*z*z; |
| 804 | value += fParamS0Par[dim][type][4]*angle*angle*angle*angle; |
| 805 | value += fParamS0Par[dim][type][5]*z*angle*angle; |
| 806 | value = TMath::Sqrt(TMath::Abs(value)); |
| 807 | return value; |
| 808 | } |
| 809 | |
| 810 | |
| 811 | |
| 812 | Float_t AliTPCClusterParam::GetError1(Int_t dim, Int_t type, Float_t z, Float_t angle){ |
| 813 | // |
| 814 | // |
| 815 | // |
| 816 | Float_t value=0; |
| 817 | Float_t length=0.75; |
| 818 | if (type==1) length=1; |
| 819 | if (type==2) length=1.5; |
| 820 | value += fParamS1[dim][0]; |
| 821 | value += fParamS1[dim][1]*z/length; |
| 822 | value += fParamS1[dim][2]*angle*angle*length; |
| 823 | value = TMath::Sqrt(TMath::Abs(value)); |
| 824 | return value; |
| 825 | } |
| 826 | |
| 827 | Float_t AliTPCClusterParam::GetErrorQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){ |
| 828 | // |
| 829 | // |
| 830 | // |
| 831 | Float_t value=0; |
| 832 | value += fParamSQ[dim][type][0]; |
| 833 | value += fParamSQ[dim][type][1]*z; |
| 834 | value += fParamSQ[dim][type][2]*angle*angle; |
| 835 | value += fParamSQ[dim][type][3]*z/Qmean; |
| 836 | value += fParamSQ[dim][type][4]*angle*angle/Qmean; |
| 837 | value = TMath::Sqrt(TMath::Abs(value)); |
| 838 | return value; |
| 839 | |
| 840 | |
| 841 | } |
| 842 | |
| 843 | Float_t AliTPCClusterParam::GetErrorQPar(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){ |
| 844 | // |
| 845 | // |
| 846 | // |
| 847 | Float_t value=0; |
| 848 | value += fParamSQPar[dim][type][0]; |
| 849 | value += fParamSQPar[dim][type][1]*z; |
| 850 | value += fParamSQPar[dim][type][2]*angle*angle; |
| 851 | value += fParamSQPar[dim][type][3]*z*z; |
| 852 | value += fParamSQPar[dim][type][4]*angle*angle*angle*angle; |
| 853 | value += fParamSQPar[dim][type][5]*z*angle*angle; |
| 854 | value += fParamSQPar[dim][type][6]*z/Qmean; |
| 855 | value += fParamSQPar[dim][type][7]*angle*angle/Qmean; |
| 856 | value = TMath::Sqrt(TMath::Abs(value)); |
| 857 | return value; |
| 858 | |
| 859 | |
| 860 | } |
| 861 | |
| 862 | Float_t AliTPCClusterParam::GetErrorQParScaled(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){ |
| 863 | // |
| 864 | // |
| 865 | // |
| 866 | Float_t value=0; |
| 867 | value += fParamSQPar[dim][type][0]; |
| 868 | value += fParamSQPar[dim][type][1]*z; |
| 869 | value += fParamSQPar[dim][type][2]*angle*angle; |
| 870 | value += fParamSQPar[dim][type][3]*z*z; |
| 871 | value += fParamSQPar[dim][type][4]*angle*angle*angle*angle; |
| 872 | value += fParamSQPar[dim][type][5]*z*angle*angle; |
| 873 | value += fParamSQPar[dim][type][6]*z/Qmean; |
| 874 | value += fParamSQPar[dim][type][7]*angle*angle/Qmean; |
| 875 | Float_t valueMean = GetError0Par(dim,type,z,angle); |
| 876 | value -= 0.35*0.35*valueMean*valueMean; |
| 877 | value = TMath::Sqrt(TMath::Abs(value)); |
| 878 | return value; |
| 879 | |
| 880 | |
| 881 | } |
| 882 | |
| 883 | Float_t AliTPCClusterParam::GetRMS0(Int_t dim, Int_t type, Float_t z, Float_t angle){ |
| 884 | // |
| 885 | // calculate mean RMS of cluster - z,angle - parameters for each pad and dimension separatelly |
| 886 | // |
| 887 | Float_t value=0; |
| 888 | value += fParamRMS0[dim][type][0]; |
| 889 | value += fParamRMS0[dim][type][1]*z; |
| 890 | value += fParamRMS0[dim][type][2]*angle*angle; |
| 891 | value = TMath::Sqrt(TMath::Abs(value)); |
| 892 | return value; |
| 893 | } |
| 894 | |
| 895 | Float_t AliTPCClusterParam::GetRMS1(Int_t dim, Int_t type, Float_t z, Float_t angle){ |
| 896 | // |
| 897 | // calculate mean RMS of cluster - z,angle - pad length scalling |
| 898 | // |
| 899 | Float_t value=0; |
| 900 | Float_t length=0.75; |
| 901 | if (type==1) length=1; |
| 902 | if (type==2) length=1.5; |
| 903 | if (type==0){ |
| 904 | value += fParamRMS1[dim][0]; |
| 905 | }else{ |
| 906 | value += fParamRMS1[dim][1]; |
| 907 | } |
| 908 | value += fParamRMS1[dim][2]*z; |
| 909 | value += fParamRMS1[dim][3]*angle*angle*length*length; |
| 910 | value = TMath::Sqrt(TMath::Abs(value)); |
| 911 | return value; |
| 912 | } |
| 913 | |
| 914 | Float_t AliTPCClusterParam::GetRMSQ(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){ |
| 915 | // |
| 916 | // calculate mean RMS of cluster - z,angle, Q dependence |
| 917 | // |
| 918 | Float_t value=0; |
| 919 | value += fParamRMSQ[dim][type][0]; |
| 920 | value += fParamRMSQ[dim][type][1]*z; |
| 921 | value += fParamRMSQ[dim][type][2]*angle*angle; |
| 922 | value += fParamRMSQ[dim][type][3]*z/Qmean; |
| 923 | value += fParamRMSQ[dim][type][4]*angle*angle/Qmean; |
| 924 | value = TMath::Sqrt(TMath::Abs(value)); |
| 925 | return value; |
| 926 | } |
| 927 | |
| 928 | Float_t AliTPCClusterParam::GetRMSSigma(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean){ |
| 929 | // |
| 930 | // calculates RMS of signal shape fluctuation |
| 931 | // |
| 932 | Float_t mean = GetRMSQ(dim,type,z,angle,Qmean); |
| 933 | Float_t value = fRMSSigmaFit[dim][type][0]; |
| 934 | value+= fRMSSigmaFit[dim][type][1]*mean; |
| 935 | return value; |
| 936 | } |
| 937 | |
| 938 | Float_t AliTPCClusterParam::GetShapeFactor(Int_t dim, Int_t type, Float_t z, Float_t angle, Float_t Qmean, Float_t rmsL, Float_t rmsM){ |
| 939 | // |
| 940 | // calculates vallue - sigma distortion contribution |
| 941 | // |
| 942 | Double_t value =0; |
| 943 | // |
| 944 | Float_t rmsMeanQ = GetRMSQ(dim,type,z,angle,Qmean); |
| 945 | if (rmsL<rmsMeanQ) return value; |
| 946 | // |
| 947 | Float_t rmsSigma = GetRMSSigma(dim,type,z,angle,Qmean); |
| 948 | // |
| 949 | if ((rmsM-rmsMeanQ)>2.0*(rmsSigma+fErrorRMSSys[dim])){ |
| 950 | //1.5 sigma cut on mean |
| 951 | value+= rmsL*rmsL+2*rmsM*rmsM-3*rmsMeanQ*rmsMeanQ; |
| 952 | }else{ |
| 953 | if ((rmsL-rmsMeanQ)>3.*(rmsSigma+fErrorRMSSys[dim])){ |
| 954 | //3 sigma cut on local |
| 955 | value+= rmsL*rmsL-rmsMeanQ*rmsMeanQ; |
| 956 | } |
| 957 | } |
| 8076baa0 |
958 | return TMath::Sqrt(TMath::Abs(value)); |
| 12ca5da1 |
959 | } |
| 960 | |
| 961 | |
| 962 | |
| 963 | void AliTPCClusterParam::FitData(TTree * tree){ |
| 964 | // |
| 965 | // make fits for error param and shape param |
| 966 | // |
| 967 | FitResol(tree); |
| 968 | FitRMS(tree); |
| 969 | |
| 970 | } |
| 971 | |
| 972 | void AliTPCClusterParam::FitResol(TTree * tree){ |
| 973 | // |
| 974 | SetInstance(this); |
| 975 | for (Int_t idir=0;idir<2; idir++){ |
| 976 | for (Int_t itype=0; itype<3; itype++){ |
| 977 | Float_t param0[10]; |
| 978 | Float_t error0[10]; |
| 979 | // model error param |
| 980 | FitResol0(tree, idir, itype,param0,error0); |
| 981 | printf("\nResol\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]); |
| 982 | printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]); |
| 983 | printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]); |
| 984 | for (Int_t ipar=0;ipar<4; ipar++){ |
| 985 | fParamS0[idir][itype][ipar] = param0[ipar]; |
| 986 | fErrorS0[idir][itype][ipar] = param0[ipar]; |
| 987 | } |
| 988 | // error param with parabolic correction |
| 989 | FitResol0Par(tree, idir, itype,param0,error0); |
| 990 | printf("\nResolPar\t%d\t%d\tchi2=%f\n",idir,itype,param0[6]); |
| 991 | printf("%f\t%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4],param0[5]); |
| 992 | printf("%f\t%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4],error0[5]); |
| 993 | for (Int_t ipar=0;ipar<7; ipar++){ |
| 994 | fParamS0Par[idir][itype][ipar] = param0[ipar]; |
| 995 | fErrorS0Par[idir][itype][ipar] = param0[ipar]; |
| 996 | } |
| 997 | // |
| 998 | FitResolQ(tree, idir, itype,param0,error0); |
| 999 | printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]); |
| 1000 | printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]); |
| 1001 | printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]); |
| 1002 | for (Int_t ipar=0;ipar<6; ipar++){ |
| 1003 | fParamSQ[idir][itype][ipar] = param0[ipar]; |
| 1004 | fErrorSQ[idir][itype][ipar] = param0[ipar]; |
| 1005 | } |
| 1006 | // |
| 1007 | FitResolQPar(tree, idir, itype,param0,error0); |
| 1008 | printf("\nResolQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[8]); |
| 1009 | printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4],param0[5],param0[6],param0[7]); |
| 1010 | printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4],error0[5],error0[6],error0[7]); |
| 1011 | for (Int_t ipar=0;ipar<9; ipar++){ |
| 1012 | fParamSQPar[idir][itype][ipar] = param0[ipar]; |
| 1013 | fErrorSQPar[idir][itype][ipar] = param0[ipar]; |
| 1014 | } |
| 1015 | } |
| 1016 | } |
| 1017 | // |
| 1018 | printf("Resol z-scaled\n"); |
| 1019 | for (Int_t idir=0;idir<2; idir++){ |
| 1020 | Float_t param0[4]; |
| 1021 | Float_t error0[4]; |
| 1022 | FitResol1(tree, idir,param0,error0); |
| 1023 | printf("\nResol\t%d\tchi2=%f\n",idir,param0[3]); |
| 1024 | printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]); |
| 1025 | printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]); |
| 1026 | for (Int_t ipar=0;ipar<4; ipar++){ |
| 1027 | fParamS1[idir][ipar] = param0[ipar]; |
| 1028 | fErrorS1[idir][ipar] = param0[ipar]; |
| 1029 | } |
| 1030 | } |
| 1031 | |
| 1032 | for (Int_t idir=0;idir<2; idir++){ |
| 1033 | printf("\nDirection %d\n",idir); |
| 1034 | printf("%d\t%f\t%f\t%f\n", -1,fParamS1[idir][0],fParamS1[idir][1],fParamS1[idir][2]); |
| 1035 | for (Int_t itype=0; itype<3; itype++){ |
| 1036 | Float_t length=0.75; |
| 1037 | if (itype==1) length=1; |
| 1038 | if (itype==2) length=1.5; |
| 1039 | printf("%d\t%f\t%f\t%f\n", itype,fParamS0[idir][itype][0],fParamS0[idir][itype][1]*TMath::Sqrt(length),fParamS0[idir][itype][2]/TMath::Sqrt(length)); |
| 1040 | } |
| 1041 | } |
| 1042 | } |
| 1043 | |
| 1044 | |
| 1045 | |
| 1046 | void AliTPCClusterParam::FitRMS(TTree * tree){ |
| 1047 | // |
| 1048 | SetInstance(this); |
| 1049 | for (Int_t idir=0;idir<2; idir++){ |
| 1050 | for (Int_t itype=0; itype<3; itype++){ |
| 1051 | Float_t param0[6]; |
| 1052 | Float_t error0[6]; |
| 1053 | FitRMS0(tree, idir, itype,param0,error0); |
| 1054 | printf("\nRMS\t%d\t%d\tchi2=%f\n",idir,itype,param0[3]); |
| 1055 | printf("%f\t%f\t%f\n", param0[0],param0[1],param0[2]); |
| 1056 | printf("%f\t%f\t%f\n", error0[0],error0[1],error0[2]); |
| 1057 | for (Int_t ipar=0;ipar<4; ipar++){ |
| 1058 | fParamRMS0[idir][itype][ipar] = param0[ipar]; |
| 1059 | fErrorRMS0[idir][itype][ipar] = param0[ipar]; |
| 1060 | } |
| 1061 | FitRMSQ(tree, idir, itype,param0,error0); |
| 1062 | printf("\nRMSQ\t%d\t%d\tchi2=%f\n",idir,itype,param0[5]); |
| 1063 | printf("%f\t%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2],param0[3],param0[4]); |
| 1064 | printf("%f\t%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2],error0[3],error0[4]); |
| 1065 | for (Int_t ipar=0;ipar<6; ipar++){ |
| 1066 | fParamRMSQ[idir][itype][ipar] = param0[ipar]; |
| 1067 | fErrorRMSQ[idir][itype][ipar] = param0[ipar]; |
| 1068 | } |
| 1069 | } |
| 1070 | } |
| 1071 | // |
| 1072 | printf("RMS z-scaled\n"); |
| 1073 | for (Int_t idir=0;idir<2; idir++){ |
| 1074 | Float_t param0[5]; |
| 1075 | Float_t error0[5]; |
| 1076 | FitRMS1(tree, idir,param0,error0); |
| 1077 | printf("\nRMS\t%d\tchi2=%f\n",idir,param0[4]); |
| 1078 | printf("%f\t%f\t%f\t%f\n", param0[0],param0[1],param0[2], param0[3]); |
| 1079 | printf("%f\t%f\t%f\t%f\n", error0[0],error0[1],error0[2], error0[3]); |
| 1080 | for (Int_t ipar=0;ipar<5; ipar++){ |
| 1081 | fParamRMS1[idir][ipar] = param0[ipar]; |
| 1082 | fErrorRMS1[idir][ipar] = param0[ipar]; |
| 1083 | } |
| 1084 | } |
| 1085 | |
| 1086 | for (Int_t idir=0;idir<2; idir++){ |
| 1087 | printf("\nDirection %d\n",idir); |
| 1088 | printf("%d\t%f\t%f\t%f\t%f\n", -1,fParamRMS1[idir][0],fParamRMS1[idir][1],fParamRMS1[idir][2], fParamRMS1[idir][3]); |
| 1089 | for (Int_t itype=0; itype<3; itype++){ |
| 1090 | Float_t length=0.75; |
| 1091 | if (itype==1) length=1; |
| 1092 | if (itype==2) length=1.5; |
| 1093 | if (itype==0) printf("%d\t%f\t\t\t%f\t%f\n", itype,fParamRMS0[idir][itype][0],fParamRMS0[idir][itype][1],fParamRMS0[idir][itype][2]/length); |
| 1094 | if (itype>0) printf("%d\t\t\t%f\t%f\t%f\n", itype,fParamRMS0[idir][itype][0],fParamRMS0[idir][itype][1],fParamRMS0[idir][itype][2]/length); |
| 1095 | } |
| 1096 | } |
| 1097 | // |
| 1098 | // Fit RMS sigma |
| 1099 | // |
| 1100 | printf("RMS fluctuation parameterization \n"); |
| 1101 | for (Int_t idir=0;idir<2; idir++){ |
| 1102 | for (Int_t itype=0; itype<3; itype++){ |
| 1103 | Float_t param0[5]; |
| 1104 | Float_t error0[5]; |
| 1105 | FitRMSSigma(tree, idir,itype,param0,error0); |
| 1106 | printf("\t%d\t%d\t%f\t%f\n", idir, itype, param0[0],param0[1]); |
| 1107 | for (Int_t ipar=0;ipar<2; ipar++){ |
| 1108 | fRMSSigmaFit[idir][itype][ipar] = param0[ipar]; |
| 1109 | } |
| 1110 | } |
| 1111 | } |
| 1112 | // |
| 1113 | // store systematic error end RMS fluctuation parameterization |
| 1114 | // |
| 1115 | TH1F hratio("hratio","hratio",100,-0.1,0.1); |
| 1116 | tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==0&&QMean>0"); |
| 1117 | fErrorRMSSys[0] = hratio.GetRMS(); |
| 1118 | tree->Draw("(RMSm-AliTPCClusterParam::SGetRMSQ(Dim,Pad,Zm,AngleM,QMean))/RMSm>>hratio","Dim==1&&QMean>0"); |
| 1119 | fErrorRMSSys[1] = hratio.GetRMS(); |
| 1120 | TH1F hratioR("hratioR","hratioR",100,0,0.2); |
| 1121 | tree->Draw("RMSs/RMSm>>hratioR","Dim==0&&QMean>0"); |
| 1122 | fRMSSigmaRatio[0][0]=hratioR.GetMean(); |
| 1123 | fRMSSigmaRatio[0][1]=hratioR.GetRMS(); |
| 1124 | tree->Draw("RMSs/RMSm>>hratioR","Dim==1&&QMean>0"); |
| 1125 | fRMSSigmaRatio[1][0]=hratioR.GetMean(); |
| 1126 | fRMSSigmaRatio[1][1]=hratioR.GetRMS(); |
| 1127 | } |
| 1128 | |
| 1129 | void AliTPCClusterParam::Test(TTree * tree, const char *output){ |
| 1130 | // |
| 1131 | // Draw standard quality histograms to output file |
| 1132 | // |
| 1133 | SetInstance(this); |
| 1134 | TFile f(output,"recreate"); |
| 1135 | f.cd(); |
| 1136 | // |
| 1137 | // 1D histograms - resolution |
| 1138 | // |
| 1139 | for (Int_t idim=0; idim<2; idim++){ |
| 1140 | for (Int_t ipad=0; ipad<3; ipad++){ |
| 1141 | char hname1[300]; |
| 1142 | char hcut1[300]; |
| 1143 | char hexp1[300]; |
| 1144 | // |
| 1145 | sprintf(hname1,"Delta0 Dir %d Pad %d",idim,ipad); |
| 1146 | sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad); |
| 1147 | sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1); |
| 1148 | TH1F his1DRel0(hname1, hname1, 100,-0.2, 0.2); |
| 1149 | sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad); |
| 1150 | tree->Draw(hexp1,hcut1,""); |
| 1151 | his1DRel0.Write(); |
| 1152 | // |
| 1153 | sprintf(hname1,"Delta0Par Dir %d Pad %d",idim,ipad); |
| 1154 | sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad); |
| 1155 | sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol>>%s",hname1); |
| 1156 | TH1F his1DRel0Par(hname1, hname1, 100,-0.2, 0.2); |
| 1157 | sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad); |
| 1158 | tree->Draw(hexp1,hcut1,""); |
| 1159 | his1DRel0Par.Write(); |
| 1160 | // |
| 1161 | } |
| 1162 | } |
| 1163 | // |
| 1164 | // 2D histograms - resolution |
| 1165 | // |
| 1166 | for (Int_t idim=0; idim<2; idim++){ |
| 1167 | for (Int_t ipad=0; ipad<3; ipad++){ |
| 1168 | char hname1[300]; |
| 1169 | char hcut1[300]; |
| 1170 | char hexp1[300]; |
| 1171 | // |
| 1172 | sprintf(hname1,"2DDelta0 Dir %d Pad %d",idim,ipad); |
| 1173 | sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad); |
| 1174 | sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1); |
| 1175 | TProfile2D profDRel0(hname1, hname1, 6,0,250,6,0,1); |
| 1176 | sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad); |
| 1177 | tree->Draw(hexp1,hcut1,""); |
| 1178 | profDRel0.Write(); |
| 1179 | // |
| 1180 | sprintf(hname1,"2DDelta0Par Dir %d Pad %d",idim,ipad); |
| 1181 | sprintf(hcut1,"Dim==%d&&QMean<0&&Pad==%d",idim,ipad); |
| 1182 | sprintf(hexp1,"(Resol-AliTPCClusterParam::SGetError0Par(Dim,Pad,Zm,AngleM))/Resol:AngleM:Zm>>%s",hname1); |
| 1183 | TProfile2D profDRel0Par(hname1, hname1,6,0,250,6,0,1); |
| 1184 | sprintf(hname1,"Dim==%d&&QMean<0&&Pad=%d",idim,ipad); |
| 1185 | tree->Draw(hexp1,hcut1,""); |
| 1186 | profDRel0Par.Write(); |
| 1187 | // |
| 1188 | } |
| 1189 | } |
| 1190 | } |
| 1191 | |
| 1192 | |
| 1193 | |
| 1194 | void AliTPCClusterParam::Print(Option_t* /*option*/) const{ |
| 1195 | // |
| 1196 | // Print param Information |
| 1197 | // |
| 1198 | |
| 1199 | // |
| 1200 | // Error parameterization |
| 1201 | // |
| 1202 | printf("\nResolution Scaled factors\n"); |
| 1203 | printf("Dir\tPad\tP0\t\tP1\t\tP2\t\tchi2\n"); |
| 8076baa0 |
1204 | printf("Y\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[0][0])),TMath::Sqrt(TMath::Abs(fParamS1[0][1])), |
| 1205 | TMath::Sqrt(TMath::Abs(fParamS1[0][2])),TMath::Sqrt(TMath::Abs(fParamS1[0][3]))); |
| 12ca5da1 |
1206 | for (Int_t ipad=0; ipad<3; ipad++){ |
| 1207 | Float_t length=0.75; |
| 1208 | if (ipad==1) length=1; |
| 1209 | if (ipad==2) length=1.5; |
| 1210 | printf("\t%d\t%f\t%f\t%f\t%f\n", ipad, |
| 1211 | TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][0])), |
| 8076baa0 |
1212 | TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][1]*length)), |
| 1213 | TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][2]/length)), |
| 1214 | TMath::Sqrt(TMath::Abs(fParamS0[0][ipad][3]))); |
| 12ca5da1 |
1215 | } |
| 1216 | for (Int_t ipad=0; ipad<3; ipad++){ |
| 1217 | Float_t length=0.75; |
| 1218 | if (ipad==1) length=1; |
| 1219 | if (ipad==2) length=1.5; |
| 1220 | printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad, |
| 1221 | TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][0])), |
| 8076baa0 |
1222 | TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][1]*length)), |
| 1223 | TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][2]/length)), |
| 1224 | TMath::Sqrt(TMath::Abs(fParamS0Par[0][ipad][6]))); |
| 12ca5da1 |
1225 | } |
| 1226 | printf("Z\tall\t%f\t%f\t%f\t%f\n", TMath::Sqrt(TMath::Abs(fParamS1[1][0])),TMath::Sqrt(fParamS1[1][1]), |
| 1227 | TMath::Sqrt(fParamS1[1][2]), TMath::Sqrt(fParamS1[1][3])); |
| 1228 | |
| 1229 | for (Int_t ipad=0; ipad<3; ipad++){ |
| 1230 | Float_t length=0.75; |
| 1231 | if (ipad==1) length=1; |
| 1232 | if (ipad==2) length=1.5; |
| 1233 | printf("\t%d\t%f\t%f\t%f\t%f\n", ipad, |
| 1234 | TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][0])), |
| 8076baa0 |
1235 | TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][1]*length)), |
| 1236 | TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][2]/length)), |
| 1237 | TMath::Sqrt(TMath::Abs(fParamS0[1][ipad][3]))); |
| 12ca5da1 |
1238 | } |
| 1239 | for (Int_t ipad=0; ipad<3; ipad++){ |
| 1240 | Float_t length=0.75; |
| 1241 | if (ipad==1) length=1; |
| 1242 | if (ipad==2) length=1.5; |
| 1243 | printf("\t%dPar\t%f\t%f\t%f\t%f\n", ipad, |
| 8076baa0 |
1244 | TMath::Sqrt(TMath::Abs(TMath::Abs(fParamS0Par[1][ipad][0]))), |
| 1245 | TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][1]*length)), |
| 1246 | TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][2]/length)), |
| 1247 | TMath::Sqrt(TMath::Abs(fParamS0Par[1][ipad][6]))); |
| 12ca5da1 |
1248 | } |
| 1249 | |
| 1250 | // |
| 1251 | // RMS scaling |
| 1252 | // |
| 1253 | printf("\n"); |
| 1254 | printf("\nRMS Scaled factors\n"); |
| 1255 | printf("Dir\tPad\tP00\t\tP01\t\tP1\t\tP2\t\tchi2\n"); |
| 8076baa0 |
1256 | printf("Y\tall\t%f\t%f\t%f\t%f\t%f\n", |
| 1257 | TMath::Sqrt(TMath::Abs(fParamRMS1[0][0])), |
| 1258 | TMath::Sqrt(TMath::Abs(fParamRMS1[0][1])), |
| 1259 | TMath::Sqrt(TMath::Abs(fParamRMS1[0][2])), |
| 1260 | TMath::Sqrt(TMath::Abs(fParamRMS1[0][3])), |
| 1261 | TMath::Sqrt(TMath::Abs(fParamRMS1[0][4]))); |
| 12ca5da1 |
1262 | for (Int_t ipad=0; ipad<3; ipad++){ |
| 1263 | Float_t length=0.75; |
| 1264 | if (ipad==1) length=1; |
| 1265 | if (ipad==2) length=1.5; |
| 1266 | if (ipad==0){ |
| 1267 | printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad, |
| 1268 | TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])), |
| 1269 | 0., |
| 8076baa0 |
1270 | TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])), |
| 1271 | TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))), |
| 1272 | TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3]))); |
| 12ca5da1 |
1273 | }else{ |
| 1274 | printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad, |
| 1275 | 0., |
| 1276 | TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][0])), |
| 8076baa0 |
1277 | TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][1])), |
| 1278 | TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][2]/(length*length))), |
| 1279 | TMath::Sqrt(TMath::Abs(fParamRMS0[0][ipad][3]))); |
| 12ca5da1 |
1280 | } |
| 1281 | } |
| 1282 | printf("\n"); |
| 8076baa0 |
1283 | printf("Z\tall\t%f\t%f\t%f\t%f\t%f\n", |
| 1284 | TMath::Sqrt(TMath::Abs(fParamRMS1[1][0])), |
| 1285 | TMath::Sqrt(TMath::Abs(fParamRMS1[1][1])), |
| 1286 | TMath::Sqrt(TMath::Abs(fParamRMS1[1][2])), |
| 1287 | TMath::Sqrt(TMath::Abs(fParamRMS1[1][3])), |
| 1288 | TMath::Sqrt(TMath::Abs(fParamRMS1[1][4]))); |
| 12ca5da1 |
1289 | for (Int_t ipad=0; ipad<3; ipad++){ |
| 1290 | Float_t length=0.75; |
| 1291 | if (ipad==1) length=1; |
| 1292 | if (ipad==2) length=1.5; |
| 1293 | if (ipad==0){ |
| 1294 | printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad, |
| 1295 | TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])), |
| 1296 | 0., |
| 8076baa0 |
1297 | TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])), |
| 1298 | TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))), |
| 1299 | TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3]))); |
| 12ca5da1 |
1300 | }else{ |
| 1301 | printf("\t%d\t%f\t%f\t%f\t%f\t%f\n", ipad, |
| 1302 | 0., |
| 1303 | TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][0])), |
| 8076baa0 |
1304 | TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][1])), |
| 1305 | TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][2]/(length*length))), |
| 1306 | TMath::Sqrt(TMath::Abs(fParamRMS0[1][ipad][3]))); |
| 12ca5da1 |
1307 | } |
| 1308 | } |
| 1309 | } |
| 1310 | |
| 1311 | |
| 1312 | |
| 1313 | |
| 1314 | |
| 0a65832b |
1315 | Float_t AliTPCClusterParam::Qnorm(Int_t ipad, Int_t itype, Float_t dr, Float_t ty, Float_t tz){ |
| 1316 | // get Q normalization |
| 1317 | // type - 0 Qtot 1 Qmax |
| 1318 | // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm) |
| 1319 | // |
| f1afff3b |
1320 | //expession formula - TString *strq0 = toolkit.FitPlane(chain,"dedxQ.fElements[2]","dr++ty++tz++dr*ty++dr*tz++ty*tz++ty^2++tz^2","IPad==0",chi2,npoints,param,covar,0,100000); |
| 1321 | |
| f1c2a4a3 |
1322 | if (fQNorm==0) return 0; |
| 0a65832b |
1323 | TVectorD * norm = (TVectorD*)fQNorm->At(3*itype+ipad); |
| 1324 | if (!norm) return 0; |
| f1afff3b |
1325 | TVectorD &no = *norm; |
| 684602c8 |
1326 | Float_t res = |
| 1327 | no[0]+ |
| f1afff3b |
1328 | no[1]*dr+ |
| 1329 | no[2]*ty+ |
| 1330 | no[3]*tz+ |
| 1331 | no[4]*dr*ty+ |
| 1332 | no[5]*dr*tz+ |
| 1333 | no[6]*ty*tz+ |
| 1334 | no[7]*dr*dr+ |
| 1335 | no[8]*ty*ty+ |
| 1336 | no[9]*tz*tz; |
| 1337 | res/=no[0]; |
| 0a65832b |
1338 | return res; |
| 1339 | } |
| 1340 | |
| 1341 | |
| 1342 | |
| 1343 | void AliTPCClusterParam::SetQnorm(Int_t ipad, Int_t itype, TVectorD * norm){ |
| 1344 | // |
| 1345 | // set normalization |
| 1346 | // |
| 1347 | // type - 0 Qtot 1 Qmax |
| 1348 | // ipad - 0 (0.75 cm) ,1 (1 cm), 2 (1.5 cm) |
| 1349 | // |
| 1350 | |
| 1351 | if (fQNorm==0) fQNorm = new TObjArray(6); |
| 1352 | fQNorm->AddAt(new TVectorD(*norm), itype*3+ipad); |
| 1353 | } |
| 236a0d03 |
1354 | |
| b17540e4 |
1355 | Float_t AliTPCClusterParam::QnormPos(Int_t ipad,Bool_t isMax, Float_t pad, Float_t time, Float_t z, Float_t sy2, Float_t sz2, Float_t qm, Float_t qt){ |
| 1356 | // |
| 1357 | // Make Q normalization as function of following parameters |
| 1358 | // Fit parameters to be used in corresponding correction function extracted in the AliTPCclaibTracksGain - Taylor expansion |
| 1359 | // 1 - dp - relative pad position |
| 1360 | // 2 - dt - relative time position |
| 1361 | // 3 - di - drift length (norm to 1); |
| 1362 | // 4 - dq0 - Tot/Max charge |
| 1363 | // 5 - dq1 - Max/Tot charge |
| 1364 | // 6 - sy - sigma y - shape |
| 1365 | // 7 - sz - sigma z - shape |
| 1366 | // |
| 1367 | |
| 1368 | //The results can be visualized using the debug streamer information of the AliTPCcalibTracksGain - |
| 1369 | // Following variable used - correspondance to the our variable conventions |
| 1370 | //chain0->SetAlias("dp","((Cl.fPad-int(Cl.fPad)-0.5)/0.5)"); |
| 1371 | Double_t dp = ((pad-int(pad)-0.5)*2.); |
| 1372 | //chain0->SetAlias("dt","((Cl.fTimeBin-int(Cl.fTimeBin)-0.5)/0.5)"); |
| 1373 | Double_t dt = ((time-int(time)-0.5)*2.); |
| 1374 | //chain0->SetAlias("di","(sqrt(1.-abs(Cl.fZ)/250.))"); |
| 1375 | Double_t di = TMath::Sqrt(1-TMath::Abs(z)/250.); |
| 1376 | //chain0->SetAlias("dq0","(0.2*(Cl.fQ+2)/(Cl.fMax+2))"); |
| 1377 | Double_t dq0 = 0.2*(qt+2.)/(qm+2.); |
| 1378 | //chain0->SetAlias("dq1","(5*(Cl.fMax+2)/(Cl.fQ+2))"); |
| 1379 | Double_t dq1 = 5.*(qm+2.)/(qt+2.); |
| 1380 | //chain0->SetAlias("sy","(0.32/sqrt(0.01^2+Cl.fSigmaY2))"); |
| 1381 | Double_t sy = 0.32/TMath::Sqrt(0.01*0.01+sy2); |
| 1382 | //chain0->SetAlias("sz","(0.32/sqrt(0.01^2+Cl.fSigmaZ2))"); |
| 1383 | Double_t sz = 0.32/TMath::Sqrt(0.01*0.01+sz2); |
| 1384 | // |
| 1385 | // |
| 1386 | // |
| 1387 | TVectorD * pvec = 0; |
| 1388 | if (isMax){ |
| 1389 | pvec = fPosQMnorm[ipad]; |
| 1390 | }else{ |
| 1391 | pvec = fPosQTnorm[ipad]; |
| 1392 | } |
| 1393 | TVectorD ¶m = *pvec; |
| 1394 | // |
| 1395 | // Eval part - in correspondance with fit part from debug streamer |
| 1396 | // |
| 1397 | Double_t result=param[0]; |
| 1398 | Int_t index =1; |
| 1399 | // |
| 1400 | result+=dp*param[index++]; //1 |
| 1401 | result+=dt*param[index++]; //2 |
| 1402 | result+=dp*dp*param[index++]; //3 |
| 1403 | result+=dt*dt*param[index++]; //4 |
| 1404 | result+=dt*dt*dt*param[index++]; //5 |
| 1405 | result+=dp*dt*param[index++]; //6 |
| 1406 | result+=dp*dt*dt*param[index++]; //7 |
| 1407 | result+=(dq0)*param[index++]; //8 |
| 1408 | result+=(dq1)*param[index++]; //9 |
| 1409 | // |
| 1410 | // |
| 1411 | result+=dp*dp*(di)*param[index++]; //10 |
| 1412 | result+=dt*dt*(di)*param[index++]; //11 |
| 1413 | result+=dp*dp*sy*param[index++]; //12 |
| 1414 | result+=dt*sz*param[index++]; //13 |
| 1415 | result+=dt*dt*sz*param[index++]; //14 |
| 1416 | result+=dt*dt*dt*sz*param[index++]; //15 |
| 1417 | // |
| 1418 | result+=dp*dp*1*sy*sz*param[index++]; //16 |
| 1419 | result+=dt*sy*sz*param[index++]; //17 |
| 1420 | result+=dt*dt*sy*sz*param[index++]; //18 |
| 1421 | result+=dt*dt*dt*sy*sz*param[index++]; //19 |
| 1422 | // |
| 1423 | result+=dp*dp*(dq0)*param[index++]; //20 |
| 1424 | result+=dt*1*(dq0)*param[index++]; //21 |
| 1425 | result+=dt*dt*(dq0)*param[index++]; //22 |
| 1426 | result+=dt*dt*dt*(dq0)*param[index++]; //23 |
| 1427 | // |
| 1428 | result+=dp*dp*(dq1)*param[index++]; //24 |
| 1429 | result+=dt*(dq1)*param[index++]; //25 |
| 1430 | result+=dt*dt*(dq1)*param[index++]; //26 |
| 1431 | result+=dt*dt*dt*(dq1)*param[index++]; //27 |
| 1432 | |
| 2e5bcb67 |
1433 | if (result<0.75) result=0.75; |
| 1434 | if (result>1.25) result=1.25; |
| 1435 | |
| b17540e4 |
1436 | return result; |
| 1437 | |
| 1438 | } |
| 236a0d03 |
1439 | |
| 1440 | |
| 1441 | |
| 236a0d03 |
1442 | |
| 236a0d03 |
1443 | |
| bf97e1c4 |
1444 | Float_t AliTPCClusterParam::PosCorrection(Int_t type, Int_t ipad, Float_t pad, Float_t time, Float_t z, Float_t /*sy2*/, Float_t /*sz2*/, Float_t /*qm*/){ |
| 2e5bcb67 |
1445 | |
| 1446 | // |
| 1447 | // Make postion correction |
| 1448 | // type - 0 - y correction |
| 1449 | // 1 - z correction |
| 1450 | // ipad - 0, 1, 2 - short, medium long pads |
| 1451 | // pad - float pad number |
| 1452 | // time - float time bin number |
| 1453 | // z - z of the cluster |
| 2e5bcb67 |
1454 | |
| 1455 | // |
| 1456 | //chainres->SetAlias("dp","(-1+(Cl.fZ>0)*2)*((Cl.fPad-int(Cl.fPad))-0.5)"); |
| 1457 | //chainres->SetAlias("dt","(-1+(Cl.fZ>0)*2)*((Cl.fTimeBin-0.66-int(Cl.fTimeBin-0.66))-0.5)"); |
| 1458 | //chainres->SetAlias("sp","(sin(dp*pi)-dp*pi)"); |
| 1459 | //chainres->SetAlias("st","(sin(dt)-dt)"); |
| 1460 | // |
| 1461 | //chainres->SetAlias("di","sqrt(1.-abs(Cl.fZ/250.))"); |
| 2e5bcb67 |
1462 | |
| 1463 | // |
| 1464 | // Derived variables |
| 1465 | // |
| 1466 | Double_t dp = (-1+(z>0)*2)*((pad-int(pad))-0.5); |
| 1467 | Double_t dt = (-1+(z>0)*2)*((time-0.66-int(time-0.66))-0.5); |
| 1468 | Double_t sp = (TMath::Sin(dp*TMath::Pi())-dp*TMath::Pi()); |
| 1469 | Double_t st = (TMath::Sin(dt)-dt); |
| 1470 | // |
| bf97e1c4 |
1471 | Double_t di = TMath::Sqrt(TMath::Abs(1.-TMath::Abs(z/250.))); |
| 2e5bcb67 |
1472 | // |
| 1473 | // |
| 1474 | // |
| 1475 | TVectorD * pvec = 0; |
| 1476 | if (type==0){ |
| 1477 | pvec = fPosYcor[ipad]; |
| 1478 | }else{ |
| 1479 | pvec = fPosZcor[ipad]; |
| 1480 | } |
| 1481 | TVectorD ¶m = *pvec; |
| 1482 | // |
| bf97e1c4 |
1483 | Double_t result=0; |
| 2e5bcb67 |
1484 | Int_t index =1; |
| 1485 | |
| 1486 | if (type==0){ |
| 1487 | // y corr |
| 1488 | result+=(dp)*param[index++]; //1 |
| 1489 | result+=(dp)*di*param[index++]; //2 |
| 2e5bcb67 |
1490 | // |
| bf97e1c4 |
1491 | result+=(sp)*param[index++]; //3 |
| 1492 | result+=(sp)*di*param[index++]; //4 |
| 2e5bcb67 |
1493 | } |
| 1494 | if (type==1){ |
| 1495 | result+=(dt)*param[index++]; //1 |
| 1496 | result+=(dt)*di*param[index++]; //2 |
| 2e5bcb67 |
1497 | // |
| bf97e1c4 |
1498 | result+=(st)*param[index++]; //3 |
| 1499 | result+=(st)*di*param[index++]; //4 |
| 2e5bcb67 |
1500 | } |
| bf97e1c4 |
1501 | if (TMath::Abs(result)>0.05) return 0; |
| 2e5bcb67 |
1502 | return result; |
| 1503 | } |
| 1504 | |
| 1505 | |
| 1506 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame