| 93096a07 |
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 | Responsible: marian.ivanov@cern.ch |
| 20 | Tools for fitting of the space point distortion parameters. |
| 21 | Functionality |
| 22 | |
| 23 | |
| 24 | 1. Creation of the distortion maps from the residual histograms |
| 25 | 2. Making fit trees |
| 26 | 3. Parameters to calculate- dO/dp |
| 27 | |
| 28 | Some functions, for the moment function present in the AliTPCPreprocesorOffline, some will be |
| 29 | extracted from the old macros |
| 30 | |
| 31 | |
| 32 | */ |
| 33 | |
| 34 | #include "Riostream.h" |
| 35 | #include <fstream> |
| 36 | #include "TMap.h" |
| 37 | #include "TGraphErrors.h" |
| 38 | #include "AliExternalTrackParam.h" |
| 39 | #include "TFile.h" |
| 40 | #include "TGraph.h" |
| 41 | #include "TMultiGraph.h" |
| 42 | #include "TCanvas.h" |
| 43 | #include "THnSparse.h" |
| 44 | #include "THnBase.h" |
| 45 | #include "TProfile.h" |
| 46 | #include "TROOT.h" |
| 47 | #include "TLegend.h" |
| 48 | #include "TPad.h" |
| 49 | #include "TH2D.h" |
| 50 | #include "TH3D.h" |
| 51 | #include "AliTPCROC.h" |
| 52 | #include "AliTPCCalROC.h" |
| 53 | #include "AliESDfriend.h" |
| 54 | #include "AliTPCcalibTime.h" |
| 55 | #include "AliSplineFit.h" |
| 56 | #include "AliCDBMetaData.h" |
| 57 | #include "AliCDBId.h" |
| 58 | #include "AliCDBManager.h" |
| 59 | #include "AliCDBStorage.h" |
| 60 | #include "AliTPCcalibBase.h" |
| 61 | #include "AliTPCcalibDB.h" |
| 62 | #include "AliTPCcalibDButil.h" |
| 63 | #include "AliRelAlignerKalman.h" |
| 64 | #include "AliTPCParamSR.h" |
| 65 | #include "AliTPCcalibTimeGain.h" |
| 66 | #include "AliTPCcalibGainMult.h" |
| 67 | #include "AliSplineFit.h" |
| 68 | #include "AliTPCComposedCorrection.h" |
| 69 | #include "AliTPCExBTwist.h" |
| 70 | #include "AliTPCCalibGlobalMisalignment.h" |
| 71 | #include "TStatToolkit.h" |
| 72 | #include "TChain.h" |
| 73 | #include "TCut.h" |
| 74 | #include "AliTrackerBase.h" |
| 75 | #include "AliTPCCorrectionFit.h" |
| 76 | #include "AliTPCLaserTrack.h" |
| 77 | #include "TDatabasePDG.h" |
| 78 | #include "AliTPCcalibAlign.h" |
| 79 | |
| 80 | ClassImp(AliTPCCorrectionFit) |
| 81 | |
| 82 | AliTPCCorrectionFit::AliTPCCorrectionFit(): |
| 83 | TNamed("TPCCorrectionFit","TPCCorrectionFit") |
| 84 | { |
| 85 | // |
| 86 | // default constructor |
| 87 | // |
| 88 | } |
| 89 | |
| 90 | AliTPCCorrectionFit::~AliTPCCorrectionFit() { |
| 91 | // |
| 92 | // Destructor |
| 93 | // |
| 94 | } |
| 95 | |
| 96 | |
| 97 | Double_t AliTPCCorrectionFit::EvalAt(Double_t phi, Double_t refX, Double_t theta, Int_t corr, Int_t ptype){ |
| 98 | // |
| 99 | // |
| 100 | // |
| 101 | Double_t sector = 9*phi/TMath::Pi(); |
| 102 | if (sector<0) sector+=18; |
| 103 | Double_t y85=AliTPCCorrection::GetCorrSector(sector,85,theta,1,corr); |
| 104 | Double_t y245=AliTPCCorrection::GetCorrSector(sector,245,theta,1,corr); |
| 105 | if (ptype==0) return y85+(y245-y85)*(refX-85.)/(245.-85.); |
| 106 | if (ptype==2) return (y245-y85)/(245.-85.); |
| 107 | return 0; |
| 108 | } |
| 109 | |
| 110 | |
| 111 | |
| 112 | |
| 113 | Double_t AliTPCCorrectionFit::EvalAtPar(Double_t phi0, Double_t snp, Double_t refX, Double_t theta, Int_t corr, Int_t ptype, Int_t nsteps){ |
| 114 | // |
| 115 | // Fit the distortion along the line with the parabolic model |
| 116 | // Parameters: |
| 117 | // phi0 - phi at the entrance of the TPC |
| 118 | // snp - local inclination angle at the entrance of the TPC |
| 119 | // refX - ref X where the distortion is evanluated |
| 120 | // theta |
| 121 | // |
| 122 | static TLinearFitter fitter(3,"pol2"); |
| 123 | fitter.ClearPoints(); |
| 124 | if (nsteps<3) nsteps=3; |
| 125 | Double_t deltaX=(245-85)/(nsteps); |
| 126 | for (Int_t istep=0; istep<(nsteps+1); istep++){ |
| 127 | // |
| 128 | Double_t localX =85.+deltaX*istep; |
| 129 | Double_t localPhi=phi0+deltaX*snp*istep; |
| 130 | Double_t sector = 9*localPhi/TMath::Pi(); |
| 131 | if (sector<0) sector+=18; |
| 132 | Double_t y=AliTPCCorrection::GetCorrSector(sector,localX,theta,1,corr); |
| 133 | Double_t dlocalX=AliTPCCorrection::GetCorrSector(sector,localX,theta,0,corr); |
| 134 | Double_t x[1]={localX-dlocalX}; |
| 135 | fitter.AddPoint(x,y); |
| 136 | } |
| 137 | fitter.Eval(); |
| 138 | Double_t par[3]; |
| 139 | par[0]=fitter.GetParameter(0); |
| 140 | par[1]=fitter.GetParameter(1); |
| 141 | par[2]=fitter.GetParameter(2); |
| 142 | |
| 143 | if (ptype==0) return par[0]+par[1]*refX+par[2]*refX*refX; |
| 144 | if (ptype==2) return par[1]+2*par[2]*refX; |
| 145 | if (ptype==4) return par[2]; |
| 146 | return 0; |
| 147 | } |
| 148 | |
| 149 | |
| 150 | |
| 151 | |
| 152 | void AliTPCCorrectionFit::CreateAlignMaps(Double_t bz, Int_t run){ |
| 153 | // |
| 154 | // Create cluster distortion map |
| 155 | // |
| 156 | TFile *falign = TFile::Open("CalibObjects.root"); |
| 157 | TObjArray * arrayAlign = (TObjArray *)falign->Get("TPCAlign"); |
| 158 | AliTPCcalibAlign * align = (AliTPCcalibAlign *)arrayAlign->FindObject("alignTPC"); |
| 159 | TTreeSRedirector * pcstream = new TTreeSRedirector("TPCAlign.root"); |
| 160 | |
| 161 | THnBase * hdY = (THnBase*)align->GetClusterDelta(0); |
| 162 | //THnBase * hdZ = (THnBase*)align->GetClusterDelta(1); |
| 163 | AliTPCCorrectionFit::MakeClusterDistortionMap(hdY,pcstream,0, bz); |
| 164 | // AliTPCCorrectionFit::MakeClusterDistortionMap(hdZ,pcstream,1, bz); |
| 165 | |
| 166 | const char * hname[5]={"dy","dz","dsnp","dtheta","d1pt"}; |
| 167 | for (Int_t ihis=0; ihis<4; ihis++){ |
| 168 | THnSparse * hisAlign =align->GetTrackletDelta(ihis); |
| 169 | AliTPCCorrectionFit::MakeDistortionMapSector(hisAlign, pcstream, hname[ihis], run, ihis,bz); |
| 170 | } |
| 171 | delete pcstream; |
| 172 | delete falign; |
| 173 | } |
| 174 | |
| 175 | |
| 176 | |
| 177 | |
| 178 | |
| 179 | void AliTPCCorrectionFit::MakeClusterDistortionMap(THnBase * hisInput,TTreeSRedirector *pcstream , Int_t ptype, Int_t dtype){ |
| 180 | // |
| 181 | // Make cluster residual map from the n-dimensional histogram |
| 182 | // hisInput supposed to have given format: |
| 183 | // 4 Dim: |
| 184 | // delta, |
| 185 | // sector |
| 186 | // localX |
| 187 | // kZ |
| 188 | // Vertex position assumed to be at (0,0,0) |
| 189 | // |
| 190 | //TTreeSRedirector *pcstream=new TTreeSRedirector(sname); |
| 191 | // |
| 192 | Int_t nbins1=hisInput->GetAxis(1)->GetNbins(); |
| 193 | Int_t nbins2=hisInput->GetAxis(2)->GetNbins(); |
| 194 | Int_t nbins3=hisInput->GetAxis(3)->GetNbins(); |
| 195 | TF1 *fgaus=0; |
| 196 | TH3F * hisResMap3D = |
| 197 | new TH3F("his3D","his3D", |
| 198 | nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(), |
| 199 | nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax(), |
| 200 | nbins3,hisInput->GetAxis(3)->GetXmin(), hisInput->GetAxis(3)->GetXmax()); |
| 201 | hisResMap3D->GetXaxis()->SetTitle("sector"); |
| 202 | hisResMap3D->GetYaxis()->SetTitle("localX"); |
| 203 | hisResMap3D->GetZaxis()->SetTitle("kZ"); |
| 204 | |
| 205 | TH2F * hisResMap2D[4] ={0,0,0,0}; |
| 206 | for (Int_t i=0; i<4; i++){ |
| 207 | hisResMap2D[i]= |
| 208 | new TH2F(Form("his2D_0%d",i),Form("his2D_0%d",i), |
| 209 | nbins1,hisInput->GetAxis(1)->GetXmin(), hisInput->GetAxis(1)->GetXmax(), |
| 210 | nbins2,hisInput->GetAxis(2)->GetXmin(), hisInput->GetAxis(2)->GetXmax()); |
| 211 | hisResMap2D[i]->GetXaxis()->SetTitle("sector"); |
| 212 | hisResMap2D[i]->GetYaxis()->SetTitle("localX"); |
| 213 | } |
| 214 | // |
| 215 | // |
| 216 | // |
| 217 | TF1 * f1= 0; |
| 218 | Int_t axis0[4]={0,1,2,3}; |
| 219 | Int_t axis1[4]={0,1,2,3}; |
| 220 | Int_t counter=0; |
| 221 | for (Int_t ibin1=1; ibin1<nbins1; ibin1+=1){ |
| 222 | // phi- sector range |
| 223 | hisInput->GetAxis(1)->SetRange(ibin1-1,ibin1+1); |
| 224 | THnBase *his1=(THnBase *)hisInput->ProjectionND(4,axis0); |
| 225 | Double_t sector=hisInput->GetAxis(1)->GetBinCenter(ibin1); |
| 226 | // |
| 227 | for (Int_t ibin2=1; ibin2<nbins2; ibin2+=1){ |
| 228 | // local x range |
| 229 | // kz fits |
| 230 | his1->GetAxis(2)->SetRange(ibin2-1,ibin2+1); |
| 231 | THnBase *his2=(THnBase *)his1->ProjectionND(4,axis1); |
| 232 | Double_t localX=hisInput->GetAxis(2)->GetBinCenter(ibin2); |
| 233 | // |
| 234 | //A side |
| 235 | his2->GetAxis(3)->SetRangeUser(0.01,0.3); |
| 236 | TH1 * hisA = his2->Projection(0); |
| 237 | Double_t meanA= hisA->GetMean(); |
| 238 | Double_t rmsA= hisA->GetRMS(); |
| 239 | Double_t entriesA= hisA->GetEntries(); |
| 240 | delete hisA; |
| 241 | //C side |
| 242 | his2->GetAxis(3)->SetRangeUser(0.01,0.3); |
| 243 | TH1 * hisC = his2->Projection(0); |
| 244 | Double_t meanC= hisC->GetMean(); |
| 245 | Double_t rmsC= hisC->GetRMS(); |
| 246 | Double_t entriesC= hisC->GetEntries(); |
| 247 | delete hisC; |
| 248 | his2->GetAxis(3)->SetRangeUser(-1.2,1.2); |
| 249 | TH2 * hisAC = his2->Projection(0,3); |
| 250 | TProfile *profAC = hisAC->ProfileX(); |
| 251 | delete hisAC; |
| 252 | profAC->Fit("pol1","QNR","QNR",0.05,1); |
| 253 | if (!f1) f1=(TF1*)gROOT->FindObject("pol1"); |
| 254 | Double_t offsetA=f1->GetParameter(0); |
| 255 | Double_t slopeA=f1->GetParameter(1); |
| 256 | Double_t offsetAE=f1->GetParError(0); |
| 257 | Double_t slopeAE=f1->GetParError(1); |
| 258 | Double_t chi2A=f1->GetChisquare()/f1->GetNumberFreeParameters(); |
| 259 | profAC->Fit("pol1","QNR","QNR",-1.1,-0.1); |
| 260 | f1=(TF1*)gROOT->FindObject("pol1"); |
| 261 | Double_t offsetC=f1->GetParameter(0); |
| 262 | Double_t slopeC=f1->GetParameter(1); |
| 263 | Double_t offsetCE=f1->GetParError(0); |
| 264 | Double_t slopeCE=f1->GetParError(1); |
| 265 | Double_t chi2C=f1->GetChisquare()/f1->GetNumberFreeParameters(); |
| 266 | if (counter%50==0) printf("%f\t%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX, entriesA+entriesC, slopeA,slopeC, chi2A, chi2C); |
| 267 | counter++; |
| 268 | (*pcstream)<<"deltaFit"<< |
| 269 | "sector="<<sector<< |
| 270 | "localX="<<localX<< |
| 271 | "meanA="<<meanA<< |
| 272 | "rmsA="<<rmsA<< |
| 273 | "entriesA="<<entriesA<< |
| 274 | "meanC="<<meanC<< |
| 275 | "rmsC="<<rmsC<< |
| 276 | "entriesC="<<entriesC<< |
| 277 | "offsetA="<<offsetA<< |
| 278 | "slopeA="<<slopeA<< |
| 279 | "offsetAE="<<offsetAE<< |
| 280 | "slopeAE="<<slopeAE<< |
| 281 | "chi2A="<<chi2A<< |
| 282 | "offsetC="<<offsetC<< |
| 283 | "slopeC="<<slopeC<< |
| 284 | "offsetCE="<<offsetCE<< |
| 285 | "slopeCE="<<slopeCE<< |
| 286 | "chi2C="<<chi2C<< |
| 287 | "\n"; |
| 288 | // |
| 289 | hisResMap2D[0]->SetBinContent(ibin1,ibin2, offsetA); |
| 290 | hisResMap2D[1]->SetBinContent(ibin1,ibin2, slopeA); |
| 291 | hisResMap2D[2]->SetBinContent(ibin1,ibin2, offsetC); |
| 292 | hisResMap2D[3]->SetBinContent(ibin1,ibin2, slopeC); |
| 293 | |
| 294 | for (Int_t ibin3=1; ibin3<nbins3; ibin3++){ |
| 295 | Double_t kZ=hisInput->GetAxis(3)->GetBinCenter(ibin3); |
| 296 | if (TMath::Abs(kZ)<0.05) continue; // crossing |
| 297 | his2->GetAxis(3)->SetRange(ibin3,ibin3); |
| 298 | if (TMath::Abs(kZ)>0.15){ |
| 299 | his2->GetAxis(3)->SetRange(ibin3,ibin3); |
| 300 | } |
| 301 | TH1 * his = his2->Projection(0); |
| 302 | Double_t mean= his->GetMean(); |
| 303 | Double_t rms= his->GetRMS(); |
| 304 | Double_t entries= his->GetEntries(); |
| 305 | //printf("%f\t%f\t%f\t%f\t%f\t%f\n", sector,localX,kZ, entries, mean,rms); |
| 306 | hisResMap3D->SetBinContent(ibin1,ibin2,ibin3, mean); |
| 307 | Double_t phi=TMath::Pi()*sector/9; |
| 308 | if (phi>TMath::Pi()) phi+=TMath::Pi(); |
| 309 | Double_t meanG=0; |
| 310 | Double_t rmsG=0; |
| 311 | if (entries>50){ |
| 312 | if (!fgaus) { |
| 313 | his->Fit("gaus","Q","goff"); |
| 314 | fgaus= (TF1*)((his->GetListOfFunctions()->FindObject("gaus"))->Clone()); |
| 315 | } |
| 316 | if (fgaus) { |
| 317 | his->Fit(fgaus,"Q","goff"); |
| 318 | meanG=fgaus->GetParameter(1); |
| 319 | rmsG=fgaus->GetParameter(2); |
| 320 | } |
| 321 | } |
| 322 | Double_t dsec=sector-Int_t(sector)-0.5; |
| 323 | Double_t snp=dsec*TMath::Pi()/9.; |
| 324 | (*pcstream)<<"delta"<< |
| 325 | "ptype="<<ptype<< |
| 326 | "dtype="<<dtype<< |
| 327 | "sector="<<sector<< |
| 328 | "dsec="<<dsec<< |
| 329 | "snp="<<snp<< |
| 330 | "phi="<<phi<< |
| 331 | "localX="<<localX<< |
| 332 | "kZ="<<kZ<< |
| 333 | "theta="<<kZ<< |
| 334 | "mean="<<mean<< |
| 335 | "rms="<<rms<< |
| 336 | "meanG="<<meanG<< |
| 337 | "rmsG="<<rmsG<< |
| 338 | "entries="<<entries<< |
| 339 | "meanA="<<meanA<< |
| 340 | "rmsA="<<rmsA<< |
| 341 | "entriesA="<<entriesA<< |
| 342 | "meanC="<<meanC<< |
| 343 | "rmsC="<<rmsC<< |
| 344 | "entriesC="<<entriesC<< |
| 345 | "offsetA="<<offsetA<< |
| 346 | "slopeA="<<slopeA<< |
| 347 | "chi2A="<<chi2A<< |
| 348 | "offsetC="<<offsetC<< |
| 349 | "slopeC="<<slopeC<< |
| 350 | "chi2C="<<chi2C<< |
| 351 | "\n"; |
| 352 | delete his; |
| 353 | } |
| 354 | delete his2; |
| 355 | } |
| 356 | delete his1; |
| 357 | } |
| 358 | hisResMap3D->Write(); |
| 359 | hisResMap2D[0]->Write(); |
| 360 | hisResMap2D[1]->Write(); |
| 361 | hisResMap2D[2]->Write(); |
| 362 | hisResMap2D[3]->Write(); |
| 363 | // delete pcstream; |
| 364 | } |
| 365 | |
| 366 | |
| 367 | |
| 368 | void AliTPCCorrectionFit::MakeDistortionMap(THnSparse * his0, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type, Int_t integ, Double_t bz){ |
| 369 | // |
| 370 | // make a distortion map out ou fthe residual histogram |
| 371 | // Results are written to the debug streamer - pcstream |
| 372 | // Parameters: |
| 373 | // his0 - input (4D) residual histogram |
| 374 | // pcstream - file to write the tree |
| 375 | // run - run number |
| 376 | // refX - track matching reference X |
| 377 | // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt |
| 378 | // THnSparse axes: |
| 379 | // OBJ: TAxis #Delta #Delta |
| 380 | // OBJ: TAxis tanTheta tan(#Theta) |
| 381 | // OBJ: TAxis phi #phi |
| 382 | // OBJ: TAxis snp snp |
| 383 | |
| 384 | // marian.ivanov@cern.ch |
| 385 | const Int_t kMinEntries=10; |
| 386 | Int_t idim[4]={0,1,2,3}; |
| 387 | // |
| 388 | // |
| 389 | // |
| 390 | Int_t nbins3=his0->GetAxis(3)->GetNbins(); |
| 391 | Int_t first3=his0->GetAxis(3)->GetFirst(); |
| 392 | Int_t last3 =his0->GetAxis(3)->GetLast(); |
| 393 | // |
| 394 | for (Int_t ibin3=first3; ibin3<last3; ibin3+=1){ // axis 3 - local angle |
| 395 | his0->GetAxis(3)->SetRange(TMath::Max(ibin3-integ,1),TMath::Min(ibin3+integ,nbins3)); |
| 396 | Double_t x3= his0->GetAxis(3)->GetBinCenter(ibin3); |
| 397 | THnSparse * his3= his0->Projection(3,idim); //projected histogram according selection 3 |
| 398 | // |
| 399 | Int_t nbins2 = his3->GetAxis(2)->GetNbins(); |
| 400 | Int_t first2 = his3->GetAxis(2)->GetFirst(); |
| 401 | Int_t last2 = his3->GetAxis(2)->GetLast(); |
| 402 | // |
| 403 | for (Int_t ibin2=first2; ibin2<last2; ibin2+=1){ // axis 2 - phi |
| 404 | his3->GetAxis(2)->SetRange(TMath::Max(ibin2-integ,1),TMath::Min(ibin2+integ,nbins2)); |
| 405 | Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2); |
| 406 | THnSparse * his2= his3->Projection(2,idim); //projected histogram according selection 2 |
| 407 | Int_t nbins1 = his2->GetAxis(1)->GetNbins(); |
| 408 | Int_t first1 = his2->GetAxis(1)->GetFirst(); |
| 409 | Int_t last1 = his2->GetAxis(1)->GetLast(); |
| 410 | for (Int_t ibin1=first1; ibin1<last1; ibin1++){ //axis 1 - theta |
| 411 | // |
| 412 | Double_t x1= his2->GetAxis(1)->GetBinCenter(ibin1); |
| 413 | his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1)); |
| 414 | if (TMath::Abs(x1)<0.1){ |
| 415 | if (x1<0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1,nbins1)); |
| 416 | if (x1>0) his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1+1,nbins1)); |
| 417 | } |
| 418 | if (TMath::Abs(x1)<0.06){ |
| 419 | his2->GetAxis(1)->SetRange(TMath::Max(ibin1,1),TMath::Min(ibin1,nbins1)); |
| 420 | } |
| 421 | TH1 * hisDelta = his2->Projection(0); |
| 422 | // |
| 423 | Double_t entries = hisDelta->GetEntries(); |
| 424 | Double_t mean=0, rms=0; |
| 425 | if (entries>kMinEntries){ |
| 426 | mean = hisDelta->GetMean(); |
| 427 | rms = hisDelta->GetRMS(); |
| 428 | } |
| 429 | Double_t sector = 9.*x2/TMath::Pi(); |
| 430 | if (sector<0) sector+=18; |
| 431 | Double_t dsec = sector-Int_t(sector)-0.5; |
| 432 | Double_t z=refX*x1; |
| 433 | (*pcstream)<<hname<< |
| 434 | "run="<<run<< |
| 435 | "bz="<<bz<< |
| 436 | "theta="<<x1<< |
| 437 | "phi="<<x2<< |
| 438 | "z="<<z<< // dummy z |
| 439 | "snp="<<x3<< |
| 440 | "entries="<<entries<< |
| 441 | "mean="<<mean<< |
| 442 | "rms="<<rms<< |
| 443 | "refX="<<refX<< // track matching refernce plane |
| 444 | "type="<<type<< // |
| 445 | "sector="<<sector<< |
| 446 | "dsec="<<dsec<< |
| 447 | "\n"; |
| 448 | delete hisDelta; |
| 449 | //printf("%f\t%f\t%f\t%f\t%f\n",x3,x2,x1, entries,mean); |
| 450 | } |
| 451 | delete his2; |
| 452 | } |
| 453 | delete his3; |
| 454 | } |
| 455 | } |
| 456 | |
| 457 | |
| 458 | |
| 459 | void AliTPCCorrectionFit::MakeDistortionMapCosmic(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Float_t refX, Int_t type){ |
| 460 | // |
| 461 | // make a distortion map out ou fthe residual histogram |
| 462 | // Results are written to the debug streamer - pcstream |
| 463 | // Parameters: |
| 464 | // his0 - input (4D) residual histogram |
| 465 | // pcstream - file to write the tree |
| 466 | // run - run number |
| 467 | // refX - track matching reference X |
| 468 | // type - 0- y 1-z,2 -snp, 3-theta, 4=1/pt |
| 469 | // marian.ivanov@cern.ch |
| 470 | // |
| 471 | // Histo axeses |
| 472 | // Collection name='TObjArray', class='TObjArray', size=16 |
| 473 | // 0. OBJ: TAxis #Delta #Delta |
| 474 | // 1. OBJ: TAxis N_{cl} N_{cl} |
| 475 | // 2. OBJ: TAxis dca_{r} (cm) dca_{r} (cm) |
| 476 | // 3. OBJ: TAxis z (cm) z (cm) |
| 477 | // 4. OBJ: TAxis sin(#phi) sin(#phi) |
| 478 | // 5. OBJ: TAxis tan(#theta) tan(#theta) |
| 479 | // 6. OBJ: TAxis 1/pt (1/GeV) 1/pt (1/GeV) |
| 480 | // 7. OBJ: TAxis pt (GeV) pt (GeV) |
| 481 | // 8. OBJ: TAxis alpha alpha |
| 482 | const Int_t kMinEntries=10; |
| 483 | // |
| 484 | // 1. make default selections |
| 485 | // |
| 486 | TH1 * hisDelta=0; |
| 487 | Int_t idim0[4]={0 , 5, 8, 3}; // delta, theta, alpha, z |
| 488 | hisInput->GetAxis(1)->SetRangeUser(110,190); //long tracks |
| 489 | hisInput->GetAxis(2)->SetRangeUser(-10,35); //tracks close to beam pipe |
| 490 | hisInput->GetAxis(4)->SetRangeUser(-0.3,0.3); //small snp at TPC entrance |
| 491 | hisInput->GetAxis(7)->SetRangeUser(3,100); //"high pt tracks" |
| 492 | hisDelta= hisInput->Projection(0); |
| 493 | hisInput->GetAxis(0)->SetRangeUser(-6.*hisDelta->GetRMS(), +6.*hisDelta->GetRMS()); |
| 494 | delete hisDelta; |
| 495 | THnSparse *his0= hisInput->Projection(4,idim0); |
| 496 | // |
| 497 | // 2. Get mean in diferent bins |
| 498 | // |
| 499 | Int_t nbins1=his0->GetAxis(1)->GetNbins(); |
| 500 | Int_t first1=his0->GetAxis(1)->GetFirst(); |
| 501 | Int_t last1 =his0->GetAxis(1)->GetLast(); |
| 502 | // |
| 503 | Double_t bz=AliTrackerBase::GetBz(); |
| 504 | Int_t idim[4]={0,1, 2, 3}; // delta, theta,alpha,z |
| 505 | // |
| 506 | for (Int_t ibin1=first1; ibin1<=last1; ibin1++){ //axis 1 - theta |
| 507 | // |
| 508 | Double_t x1= his0->GetAxis(1)->GetBinCenter(ibin1); |
| 509 | his0->GetAxis(1)->SetRange(TMath::Max(ibin1-1,1),TMath::Min(ibin1+1,nbins1)); |
| 510 | // |
| 511 | THnSparse * his1 = his0->Projection(4,idim); // projected histogram according range1 |
| 512 | Int_t nbins3 = his1->GetAxis(3)->GetNbins(); |
| 513 | Int_t first3 = his1->GetAxis(3)->GetFirst(); |
| 514 | Int_t last3 = his1->GetAxis(3)->GetLast(); |
| 515 | // |
| 516 | for (Int_t ibin3=first3-1; ibin3<=last3; ibin3+=1){ // axis 3 - z at "vertex" |
| 517 | his1->GetAxis(3)->SetRange(TMath::Max(ibin3-1,1),TMath::Min(ibin3+1,nbins3)); |
| 518 | Double_t x3= his1->GetAxis(3)->GetBinCenter(ibin3); |
| 519 | if (ibin3<first3) { |
| 520 | his1->GetAxis(3)->SetRangeUser(-1,1); |
| 521 | x3=0; |
| 522 | } |
| 523 | THnSparse * his3= his1->Projection(4,idim); //projected histogram according selection 3 |
| 524 | Int_t nbins2 = his3->GetAxis(2)->GetNbins(); |
| 525 | Int_t first2 = his3->GetAxis(2)->GetFirst(); |
| 526 | Int_t last2 = his3->GetAxis(2)->GetLast(); |
| 527 | // |
| 528 | for (Int_t ibin2=first2; ibin2<=last2; ibin2+=1){ |
| 529 | his3->GetAxis(2)->SetRange(TMath::Max(ibin2-1,1),TMath::Min(ibin2+1,nbins2)); |
| 530 | Double_t x2= his3->GetAxis(2)->GetBinCenter(ibin2); |
| 531 | hisDelta = his3->Projection(0); |
| 532 | // |
| 533 | Double_t entries = hisDelta->GetEntries(); |
| 534 | Double_t mean=0, rms=0; |
| 535 | if (entries>kMinEntries){ |
| 536 | mean = hisDelta->GetMean(); |
| 537 | rms = hisDelta->GetRMS(); |
| 538 | } |
| 539 | Double_t sector = 9.*x2/TMath::Pi(); |
| 540 | if (sector<0) sector+=18; |
| 541 | Double_t dsec = sector-Int_t(sector)-0.5; |
| 542 | Double_t snp=0; // dummy snp - equal 0 |
| 543 | (*pcstream)<<hname<< |
| 544 | "run="<<run<< |
| 545 | "bz="<<bz<< // magnetic field |
| 546 | "theta="<<x1<< // theta |
| 547 | "phi="<<x2<< // phi (alpha) |
| 548 | "z="<<x3<< // z at "vertex" |
| 549 | "snp="<<snp<< // dummy snp |
| 550 | "entries="<<entries<< // entries in bin |
| 551 | "mean="<<mean<< // mean |
| 552 | "rms="<<rms<< |
| 553 | "refX="<<refX<< // track matching refernce plane |
| 554 | "type="<<type<< // parameter type |
| 555 | "sector="<<sector<< // sector |
| 556 | "dsec="<<dsec<< // dummy delta sector |
| 557 | "\n"; |
| 558 | delete hisDelta; |
| 559 | printf("%f\t%f\t%f\t%f\t%f\n",x1,x3,x2, entries,mean); |
| 560 | } |
| 561 | delete his3; |
| 562 | } |
| 563 | delete his1; |
| 564 | } |
| 565 | delete his0; |
| 566 | } |
| 567 | |
| 568 | |
| 569 | |
| 570 | void AliTPCCorrectionFit::MakeDistortionMapSector(THnSparse * hisInput, TTreeSRedirector * const pcstream, const char* hname, Int_t run, Int_t type, Double_t bz){ |
| 571 | // |
| 572 | // make a distortion map out of the residual histogram |
| 573 | // Results are written to the debug streamer - pcstream |
| 574 | // Parameters: |
| 575 | // his0 - input (4D) residual histogram |
| 576 | // pcstream - file to write the tree |
| 577 | // run - run number |
| 578 | // type - 0- y 1-z,2 -snp, 3-theta |
| 579 | // bz - magnetic field |
| 580 | // marian.ivanov@cern.ch |
| 581 | |
| 582 | //Collection name='TObjArray', class='TObjArray', size=16 |
| 583 | //0 OBJ: TAxis delta delta |
| 584 | //1 OBJ: TAxis phi phi |
| 585 | //2 OBJ: TAxis localX localX |
| 586 | //3 OBJ: TAxis kY kY |
| 587 | //4 OBJ: TAxis kZ kZ |
| 588 | //5 OBJ: TAxis is1 is1 |
| 589 | //6 OBJ: TAxis is0 is0 |
| 590 | //7. OBJ: TAxis z z |
| 591 | //8. OBJ: TAxis IsPrimary IsPrimary |
| 592 | |
| 593 | const Int_t kMinEntries=10; |
| 594 | THnSparse * hisSector0=0; |
| 595 | TH1 * htemp=0; // histogram to calculate mean value of parameter |
| 596 | // Double_t bz=AliTrackerBase::GetBz(); |
| 597 | |
| 598 | // |
| 599 | // Loop over pair of sector: |
| 600 | // isPrim - 8 ==> 8 |
| 601 | // isec0 - 6 ==> 7 |
| 602 | // isec1 - 5 ==> 6 |
| 603 | // refX - 2 ==> 5 |
| 604 | // |
| 605 | // phi - 1 ==> 4 |
| 606 | // z - 7 ==> 3 |
| 607 | // snp - 3 ==> 2 |
| 608 | // theta- 4 ==> 1 |
| 609 | // 0 ==> 0; |
| 610 | for (Int_t isec0=0; isec0<72; isec0++){ |
| 611 | Int_t index0[9]={0, 4, 3, 7, 1, 2, 5, 6,8}; //regroup indeces |
| 612 | // |
| 613 | //hisInput->GetAxis(8)->SetRangeUser(-0.1,0.4); // select secondaries only ? - get out later ? |
| 614 | hisInput->GetAxis(6)->SetRangeUser(isec0-0.1,isec0+0.1); |
| 615 | hisSector0=hisInput->Projection(7,index0); |
| 616 | // |
| 617 | // |
| 618 | for (Int_t isec1=isec0+1; isec1<72; isec1++){ |
| 619 | //if (isec1!=isec0+36) continue; |
| 620 | if ( TMath::Abs((isec0%18)-(isec1%18))>1.5 && TMath::Abs((isec0%18)-(isec1%18))<16.5) continue; |
| 621 | printf("Sectors %d\t%d\n",isec1,isec0); |
| 622 | hisSector0->GetAxis(6)->SetRangeUser(isec1-0.1,isec1+0.1); |
| 623 | TH1 * hisX=hisSector0->Projection(5); |
| 624 | Double_t refX= hisX->GetMean(); |
| 625 | delete hisX; |
| 626 | TH1 *hisDelta=hisSector0->Projection(0); |
| 627 | Double_t dmean = hisDelta->GetMean(); |
| 628 | Double_t drms = hisDelta->GetRMS(); |
| 629 | hisSector0->GetAxis(0)->SetRangeUser(dmean-5.*drms, dmean+5.*drms); |
| 630 | delete hisDelta; |
| 631 | // |
| 632 | // 1. make default selections |
| 633 | // |
| 634 | Int_t idim0[5]={0 , 1, 2, 3, 4}; // {delta, theta, snp, z, phi } |
| 635 | THnBase *hisSector1= hisSector0->ProjectionND(5,idim0); |
| 636 | // |
| 637 | // 2. Get mean in diferent bins |
| 638 | // |
| 639 | Int_t idim[5]={0, 1, 2, 3, 4}; // {delta, theta-1,snp-2 ,z-3, phi-4} |
| 640 | // |
| 641 | // Int_t nbinsPhi=hisSector1->GetAxis(4)->GetNbins(); |
| 642 | Int_t firstPhi=hisSector1->GetAxis(4)->GetFirst(); |
| 643 | Int_t lastPhi =hisSector1->GetAxis(4)->GetLast(); |
| 644 | // |
| 645 | for (Int_t ibinPhi=firstPhi; ibinPhi<=lastPhi; ibinPhi+=2){ //axis 4 - phi |
| 646 | // |
| 647 | // Phi loop |
| 648 | // |
| 649 | Double_t xPhi= hisSector1->GetAxis(4)->GetBinCenter(ibinPhi); |
| 650 | Double_t psec = (9*xPhi/TMath::Pi()); |
| 651 | if (psec<0) psec+=18; |
| 652 | Bool_t isOK0=kFALSE; |
| 653 | Bool_t isOK1=kFALSE; |
| 654 | if (TMath::Abs(psec-isec0%18-0.5)<1. || TMath::Abs(psec-isec0%18-17.5)<1.) isOK0=kTRUE; |
| 655 | if (TMath::Abs(psec-isec1%18-0.5)<1. || TMath::Abs(psec-isec1%18-17.5)<1.) isOK1=kTRUE; |
| 656 | if (!isOK0) continue; |
| 657 | if (!isOK1) continue; |
| 658 | // |
| 659 | hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-2,firstPhi),TMath::Min(ibinPhi+2,lastPhi)); |
| 660 | if (isec1!=isec0+36) { |
| 661 | hisSector1->GetAxis(4)->SetRange(TMath::Max(ibinPhi-3,firstPhi),TMath::Min(ibinPhi+3,lastPhi)); |
| 662 | } |
| 663 | // |
| 664 | htemp = hisSector1->Projection(4); |
| 665 | xPhi=htemp->GetMean(); |
| 666 | delete htemp; |
| 667 | THnBase * hisPhi = hisSector1->ProjectionND(4,idim); |
| 668 | //Int_t nbinsZ = hisPhi->GetAxis(3)->GetNbins(); |
| 669 | Int_t firstZ = hisPhi->GetAxis(3)->GetFirst(); |
| 670 | Int_t lastZ = hisPhi->GetAxis(3)->GetLast(); |
| 671 | // |
| 672 | for (Int_t ibinZ=firstZ; ibinZ<=lastZ; ibinZ+=2){ // axis 3 - z |
| 673 | // |
| 674 | // Z loop |
| 675 | // |
| 676 | hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ,firstZ),TMath::Min(ibinZ,lastZ)); |
| 677 | if (isec1!=isec0+36) { |
| 678 | hisPhi->GetAxis(3)->SetRange(TMath::Max(ibinZ-1,firstZ),TMath::Min(ibinZ-1,lastZ)); |
| 679 | } |
| 680 | htemp = hisPhi->Projection(3); |
| 681 | Double_t xZ= htemp->GetMean(); |
| 682 | delete htemp; |
| 683 | THnBase * hisZ= hisPhi->ProjectionND(3,idim); |
| 684 | //projected histogram according selection 3 -z |
| 685 | // |
| 686 | // |
| 687 | //Int_t nbinsSnp = hisZ->GetAxis(2)->GetNbins(); |
| 688 | Int_t firstSnp = hisZ->GetAxis(2)->GetFirst(); |
| 689 | Int_t lastSnp = hisZ->GetAxis(2)->GetLast(); |
| 690 | for (Int_t ibinSnp=firstSnp; ibinSnp<=lastSnp; ibinSnp+=2){ // axis 2 - snp |
| 691 | // |
| 692 | // Snp loop |
| 693 | // |
| 694 | hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-1,firstSnp),TMath::Min(ibinSnp+1,lastSnp)); |
| 695 | if (isec1!=isec0+36) { |
| 696 | hisZ->GetAxis(2)->SetRange(TMath::Max(ibinSnp-2,firstSnp),TMath::Min(ibinSnp+2,lastSnp)); |
| 697 | } |
| 698 | htemp = hisZ->Projection(2); |
| 699 | Double_t xSnp= htemp->GetMean(); |
| 700 | delete htemp; |
| 701 | THnBase * hisSnp= hisZ->ProjectionND(2,idim); |
| 702 | //projected histogram according selection 2 - snp |
| 703 | |
| 704 | //Int_t nbinsTheta = hisSnp->GetAxis(1)->GetNbins(); |
| 705 | Int_t firstTheta = hisSnp->GetAxis(1)->GetFirst(); |
| 706 | Int_t lastTheta = hisSnp->GetAxis(1)->GetLast(); |
| 707 | // |
| 708 | for (Int_t ibinTheta=firstTheta; ibinTheta<=lastTheta; ibinTheta+=2){ // axis1 theta |
| 709 | |
| 710 | |
| 711 | hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-2,firstTheta),TMath::Min(ibinTheta+2,lastTheta)); |
| 712 | if (isec1!=isec0+36) { |
| 713 | hisSnp->GetAxis(1)->SetRange(TMath::Max(ibinTheta-3,firstTheta),TMath::Min(ibinTheta+3,lastTheta)); |
| 714 | } |
| 715 | htemp = hisSnp->Projection(1); |
| 716 | Double_t xTheta=htemp->GetMean(); |
| 717 | delete htemp; |
| 718 | hisDelta = hisSnp->Projection(0); |
| 719 | // |
| 720 | Double_t entries = hisDelta->GetEntries(); |
| 721 | Double_t mean=0, rms=0; |
| 722 | if (entries>kMinEntries){ |
| 723 | mean = hisDelta->GetMean(); |
| 724 | rms = hisDelta->GetRMS(); |
| 725 | } |
| 726 | Double_t sector = 9.*xPhi/TMath::Pi(); |
| 727 | if (sector<0) sector+=18; |
| 728 | Double_t dsec = sector-Int_t(sector)-0.5; |
| 729 | Int_t dtype=1; // TPC alignment type |
| 730 | (*pcstream)<<hname<< |
| 731 | "run="<<run<< |
| 732 | "bz="<<bz<< // magnetic field |
| 733 | "ptype="<<type<< // parameter type |
| 734 | "dtype="<<dtype<< // parameter type |
| 735 | "isec0="<<isec0<< // sector 0 |
| 736 | "isec1="<<isec1<< // sector 1 |
| 737 | "sector="<<sector<< // sector as float |
| 738 | "dsec="<<dsec<< // delta sector |
| 739 | // |
| 740 | "theta="<<xTheta<< // theta |
| 741 | "phi="<<xPhi<< // phi (alpha) |
| 742 | "z="<<xZ<< // z |
| 743 | "snp="<<xSnp<< // snp |
| 744 | |
| 745 | // |
| 746 | "entries="<<entries<< // entries in bin |
| 747 | "mean="<<mean<< // mean |
| 748 | "rms="<<rms<< // rms |
| 749 | "refX="<<refX<< // track matching reference plane |
| 750 | "\n"; |
| 751 | delete hisDelta; |
| 752 | //printf("%d\t%d\t%f\t%f\t%f\t%f\t%f\t%f\n",isec0, isec1, xPhi,xZ,xSnp, xTheta, entries,mean); |
| 753 | // |
| 754 | }//ibinTheta |
| 755 | delete hisSnp; |
| 756 | } //ibinSnp |
| 757 | delete hisZ; |
| 758 | }//ibinZ |
| 759 | delete hisPhi; |
| 760 | }//ibinPhi |
| 761 | delete hisSector1; |
| 762 | }//isec1 |
| 763 | delete hisSector0; |
| 764 | }//isec0 |
| 765 | } |
| 766 | |
| 767 | |
| 768 | |
| 769 | |
| 770 | |
| 771 | // void AliTPCCorrectionFit::MakeLaserDistortionTree(TTree* tree, TObjArray */*corrArray*/, Int_t /*itype*/){ |
| 772 | // // |
| 773 | // // Make a laser fit tree for global minimization |
| 774 | // // |
| 775 | // AliTPCcalibDB* calib=AliTPCcalibDB::Instance(); |
| 776 | // AliTPCCorrection * correction = calib->GetTPCComposedCorrection(); |
| 777 | // if (!correction) correction = calib->GetTPCComposedCorrection(AliTrackerBase::GetBz()); |
| 778 | // correction->AddVisualCorrection(correction,0); //register correction |
| 779 | |
| 780 | // // AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; |
| 781 | // //AliTPCParam *param = AliTPCcalibDB::Instance()->GetParameters(); |
| 782 | // // |
| 783 | // const Double_t cutErrY=0.05; |
| 784 | // const Double_t kSigmaCut=4; |
| 785 | // // const Double_t cutErrZ=0.03; |
| 786 | // const Double_t kEpsilon=0.00000001; |
| 787 | // // const Double_t kMaxDist=1.; // max distance - space correction |
| 788 | // TVectorD *vecdY=0; |
| 789 | // TVectorD *vecdZ=0; |
| 790 | // TVectorD *veceY=0; |
| 791 | // TVectorD *veceZ=0; |
| 792 | // AliTPCLaserTrack *ltr=0; |
| 793 | // AliTPCLaserTrack::LoadTracks(); |
| 794 | // tree->SetBranchAddress("dY.",&vecdY); |
| 795 | // tree->SetBranchAddress("dZ.",&vecdZ); |
| 796 | // tree->SetBranchAddress("eY.",&veceY); |
| 797 | // tree->SetBranchAddress("eZ.",&veceZ); |
| 798 | // tree->SetBranchAddress("LTr.",<r); |
| 799 | // Int_t entries= tree->GetEntries(); |
| 800 | // TTreeSRedirector *pcstream= new TTreeSRedirector("distortionLaser_0.root"); |
| 801 | // Double_t bz=AliTrackerBase::GetBz(); |
| 802 | // // |
| 803 | // // Double_t globalXYZ[3]; |
| 804 | // //Double_t globalXYZCorr[3]; |
| 805 | // for (Int_t ientry=0; ientry<entries; ientry++){ |
| 806 | // tree->GetEntry(ientry); |
| 807 | // if (!ltr->GetVecGX()){ |
| 808 | // ltr->UpdatePoints(); |
| 809 | // } |
| 810 | // // |
| 811 | // TVectorD fit10(5); |
| 812 | // TVectorD fit5(5); |
| 813 | // printf("Entry\t%d\n",ientry); |
| 814 | // for (Int_t irow0=0; irow0<158; irow0+=1){ |
| 815 | // // |
| 816 | // TLinearFitter fitter10(4,"hyp3"); |
| 817 | // TLinearFitter fitter5(2,"hyp1"); |
| 818 | // Int_t sector= (Int_t)(*ltr->GetVecSec())[irow0]; |
| 819 | // if (sector<0) continue; |
| 820 | // //if (TMath::Abs(vecdY->GetMatrixArray()[irow0])<kEpsilon) continue; |
| 821 | |
| 822 | // Double_t refX= (*ltr->GetVecLX())[irow0]; |
| 823 | // Int_t firstRow1 = TMath::Max(irow0-10,0); |
| 824 | // Int_t lastRow1 = TMath::Min(irow0+10,158); |
| 825 | // Double_t padWidth=(irow0<64)?0.4:0.6; |
| 826 | // // make long range fit |
| 827 | // for (Int_t irow1=firstRow1; irow1<=lastRow1; irow1++){ |
| 828 | // if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue; |
| 829 | // if (veceY->GetMatrixArray()[irow1]>cutErrY) continue; |
| 830 | // if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue; |
| 831 | // Double_t idealX= (*ltr->GetVecLX())[irow1]; |
| 832 | // Double_t idealY= (*ltr->GetVecLY())[irow1]; |
| 833 | // // Double_t idealZ= (*ltr->GetVecLZ())[irow1]; |
| 834 | // Double_t gx= (*ltr->GetVecGX())[irow1]; |
| 835 | // Double_t gy= (*ltr->GetVecGY())[irow1]; |
| 836 | // Double_t gz= (*ltr->GetVecGZ())[irow1]; |
| 837 | // Double_t measY=(*vecdY)[irow1]+idealY; |
| 838 | // Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0); |
| 839 | // // deltaR = R distorted -R ideal |
| 840 | // Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)}; |
| 841 | // fitter10.AddPoint(xxx,measY,1); |
| 842 | // } |
| 843 | // Bool_t isOK=kTRUE; |
| 844 | // Double_t rms10=0;//TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4)); |
| 845 | // Double_t mean10 =0;// fitter10.GetParameter(0); |
| 846 | // Double_t slope10 =0;// fitter10.GetParameter(0); |
| 847 | // Double_t cosPart10 = 0;// fitter10.GetParameter(2); |
| 848 | // Double_t sinPart10 =0;// fitter10.GetParameter(3); |
| 849 | |
| 850 | // if (fitter10.GetNpoints()>10){ |
| 851 | // fitter10.Eval(); |
| 852 | // rms10=TMath::Sqrt(fitter10.GetChisquare()/(fitter10.GetNpoints()-4)); |
| 853 | // mean10 = fitter10.GetParameter(0); |
| 854 | // slope10 = fitter10.GetParameter(1); |
| 855 | // cosPart10 = fitter10.GetParameter(2); |
| 856 | // sinPart10 = fitter10.GetParameter(3); |
| 857 | // // |
| 858 | // // make short range fit |
| 859 | // // |
| 860 | // for (Int_t irow1=firstRow1+5; irow1<=lastRow1-5; irow1++){ |
| 861 | // if (TMath::Abs((*ltr->GetVecSec())[irow1]-sector)>kEpsilon) continue; |
| 862 | // if (veceY->GetMatrixArray()[irow1]>cutErrY) continue; |
| 863 | // if (TMath::Abs(vecdY->GetMatrixArray()[irow1])<kEpsilon) continue; |
| 864 | // Double_t idealX= (*ltr->GetVecLX())[irow1]; |
| 865 | // Double_t idealY= (*ltr->GetVecLY())[irow1]; |
| 866 | // // Double_t idealZ= (*ltr->GetVecLZ())[irow1]; |
| 867 | // Double_t gx= (*ltr->GetVecGX())[irow1]; |
| 868 | // Double_t gy= (*ltr->GetVecGY())[irow1]; |
| 869 | // Double_t gz= (*ltr->GetVecGZ())[irow1]; |
| 870 | // Double_t measY=(*vecdY)[irow1]+idealY; |
| 871 | // Double_t deltaR = GetCorrXYZ(gx, gy, gz, 0,0); |
| 872 | // // deltaR = R distorted -R ideal |
| 873 | // Double_t expY= mean10+slope10*(idealX+deltaR-refX); |
| 874 | // if (TMath::Abs(measY-expY)>kSigmaCut*rms10) continue; |
| 875 | // // |
| 876 | // Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth); |
| 877 | // Double_t xxx[4]={idealX+deltaR-refX,TMath::Cos(idealY/padWidth), TMath::Sin(idealY/padWidth)}; |
| 878 | // fitter5.AddPoint(xxx,measY-corr,1); |
| 879 | // } |
| 880 | // }else{ |
| 881 | // isOK=kFALSE; |
| 882 | // } |
| 883 | // if (fitter5.GetNpoints()<8) isOK=kFALSE; |
| 884 | |
| 885 | // Double_t rms5=0;//TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4)); |
| 886 | // Double_t offset5 =0;// fitter5.GetParameter(0); |
| 887 | // Double_t slope5 =0;// fitter5.GetParameter(0); |
| 888 | // if (isOK){ |
| 889 | // fitter5.Eval(); |
| 890 | // rms5=TMath::Sqrt(fitter5.GetChisquare()/(fitter5.GetNpoints()-4)); |
| 891 | // offset5 = fitter5.GetParameter(0); |
| 892 | // slope5 = fitter5.GetParameter(0); |
| 893 | // } |
| 894 | // // |
| 895 | // Double_t dtype=5; |
| 896 | // Double_t ptype=0; |
| 897 | // Double_t phi =(*ltr->GetVecPhi())[irow0]; |
| 898 | // Double_t theta =ltr->GetTgl(); |
| 899 | // Double_t mean=(vecdY)->GetMatrixArray()[irow0]; |
| 900 | // Double_t gx=0,gy=0,gz=0; |
| 901 | // Double_t snp = (*ltr->GetVecP2())[irow0]; |
| 902 | // Int_t bundle= ltr->GetBundle(); |
| 903 | // Int_t id= ltr->GetId(); |
| 904 | // // Double_t rms = err->GetMatrixArray()[irow]; |
| 905 | // // |
| 906 | // gx = (*ltr->GetVecGX())[irow0]; |
| 907 | // gy = (*ltr->GetVecGY())[irow0]; |
| 908 | // gz = (*ltr->GetVecGZ())[irow0]; |
| 909 | // Double_t dRrec = GetCorrXYZ(gx, gy, gz, 0,0); |
| 910 | // fitter10.GetParameters(fit10); |
| 911 | // fitter5.GetParameters(fit5); |
| 912 | // Double_t idealY= (*ltr->GetVecLY())[irow0]; |
| 913 | // Double_t measY=(*vecdY)[irow0]+idealY; |
| 914 | // Double_t corr=cosPart10*TMath::Cos(idealY/padWidth)+sinPart10*TMath::Sin(idealY/padWidth); |
| 915 | // if (TMath::Max(rms5,rms10)>0.06) isOK=kFALSE; |
| 916 | // // |
| 917 | // (*pcstream)<<"fitFull"<< // dumpe also intermediate results |
| 918 | // "bz="<<bz<< // magnetic filed used |
| 919 | // "dtype="<<dtype<< // detector match type |
| 920 | // "ptype="<<ptype<< // parameter type |
| 921 | // "theta="<<theta<< // theta |
| 922 | // "phi="<<phi<< // phi |
| 923 | // "snp="<<snp<< // snp |
| 924 | // "sector="<<sector<< |
| 925 | // "bundle="<<bundle<< |
| 926 | // // // "dsec="<<dsec<< |
| 927 | // "refX="<<refX<< // reference radius |
| 928 | // "gx="<<gx<< // global position |
| 929 | // "gy="<<gy<< // global position |
| 930 | // "gz="<<gz<< // global position |
| 931 | // "dRrec="<<dRrec<< // delta Radius in reconstruction |
| 932 | // "id="<<id<< //bundle |
| 933 | // "rms10="<<rms10<< |
| 934 | // "rms5="<<rms5<< |
| 935 | // "fit10.="<<&fit10<< |
| 936 | // "fit5.="<<&fit5<< |
| 937 | // "measY="<<measY<< |
| 938 | // "mean="<<mean<< |
| 939 | // "idealY="<<idealY<< |
| 940 | // "corr="<<corr<< |
| 941 | // "isOK="<<isOK<< |
| 942 | // "\n"; |
| 943 | // } |
| 944 | // } |
| 945 | // delete pcstream; |
| 946 | // } |
| 947 | |
| 948 | |
| 949 | void AliTPCCorrectionFit::MakeTrackDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){ |
| 950 | // |
| 951 | // Make a fit tree: |
| 952 | // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX) |
| 953 | // calculates partial distortions |
| 954 | // Partial distortion is stored in the resulting tree |
| 955 | // Output is storred in the file distortion_<dettype>_<partype>.root |
| 956 | // Partial distortion is stored with the name given by correction name |
| 957 | // |
| 958 | // |
| 959 | // Parameters of function: |
| 960 | // input - input tree |
| 961 | // dtype - distortion type 0 - ITSTPC, 1 -TPCTRD, 2 - TPCvertex , 3 - TPC-TOF, 4 - TPCTPC track crossing |
| 962 | // ppype - parameter type |
| 963 | // corrArray - array with partial corrections |
| 964 | // step - skipe entries - if 1 all entries processed - it is slow |
| 965 | // debug 0 if debug on also space points dumped - it is slow |
| 966 | |
| 967 | const Double_t kMaxSnp = 0.85; |
| 968 | const Double_t kcutSnp=0.25; |
| 969 | const Double_t kcutTheta=1.; |
| 970 | const Double_t kRadiusTPC=85; |
| 971 | // AliTPCROC *tpcRoc =AliTPCROC::Instance(); |
| 972 | // |
| 973 | const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); |
| 974 | // const Double_t kB2C=-0.299792458e-3; |
| 975 | const Int_t kMinEntries=20; |
| 976 | Double_t phi,theta, snp, mean,rms, entries,sector,dsec; |
| 977 | Float_t refX; |
| 978 | Int_t run; |
| 979 | tinput->SetBranchAddress("run",&run); |
| 980 | tinput->SetBranchAddress("theta",&theta); |
| 981 | tinput->SetBranchAddress("phi", &phi); |
| 982 | tinput->SetBranchAddress("snp",&snp); |
| 983 | tinput->SetBranchAddress("mean",&mean); |
| 984 | tinput->SetBranchAddress("rms",&rms); |
| 985 | tinput->SetBranchAddress("entries",&entries); |
| 986 | tinput->SetBranchAddress("sector",§or); |
| 987 | tinput->SetBranchAddress("dsec",&dsec); |
| 988 | tinput->SetBranchAddress("refX",&refX); |
| 989 | TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortion%d_%d_%d.root",dtype,ptype,offset)); |
| 990 | // |
| 991 | Int_t nentries=tinput->GetEntries(); |
| 992 | Int_t ncorr=corrArray->GetEntries(); |
| 993 | Double_t corrections[100]={0}; // |
| 994 | Double_t tPar[5]; |
| 995 | Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; |
| 996 | Int_t dir=0; |
| 997 | if (dtype==5 || dtype==6) dtype=4; |
| 998 | if (dtype==0) { dir=-1;} |
| 999 | if (dtype==1) { dir=1;} |
| 1000 | if (dtype==2) { dir=-1;} |
| 1001 | if (dtype==3) { dir=1;} |
| 1002 | if (dtype==4) { dir=-1;} |
| 1003 | // |
| 1004 | for (Int_t ientry=offset; ientry<nentries; ientry+=step){ |
| 1005 | tinput->GetEntry(ientry); |
| 1006 | if (TMath::Abs(snp)>kMaxSnp) continue; |
| 1007 | tPar[0]=0; |
| 1008 | tPar[1]=theta*refX; |
| 1009 | if (dtype==2) tPar[1]=theta*kRadiusTPC; |
| 1010 | tPar[2]=snp; |
| 1011 | tPar[3]=theta; |
| 1012 | tPar[4]=(gRandom->Rndm()-0.5)*0.02; // should be calculated - non equal to 0 |
| 1013 | if (dtype==4){ |
| 1014 | // tracks crossing CE |
| 1015 | tPar[1]=0; // track at the CE |
| 1016 | //if (TMath::Abs(theta) <0.05) continue; // deep cross |
| 1017 | } |
| 1018 | |
| 1019 | if (TMath::Abs(snp) >kcutSnp) continue; |
| 1020 | if (TMath::Abs(theta) >kcutTheta) continue; |
| 1021 | printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms); |
| 1022 | Double_t bz=AliTrackerBase::GetBz(); |
| 1023 | if (dtype !=4) { //exclude TPC - for TPC mainly non primary tracks |
| 1024 | if (dtype!=2 && TMath::Abs(bz)>0.1 ) tPar[4]=snp/(refX*bz*kB2C*2); |
| 1025 | |
| 1026 | if (dtype==2 && TMath::Abs(bz)>0.1 ) { |
| 1027 | tPar[4]=snp/(kRadiusTPC*bz*kB2C*2);// |
| 1028 | // snp at the TPC inner radius in case the vertex match used |
| 1029 | } |
| 1030 | } |
| 1031 | // |
| 1032 | tPar[4]+=(gRandom->Rndm()-0.5)*0.02; |
| 1033 | AliExternalTrackParam track(refX,phi,tPar,cov); |
| 1034 | Double_t xyz[3]; |
| 1035 | track.GetXYZ(xyz); |
| 1036 | Int_t id=0; |
| 1037 | Double_t pt=1./tPar[4]; |
| 1038 | Double_t dRrec=0; // dummy value - needed for points - e.g for laser |
| 1039 | //if (ptype==4 &&bz<0) mean*=-1; // interpret as curvature -- COMMENTED out - in lookup signed 1/pt used |
| 1040 | Double_t refXD=refX; |
| 1041 | (*pcstream)<<"fit"<< |
| 1042 | "run="<<run<< // run number |
| 1043 | "bz="<<bz<< // magnetic filed used |
| 1044 | "dtype="<<dtype<< // detector match type |
| 1045 | "ptype="<<ptype<< // parameter type |
| 1046 | "theta="<<theta<< // theta |
| 1047 | "phi="<<phi<< // phi |
| 1048 | "snp="<<snp<< // snp |
| 1049 | "mean="<<mean<< // mean dist value |
| 1050 | "rms="<<rms<< // rms |
| 1051 | "sector="<<sector<< |
| 1052 | "dsec="<<dsec<< |
| 1053 | "refX="<<refXD<< // referece X as double |
| 1054 | "gx="<<xyz[0]<< // global position at reference |
| 1055 | "gy="<<xyz[1]<< // global position at reference |
| 1056 | "gz="<<xyz[2]<< // global position at reference |
| 1057 | "dRrec="<<dRrec<< // delta Radius in reconstruction |
| 1058 | "pt="<<pt<< // pt |
| 1059 | "id="<<id<< // track id |
| 1060 | "entries="<<entries;// number of entries in bin |
| 1061 | // |
| 1062 | Bool_t isOK=kTRUE; |
| 1063 | if (entries<kMinEntries) isOK=kFALSE; |
| 1064 | // |
| 1065 | if (dtype!=4) for (Int_t icorr=0; icorr<ncorr; icorr++) { |
| 1066 | AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr); |
| 1067 | corrections[icorr]=0; |
| 1068 | if (entries>kMinEntries){ |
| 1069 | AliExternalTrackParam trackIn(refX,phi,tPar,cov); |
| 1070 | AliExternalTrackParam *trackOut = 0; |
| 1071 | if (debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,pcstream); |
| 1072 | if (!debug) trackOut=corr->FitDistortedTrack(trackIn, refX, dir,0); |
| 1073 | if (dtype==0) {dir= -1;} |
| 1074 | if (dtype==1) {dir= 1;} |
| 1075 | if (dtype==2) {dir= -1;} |
| 1076 | if (dtype==3) {dir= 1;} |
| 1077 | // |
| 1078 | if (trackOut){ |
| 1079 | if (!AliTrackerBase::PropagateTrackTo(&trackIn,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE; |
| 1080 | if (!trackOut->Rotate(trackIn.GetAlpha())) isOK=kFALSE; |
| 1081 | if (!AliTrackerBase::PropagateTrackTo(trackOut,trackIn.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE; |
| 1082 | // trackOut->PropagateTo(trackIn.GetX(),AliTrackerBase::GetBz()); |
| 1083 | // |
| 1084 | corrections[icorr]= trackOut->GetParameter()[ptype]-trackIn.GetParameter()[ptype]; |
| 1085 | delete trackOut; |
| 1086 | }else{ |
| 1087 | corrections[icorr]=0; |
| 1088 | isOK=kFALSE; |
| 1089 | } |
| 1090 | //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out |
| 1091 | } |
| 1092 | (*pcstream)<<"fit"<< |
| 1093 | Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value |
| 1094 | } |
| 1095 | |
| 1096 | if (dtype==4) for (Int_t icorr=0; icorr<ncorr; icorr++) { |
| 1097 | // |
| 1098 | // special case of the TPC tracks crossing the CE |
| 1099 | // |
| 1100 | AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr); |
| 1101 | corrections[icorr]=0; |
| 1102 | if (entries>kMinEntries){ |
| 1103 | AliExternalTrackParam trackIn0(refX,phi,tPar,cov); //Outer - direction to vertex |
| 1104 | AliExternalTrackParam trackIn1(refX,phi,tPar,cov); //Inner - direction magnet |
| 1105 | AliExternalTrackParam *trackOut0 = 0; |
| 1106 | AliExternalTrackParam *trackOut1 = 0; |
| 1107 | // |
| 1108 | if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream); |
| 1109 | if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0); |
| 1110 | if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream); |
| 1111 | if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0); |
| 1112 | // |
| 1113 | if (trackOut0 && trackOut1){ |
| 1114 | if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE; |
| 1115 | if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; |
| 1116 | if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE; |
| 1117 | if (!AliTrackerBase::PropagateTrackTo(trackOut0,trackIn0.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE; |
| 1118 | // |
| 1119 | if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,5,kTRUE,kMaxSnp)) isOK=kFALSE; |
| 1120 | if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE; |
| 1121 | if (!AliTrackerBase::PropagateTrackTo(&trackIn1,trackIn0.GetX(),kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; |
| 1122 | if (!trackOut1->Rotate(trackIn1.GetAlpha())) isOK=kFALSE; |
| 1123 | if (!AliTrackerBase::PropagateTrackTo(trackOut1,trackIn1.GetX(),kMass,5,kFALSE,kMaxSnp)) isOK=kFALSE; |
| 1124 | // |
| 1125 | corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]); |
| 1126 | corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]); |
| 1127 | if (isOK) |
| 1128 | if ((TMath::Abs(trackOut0->GetX()-trackOut1->GetX())>0.1)|| |
| 1129 | (TMath::Abs(trackOut0->GetX()-trackIn1.GetX())>0.1)|| |
| 1130 | (TMath::Abs(trackOut0->GetAlpha()-trackOut1->GetAlpha())>0.00001)|| |
| 1131 | (TMath::Abs(trackOut0->GetAlpha()-trackIn1.GetAlpha())>0.00001)|| |
| 1132 | (TMath::Abs(trackIn0.GetTgl()-trackIn1.GetTgl())>0.0001)|| |
| 1133 | (TMath::Abs(trackIn0.GetSnp()-trackIn1.GetSnp())>0.0001) |
| 1134 | ){ |
| 1135 | isOK=kFALSE; |
| 1136 | } |
| 1137 | delete trackOut0; |
| 1138 | delete trackOut1; |
| 1139 | }else{ |
| 1140 | corrections[icorr]=0; |
| 1141 | isOK=kFALSE; |
| 1142 | } |
| 1143 | // |
| 1144 | //if (ptype==4 &&bz<0) corrections[icorr]*=-1; // interpret as curvature - commented out no in lookup |
| 1145 | } |
| 1146 | (*pcstream)<<"fit"<< |
| 1147 | Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value |
| 1148 | } |
| 1149 | // |
| 1150 | (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n"; |
| 1151 | } |
| 1152 | |
| 1153 | |
| 1154 | delete pcstream; |
| 1155 | } |
| 1156 | |
| 1157 | |
| 1158 | |
| 1159 | void AliTPCCorrectionFit::MakeSectorDistortionTree(TTree *tinput, Int_t dtype, Int_t ptype, const TObjArray * corrArray, Int_t step, Int_t offset, Bool_t debug ){ |
| 1160 | // |
| 1161 | // Make a fit tree: |
| 1162 | // For each partial correction (specified in array) and given track topology (phi, theta, snp, refX) |
| 1163 | // calculates partial distortions |
| 1164 | // Partial distortion is stored in the resulting tree |
| 1165 | // Output is storred in the file distortion_<dettype>_<partype>.root |
| 1166 | // Partial distortion is stored with the name given by correction name |
| 1167 | // |
| 1168 | // |
| 1169 | // Parameters of function: |
| 1170 | // input - input tree |
| 1171 | // dtype - distortion type 10 - IROC-OROC |
| 1172 | // ppype - parameter type |
| 1173 | // corrArray - array with partial corrections |
| 1174 | // step - skipe entries - if 1 all entries processed - it is slow |
| 1175 | // debug 0 if debug on also space points dumped - it is slow |
| 1176 | |
| 1177 | const Double_t kMaxSnp = 0.8; |
| 1178 | const Int_t kMinEntries=200; |
| 1179 | // AliTPCROC *tpcRoc =AliTPCROC::Instance(); |
| 1180 | // |
| 1181 | const Double_t kMass = TDatabasePDG::Instance()->GetParticle("pi+")->Mass(); |
| 1182 | // const Double_t kB2C=-0.299792458e-3; |
| 1183 | Double_t phi,theta, snp, mean,rms, entries,sector,dsec,globalZ; |
| 1184 | Int_t isec1, isec0; |
| 1185 | Double_t refXD; |
| 1186 | Float_t refX; |
| 1187 | Int_t run; |
| 1188 | tinput->SetBranchAddress("run",&run); |
| 1189 | tinput->SetBranchAddress("theta",&theta); |
| 1190 | tinput->SetBranchAddress("phi", &phi); |
| 1191 | tinput->SetBranchAddress("snp",&snp); |
| 1192 | tinput->SetBranchAddress("mean",&mean); |
| 1193 | tinput->SetBranchAddress("rms",&rms); |
| 1194 | tinput->SetBranchAddress("entries",&entries); |
| 1195 | tinput->SetBranchAddress("sector",§or); |
| 1196 | tinput->SetBranchAddress("dsec",&dsec); |
| 1197 | tinput->SetBranchAddress("refX",&refXD); |
| 1198 | tinput->SetBranchAddress("z",&globalZ); |
| 1199 | tinput->SetBranchAddress("isec0",&isec0); |
| 1200 | tinput->SetBranchAddress("isec1",&isec1); |
| 1201 | TTreeSRedirector *pcstream = new TTreeSRedirector(Form("distortionSector%d_%d_%d.root",dtype,ptype,offset)); |
| 1202 | // |
| 1203 | Int_t nentries=tinput->GetEntries(); |
| 1204 | Int_t ncorr=corrArray->GetEntries(); |
| 1205 | Double_t corrections[100]={0}; // |
| 1206 | Double_t tPar[5]; |
| 1207 | Double_t cov[15]={0,0,0,0,0,0,0,0,0,0,0,0,0,0}; |
| 1208 | Int_t dir=0; |
| 1209 | // |
| 1210 | for (Int_t ientry=offset; ientry<nentries; ientry+=step){ |
| 1211 | tinput->GetEntry(ientry); |
| 1212 | refX=refXD; |
| 1213 | Int_t id=-1; |
| 1214 | if (TMath::Abs(TMath::Abs(isec0%18)-TMath::Abs(isec1%18))==0) id=1; // IROC-OROC - opposite side |
| 1215 | if (TMath::Abs(TMath::Abs(isec0%36)-TMath::Abs(isec1%36))==0) id=2; // IROC-OROC - same side |
| 1216 | if (dtype==10 && id==-1) continue; |
| 1217 | // |
| 1218 | dir=-1; |
| 1219 | tPar[0]=0; |
| 1220 | tPar[1]=globalZ; |
| 1221 | tPar[2]=snp; |
| 1222 | tPar[3]=theta; |
| 1223 | tPar[4]=(gRandom->Rndm()-0.1)*0.2; // |
| 1224 | Double_t pt=1./tPar[4]; |
| 1225 | // |
| 1226 | printf("%f\t%f\t%f\t%f\t%f\t%f\n",entries, sector,theta,snp, mean,rms); |
| 1227 | Double_t bz=AliTrackerBase::GetBz(); |
| 1228 | AliExternalTrackParam track(refX,phi,tPar,cov); |
| 1229 | Double_t xyz[3],xyzIn[3],xyzOut[3]; |
| 1230 | track.GetXYZ(xyz); |
| 1231 | track.GetXYZAt(85,bz,xyzIn); |
| 1232 | track.GetXYZAt(245,bz,xyzOut); |
| 1233 | Double_t phiIn = TMath::ATan2(xyzIn[1],xyzIn[0]); |
| 1234 | Double_t phiOut = TMath::ATan2(xyzOut[1],xyzOut[0]); |
| 1235 | Double_t phiRef = TMath::ATan2(xyz[1],xyz[0]); |
| 1236 | Int_t sectorRef = TMath::Nint(9.*phiRef/TMath::Pi()-0.5); |
| 1237 | Int_t sectorIn = TMath::Nint(9.*phiIn/TMath::Pi()-0.5); |
| 1238 | Int_t sectorOut = TMath::Nint(9.*phiOut/TMath::Pi()-0.5); |
| 1239 | // |
| 1240 | Bool_t isOK=kTRUE; |
| 1241 | if (sectorIn!=sectorOut) isOK=kFALSE; // requironment - cluster in the same sector |
| 1242 | if (sectorIn!=sectorRef) isOK=kFALSE; // requironment - cluster in the same sector |
| 1243 | if (entries<kMinEntries/(1+TMath::Abs(globalZ/100.))) isOK=kFALSE; // requironment - minimal amount of tracks in bin |
| 1244 | // Do downscale |
| 1245 | if (TMath::Abs(theta)>1) isOK=kFALSE; |
| 1246 | // |
| 1247 | Double_t dRrec=0; // dummy value - needed for points - e.g for laser |
| 1248 | // |
| 1249 | (*pcstream)<<"fit"<< |
| 1250 | "run="<<run<< //run |
| 1251 | "bz="<<bz<< // magnetic filed used |
| 1252 | "dtype="<<dtype<< // detector match type |
| 1253 | "ptype="<<ptype<< // parameter type |
| 1254 | "theta="<<theta<< // theta |
| 1255 | "phi="<<phi<< // phi |
| 1256 | "snp="<<snp<< // snp |
| 1257 | "mean="<<mean<< // mean dist value |
| 1258 | "rms="<<rms<< // rms |
| 1259 | "sector="<<sector<< |
| 1260 | "dsec="<<dsec<< |
| 1261 | "refX="<<refXD<< // referece X |
| 1262 | "gx="<<xyz[0]<< // global position at reference |
| 1263 | "gy="<<xyz[1]<< // global position at reference |
| 1264 | "gz="<<xyz[2]<< // global position at reference |
| 1265 | "dRrec="<<dRrec<< // delta Radius in reconstruction |
| 1266 | "pt="<<pt<< //pt |
| 1267 | "id="<<id<< // track id |
| 1268 | "entries="<<entries;// number of entries in bin |
| 1269 | // |
| 1270 | AliExternalTrackParam *trackOut0 = 0; |
| 1271 | AliExternalTrackParam *trackOut1 = 0; |
| 1272 | AliExternalTrackParam *ptrackIn0 = 0; |
| 1273 | AliExternalTrackParam *ptrackIn1 = 0; |
| 1274 | |
| 1275 | for (Int_t icorr=0; icorr<ncorr; icorr++) { |
| 1276 | // |
| 1277 | // special case of the TPC tracks crossing the CE |
| 1278 | // |
| 1279 | AliTPCCorrection *corr = (AliTPCCorrection*)corrArray->At(icorr); |
| 1280 | corrections[icorr]=0; |
| 1281 | if (entries>kMinEntries &&isOK){ |
| 1282 | AliExternalTrackParam trackIn0(refX,phi,tPar,cov); |
| 1283 | AliExternalTrackParam trackIn1(refX,phi,tPar,cov); |
| 1284 | ptrackIn1=&trackIn0; |
| 1285 | ptrackIn0=&trackIn1; |
| 1286 | // |
| 1287 | if (debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,pcstream); |
| 1288 | if (!debug) trackOut0=corr->FitDistortedTrack(trackIn0, refX, dir,0); |
| 1289 | if (debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,pcstream); |
| 1290 | if (!debug) trackOut1=corr->FitDistortedTrack(trackIn1, refX, -dir,0); |
| 1291 | // |
| 1292 | if (trackOut0 && trackOut1){ |
| 1293 | // |
| 1294 | if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kTRUE,kMaxSnp)) isOK=kFALSE; |
| 1295 | if (!AliTrackerBase::PropagateTrackTo(&trackIn0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; |
| 1296 | // rotate all tracks to the same frame |
| 1297 | if (!trackOut0->Rotate(trackIn0.GetAlpha())) isOK=kFALSE; |
| 1298 | if (!trackIn1.Rotate(trackIn0.GetAlpha())) isOK=kFALSE; |
| 1299 | if (!trackOut1->Rotate(trackIn0.GetAlpha())) isOK=kFALSE; |
| 1300 | // |
| 1301 | if (!AliTrackerBase::PropagateTrackTo(trackOut0,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; |
| 1302 | if (!AliTrackerBase::PropagateTrackTo(&trackIn1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; |
| 1303 | if (!AliTrackerBase::PropagateTrackTo(trackOut1,refX,kMass,1,kFALSE,kMaxSnp)) isOK=kFALSE; |
| 1304 | // |
| 1305 | corrections[icorr] = (trackOut0->GetParameter()[ptype]-trackIn0.GetParameter()[ptype]); |
| 1306 | corrections[icorr]-= (trackOut1->GetParameter()[ptype]-trackIn1.GetParameter()[ptype]); |
| 1307 | (*pcstream)<<"fitDebug"<< // just to debug the correction |
| 1308 | "mean="<<mean<< |
| 1309 | "pIn0.="<<ptrackIn0<< |
| 1310 | "pIn1.="<<ptrackIn1<< |
| 1311 | "pOut0.="<<trackOut0<< |
| 1312 | "pOut1.="<<trackOut1<< |
| 1313 | "refX="<<refXD<< |
| 1314 | "\n"; |
| 1315 | delete trackOut0; |
| 1316 | delete trackOut1; |
| 1317 | }else{ |
| 1318 | corrections[icorr]=0; |
| 1319 | isOK=kFALSE; |
| 1320 | } |
| 1321 | } |
| 1322 | (*pcstream)<<"fit"<< |
| 1323 | Form("%s=",corr->GetName())<<corrections[icorr]; // dump correction value |
| 1324 | } |
| 1325 | // |
| 1326 | (*pcstream)<<"fit"<<"isOK="<<isOK<<"\n"; |
| 1327 | } |
| 1328 | delete pcstream; |
| 1329 | } |
| 1330 | |
| 1331 | |
git://git.uio.no / u / mrichter / AliRoot.git / blame