]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibCosmic.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
M AliTPCcalibAlign.cxx - filtering the debug output
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibCosmic.cxx
CommitLineData
f7f33dec 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
54b76c13 16/*
108953e9 17 Comments to be written here:
54b76c13 18 1. What do we calibrate.
19 2. How to interpret results
20 3. Simple example
21 4. Analysis using debug streamers.
22
23
24
25 3.Simple example
26 // To make cosmic scan the user interaction neccessary
27 //
28 .x ~/UliStyle.C
29 gSystem->Load("libANALYSIS");
30 gSystem->Load("libTPCcalib");
31 TFile fcalib("CalibObjects.root");
32 TObjArray * array = (TObjArray*)fcalib.Get("TPCCalib");
33 AliTPCcalibCosmic * cosmic = ( AliTPCcalibCosmic *)array->FindObject("cosmicTPC");
34
35
3e55050f 36 //
37 //
38 gSystem->AddIncludePath("-I$ALICE_ROOT/TPC/macros");
39 gROOT->LoadMacro("$ALICE_ROOT/TPC/macros/AliXRDPROOFtoolkit.cxx+")
40 AliXRDPROOFtoolkit tool;
a390f11f 41 TChain * chainCosmic = tool.MakeChainRandom("cosmic.txt","Track0",0,10000);
42 TChain * chainBudget = tool.MakeChainRandom("cosmic.txt","material",0,1000);
43 TCut cutptV="abs(1/pt0V-1/pt1V)<0.1";
44 TCut cutptI="abs(1/pt0In-1/pt1In)<0.5";
45 TCut cutncl="nclmin>120";
46 TCut cutDz="abs(p0.fP[1])<50";
47 TCut cutDr="abs(p0.fP[0])<50";
48 //
49 chainBudget->Draw(">>listB",cutptV+cutptI+cutncl+cutDr+cutDz,"entryList");
50 TEntryList *elistB = (TEntryList*)gDirectory->Get("listB");
51 chainBudget->SetEntryList(elistB);
52
53 chainBudget->SetAlias("dptrel","(pt0V-pt1V)/((pt0V+pt1V)*0.5)");
54 chainBudget->SetAlias("dptInrel","(pt0In-pt1In)/((pt0In+pt1In)*0.5)");
55 chainBudget->SetAlias("ptcorr","(pt0In-pt0V)/(pt0V)+(pt1V-pt1In)/(pt1In)");
54b76c13 56*/
57
58
59
f7f33dec 60#include "Riostream.h"
61#include "TChain.h"
62#include "TTree.h"
63#include "TH1F.h"
64#include "TH2F.h"
65#include "TList.h"
66#include "TMath.h"
67#include "TCanvas.h"
68#include "TFile.h"
54b76c13 69#include "TF1.h"
91fd44c9 70#include "THnSparse.h"
9963b5e2 71#include "TDatabasePDG.h"
f7f33dec 72
54b76c13 73#include "AliTPCclusterMI.h"
f7f33dec 74#include "AliTPCseed.h"
75#include "AliESDVertex.h"
76#include "AliESDEvent.h"
77#include "AliESDfriend.h"
78#include "AliESDInputHandler.h"
54b76c13 79#include "AliAnalysisManager.h"
f7f33dec 80
81#include "AliTracker.h"
f7a1cc68 82#include "AliMagF.h"
54b76c13 83#include "AliTPCCalROC.h"
f7f33dec 84
85#include "AliLog.h"
86
87#include "AliTPCcalibCosmic.h"
f7f33dec 88#include "TTreeStream.h"
89#include "AliTPCTracklet.h"
3326b323 90//#include "AliESDcosmic.h"
15e48021 91
f7f33dec 92
93ClassImp(AliTPCcalibCosmic)
94
95
96AliTPCcalibCosmic::AliTPCcalibCosmic()
97 :AliTPCcalibBase(),
9b27d39b 98 fHistNTracks(0),
99 fClusters(0),
100 fModules(0),
101 fHistPt(0),
9b27d39b 102 fDeDx(0),
54b76c13 103 fDeDxMIP(0),
104 fMIPvalue(1),
9b27d39b 105 fCutMaxD(5), // maximal distance in rfi ditection
a6dc0cf6 106 fCutMaxDz(40), // maximal distance in z ditection
9b27d39b 107 fCutTheta(0.03), // maximal distan theta
108 fCutMinDir(-0.99) // direction vector products
f7f33dec 109{
54b76c13 110 AliInfo("Default Constructor");
91fd44c9 111 for (Int_t ihis=0; ihis<6;ihis++){
112 fHistoDelta[ihis]=0;
113 fHistoPull[ihis]=0;
8a92e133 114 }
115 for (Int_t ihis=0; ihis<4;ihis++){
116 fHistodEdxMax[ihis] =0;
117 fHistodEdxTot[ihis] =0;
91fd44c9 118 }
f7f33dec 119}
120
121
122AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
123 :AliTPCcalibBase(),
9b27d39b 124 fHistNTracks(0),
125 fClusters(0),
126 fModules(0),
127 fHistPt(0),
9b27d39b 128 fDeDx(0),
54b76c13 129 fDeDxMIP(0),
130 fMIPvalue(1),
a6dc0cf6 131 fCutMaxD(5), // maximal distance in rfi ditection
132 fCutMaxDz(40), // maximal distance in z ditection
9b27d39b 133 fCutTheta(0.03), // maximal distan theta
134 fCutMinDir(-0.99) // direction vector products
f7f33dec 135{
136 SetName(name);
137 SetTitle(title);
54b76c13 138
139 fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
140 fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
141 fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
142 fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
143 fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
f7f33dec 144 BinLogX(fDeDx);
54b76c13 145 fDeDxMIP = new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
91fd44c9 146 Init();
9b27d39b 147 AliInfo("Non Default Constructor");
bca20570 148 //
f7f33dec 149}
150
151AliTPCcalibCosmic::~AliTPCcalibCosmic(){
152 //
153 //
154 //
91fd44c9 155 for (Int_t ihis=0; ihis<6;ihis++){
156 delete fHistoDelta[ihis];
157 delete fHistoPull[ihis];
91fd44c9 158 }
8a92e133 159 for (Int_t ihis=0; ihis<4;ihis++){
160 delete fHistodEdxTot[ihis];
161 delete fHistodEdxMax[ihis];
162 }
163
164 delete fHistNTracks; // histogram showing number of ESD tracks per event
165 delete fClusters; // histogram showing the number of clusters per track
166 delete fModules; // 2d histogram of tracks which are propagated to the ACORDE scintillator array
167 delete fHistPt; // Pt histogram of reconstructed tracks
168 delete fDeDx; // dEdx spectrum showing the different particle types
169 delete fDeDxMIP; // TPC signal close to the MIP region of muons 0.4 < p < 0.45 GeV
f7f33dec 170}
171
172
91fd44c9 173void AliTPCcalibCosmic::Init(){
174 //
175 // init component
176 // Make performance histograms
177 //
178
179 // tracking performance bins
180 // 0 - delta of interest
181 // 1 - min (track0, track1) number of clusters
182 // 2 - R - vertex radius
183 // 3 - P1 - mean z
184 // 4 - P2 - snp(phi) at inner wall of TPC
185 // 5 - P3 - tan(theta) at inner wall of TPC
186 // 6 - P4 - 1/pt mean
187 // 7 - pt - pt mean
188 // 8 - alpha
189
190 Double_t xminTrack[9], xmaxTrack[9];
191 Int_t binsTrack[9];
192 TString axisName[9];
193 //
194 binsTrack[0] =100;
195 axisName[0] ="#Delta";
196 //
197 binsTrack[1] =8;
198 xminTrack[1] =80; xmaxTrack[1]=160;
199 axisName[1] ="N_{cl}";
200 //
201 binsTrack[2] =10;
202 xminTrack[2] =0; xmaxTrack[2]=90; //
203 axisName[2] ="dca_{r} (cm)";
204 //
205 binsTrack[3] =25;
206 xminTrack[3] =-250; xmaxTrack[3]=250; //
207 axisName[3] ="z (cm)";
208 //
209 binsTrack[4] =10;
210 xminTrack[4] =-0.8; xmaxTrack[4]=0.8; //
211 axisName[4] ="sin(#phi)";
212 //
213 binsTrack[5] =10;
214 xminTrack[5] =-1; xmaxTrack[5]=1; //
8a92e133 215 axisName[5] ="tan(#theta)";
91fd44c9 216 //
a390f11f 217 binsTrack[6] =40;
218 xminTrack[6] =-2; xmaxTrack[6]=2; //
91fd44c9 219 axisName[6] ="1/pt (1/GeV)";
220 //
8a92e133 221 binsTrack[7] =40;
91fd44c9 222 xminTrack[7] =0.2; xmaxTrack[7]=50; //
8a92e133 223 axisName[7] ="pt (GeV)";
91fd44c9 224 //
a390f11f 225 binsTrack[8] =18;
91fd44c9 226 xminTrack[8] =0; xmaxTrack[8]=TMath::Pi(); //
227 axisName[8] ="alpha";
228 //
229 // delta y
230 xminTrack[0] =-1; xmaxTrack[0]=1; //
231 fHistoDelta[0] = new THnSparseS("#Delta_{Y} (cm)","#Delta_{Y} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
232 xminTrack[0] =-5; xmaxTrack[0]=5; //
233 fHistoPull[0] = new THnSparseS("#Delta_{Y} (unit)","#Delta_{Y} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
234 //
235 // delta z
236 xminTrack[0] =-1; xmaxTrack[0]=1; //
237 fHistoDelta[1] = new THnSparseS("#Delta_{Z} (cm)","#Delta_{Z} (cm)", 9, binsTrack,xminTrack, xmaxTrack);
238 xminTrack[0] =-5; xmaxTrack[0]=5; //
239 fHistoPull[1] = new THnSparseS("#Delta_{Z} (unit)","#Delta_{Z} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
240 //
241 // delta P2
242 xminTrack[0] =-10; xmaxTrack[0]=10; //
243 fHistoDelta[2] = new THnSparseS("#Delta_{#phi} (mrad)","#Delta_{#phi} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
244 xminTrack[0] =-5; xmaxTrack[0]=5; //
245 fHistoPull[2] = new THnSparseS("#Delta_{#phi} (unit)","#Delta_{#phi} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
246 //
247 // delta P3
248 xminTrack[0] =-10; xmaxTrack[0]=10; //
249 fHistoDelta[3] = new THnSparseS("#Delta_{#theta} (mrad)","#Delta_{#theta} (mrad)", 9, binsTrack,xminTrack, xmaxTrack);
250 xminTrack[0] =-5; xmaxTrack[0]=5; //
251 fHistoPull[3] = new THnSparseS("#Delta_{#theta} (unit)","#Delta_{#theta} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
252 //
253 // delta P4
254 xminTrack[0] =-0.2; xmaxTrack[0]=0.2; //
255 fHistoDelta[4] = new THnSparseS("#Delta_{1/pt} (1/GeV)","#Delta_{1/pt} (1/GeV)", 9, binsTrack,xminTrack, xmaxTrack);
256 xminTrack[0] =-5; xmaxTrack[0]=5; //
257 fHistoPull[4] = new THnSparseS("#Delta_{1/pt} (unit)","#Delta_{1/pt} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
258
259 //
260 // delta Pt
261 xminTrack[0] =-0.5; xmaxTrack[0]=0.5; //
262 fHistoDelta[5] = new THnSparseS("#Delta_{pt}/p_{t}","#Delta_{pt}/p_{t}", 9, binsTrack,xminTrack, xmaxTrack);
263 xminTrack[0] =-5; xmaxTrack[0]=5; //
264 fHistoPull[5] = new THnSparseS("#Delta_{pt}/p_{t} (unit)","#Delta_{pt}/p_{t} (unit)", 9, binsTrack,xminTrack, xmaxTrack);
265 //
8a92e133 266
267 for (Int_t idedx=0;idedx<4;idedx++){
268 xminTrack[0] =0.5; xmaxTrack[0]=1.5; //
269 binsTrack[1] =40;
270 xminTrack[1] =10; xmaxTrack[1]=160;
271
272 fHistodEdxMax[idedx] = new THnSparseS(Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
273 Form("dEdx_{MaxUp}/dEdx_{MaxDown}_Pad%d",idedx),
274 9, binsTrack,xminTrack, xmaxTrack);
275 fHistodEdxTot[idedx] = new THnSparseS(Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
276 Form("dEdx_{TotUp}/dEdx_{TotDown}_Pad%d",idedx),
277 9, binsTrack,xminTrack, xmaxTrack);
278 }
279
280
281
91fd44c9 282 for (Int_t ivar=0;ivar<6;ivar++){
8a92e133 283 for (Int_t ivar2=0;ivar2<9;ivar2++){
91fd44c9 284 fHistoDelta[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
285 fHistoDelta[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
286 fHistoPull[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
287 fHistoPull[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
288 BinLogX(fHistoDelta[ivar],7);
289 BinLogX(fHistoPull[ivar],7);
8a92e133 290 if (ivar<4){
291 fHistodEdxMax[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
292 fHistodEdxMax[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
293 fHistodEdxTot[ivar]->GetAxis(ivar2)->SetName(axisName[ivar2].Data());
294 fHistodEdxTot[ivar]->GetAxis(ivar2)->SetTitle(axisName[ivar2].Data());
295 BinLogX(fHistodEdxMax[ivar],7);
296 BinLogX(fHistodEdxTot[ivar],7);
297 }
91fd44c9 298 }
299 }
300}
301
302
303void AliTPCcalibCosmic::Add(const AliTPCcalibCosmic* cosmic){
304 //
305 //
306 //
307 for (Int_t ivar=0; ivar<6;ivar++){
308 if (fHistoDelta[ivar] && cosmic->fHistoDelta[ivar]){
309 fHistoDelta[ivar]->Add(cosmic->fHistoDelta[ivar]);
310 }
311 if (fHistoPull[ivar] && cosmic->fHistoPull[ivar]){
312 fHistoPull[ivar]->Add(cosmic->fHistoPull[ivar]);
313 }
314 }
8a92e133 315 for (Int_t ivar=0; ivar<4;ivar++){
316 if (fHistodEdxMax[ivar] && cosmic->fHistodEdxMax[ivar]){
317 fHistodEdxMax[ivar]->Add(cosmic->fHistodEdxMax[ivar]);
318 }
319 if (fHistodEdxTot[ivar] && cosmic->fHistodEdxTot[ivar]){
320 fHistodEdxTot[ivar]->Add(cosmic->fHistodEdxTot[ivar]);
321 }
322 }
91fd44c9 323}
324
9b27d39b 325
326
327
f7f33dec 328void AliTPCcalibCosmic::Process(AliESDEvent *event) {
9b27d39b 329 //
330 //
331 //
f7f33dec 332 if (!event) {
333 Printf("ERROR: ESD not available");
334 return;
9b27d39b 335 }
76273318 336 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
337 if (!esdFriend) {
338 Printf("ERROR: esdFriend not available");
f7f33dec 339 return;
340 }
2acad464 341
108953e9 342
54b76c13 343 FindPairs(event); // nearly everything takes place in find pairs...
344
2acad464 345 if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
9b27d39b 346 Int_t ntracks=event->GetNumberOfTracks();
347 fHistNTracks->Fill(ntracks);
9b27d39b 348 if (ntracks==0) return;
cc65e4f5 349 // AliESDcosmic cosmicESD;
350// TTreeSRedirector * cstream = GetDebugStreamer();
351// cosmicESD.SetDebugStreamer(cstream);
352// cosmicESD.ProcessEvent(event);
353// if (cstream) cosmicESD.DumpToTree();
15e48021 354
355
54b76c13 356}
9b27d39b 357
f7f33dec 358
3326b323 359void AliTPCcalibCosmic::FillHistoPerformance(const AliExternalTrackParam *par0, const AliExternalTrackParam *par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam */*inner1*/, AliTPCseed *seed0, AliTPCseed *seed1, const AliExternalTrackParam *param0Combined ){
91fd44c9 360 //
a390f11f 361 // par0,par1 - parameter of tracks at DCA 0
362 // inner0,inner1 - parameter of tracks at the TPC entrance
363 // seed0, seed1 - detailed track information
364 // param0Combined - Use combined track parameters for binning
91fd44c9 365 //
8a92e133 366 Int_t kMinCldEdx =20;
367 Int_t ncl0 = seed0->GetNumberOfClusters();
368 Int_t ncl1 = seed1->GetNumberOfClusters();
369
91fd44c9 370 const Double_t kpullCut = 10;
371 Double_t x[9];
372 Double_t xyz0[3];
373 Double_t xyz1[3];
374 par0->GetXYZ(xyz0);
375 par1->GetXYZ(xyz1);
376 Double_t radius0 = TMath::Sqrt(xyz0[0]*xyz0[0]+xyz0[1]*xyz0[1]);
377 Double_t radius1 = TMath::Sqrt(xyz1[0]*xyz1[0]+xyz1[1]*xyz1[1]);
378 inner0->GetXYZ(xyz0);
379 Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
380 // bin parameters
381 x[1] = TMath::Min(ncl0,ncl1);
382 x[2] = (radius0+radius1)*0.5;
a390f11f 383 x[3] = param0Combined->GetZ();
384 x[4] = inner0->GetSnp();
385 x[5] = param0Combined->GetTgl();
386 x[6] = param0Combined->GetSigned1Pt();
387 x[7] = param0Combined->Pt();
91fd44c9 388 x[8] = alpha;
389 // deltas
390 Double_t delta[6];
391 Double_t sigma[6];
392 delta[0] = (par0->GetY()+par1->GetY());
393 delta[1] = (par0->GetZ()-par1->GetZ());
394 delta[2] = (par0->GetAlpha()-par1->GetAlpha()-TMath::Pi());
395 delta[3] = (par0->GetTgl()+par1->GetTgl());
396 delta[4] = (par0->GetParameter()[4]+par1->GetParameter()[4]);
397 delta[5] = (par0->Pt()-par1->Pt())/((par0->Pt()+par1->Pt())*0.5);
398 //
399 sigma[0] = TMath::Sqrt(par0->GetSigmaY2()+par1->GetSigmaY2());
400 sigma[1] = TMath::Sqrt(par0->GetSigmaZ2()+par1->GetSigmaZ2());
401 sigma[2] = TMath::Sqrt(par0->GetSigmaSnp2()+par1->GetSigmaSnp2());
402 sigma[3] = TMath::Sqrt(par0->GetSigmaTgl2()+par1->GetSigmaTgl2());
403 sigma[4] = TMath::Sqrt(par0->GetSigma1Pt2()+par1->GetSigma1Pt2());
404 sigma[5] = sigma[4]*((par0->Pt()+par1->Pt())*0.5);
405 //
406 Bool_t isOK = kTRUE;
407 for (Int_t ivar=0;ivar<6;ivar++){
408 if (sigma[ivar]==0) isOK=kFALSE;
409 x[0]= delta[ivar]/sigma[ivar];
410 if (TMath::Abs(x[0])>kpullCut) isOK = kFALSE;
411 }
412 //
413
414 if (isOK) for (Int_t ivar=0;ivar<6;ivar++){
415 x[0]= delta[ivar]/TMath::Sqrt(2);
416 if (ivar==2 || ivar ==3) x[0]*=1000;
417 fHistoDelta[ivar]->Fill(x);
418 if (sigma[ivar]>0){
419 x[0]= delta[ivar]/sigma[ivar];
420 fHistoPull[ivar]->Fill(x);
421 }
422 }
8a92e133 423
424 //
425 // Fill dedx performance
426 //
427 for (Int_t ipad=0; ipad<4;ipad++){
428 //
429 //
430 //
431 Int_t row0=0;
432 Int_t row1=160;
433 if (ipad==0) row1=63;
434 if (ipad==1) {row0=63; row1=63+64;}
435 if (ipad==2) {row0=128;}
436 Int_t nclUp = TMath::Nint(seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
437 Int_t nclDown = TMath::Nint(seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1,2));
438 Int_t minCl = TMath::Min(nclUp,nclDown);
439 if (minCl<kMinCldEdx) continue;
440 x[1] = minCl;
441 //
442 Float_t dEdxTotUp = seed0->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
443 Float_t dEdxTotDown = seed1->CookdEdxAnalytical(0.01,0.7,0,row0,row1);
444 Float_t dEdxMaxUp = seed0->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
445 Float_t dEdxMaxDown = seed1->CookdEdxAnalytical(0.01,0.7,1,row0,row1);
446 //
447 if (dEdxTotDown<=0) continue;
448 if (dEdxMaxDown<=0) continue;
449 x[0]=dEdxTotUp/dEdxTotDown;
450 fHistodEdxTot[ipad]->Fill(x);
451 x[0]=dEdxMaxUp/dEdxMaxDown;
452 fHistodEdxMax[ipad]->Fill(x);
453 }
454
455
91fd44c9 456
457}
458
3326b323 459void AliTPCcalibCosmic::MaterialBudgetDump(AliExternalTrackParam *const par0, AliExternalTrackParam *const par1, const AliExternalTrackParam *inner0, const AliExternalTrackParam *inner1, AliTPCseed *const seed0, AliTPCseed *const seed1, AliExternalTrackParam *const param0Combined, AliExternalTrackParam *const param1Combined){
a390f11f 460 //
461 // matrial budget AOD dump
462 //
463 // par0,par1 - parameter of tracks at DCA 0
464 // inner0,inner1 - parameter of tracks at the TPC entrance
465 // seed0, seed1 - detailed track information
466 // param0Combined - Use combined track parameters for binning
467 // param1Combined -
468 Double_t p0In = inner0->GetP();
469 Double_t p1In = inner1->GetP();
470 Double_t p0V = par0->GetP();
471 Double_t p1V = par1->GetP();
472 //
473 Double_t pt0In = inner0->Pt();
474 Double_t pt1In = inner1->Pt();
475 Double_t pt0V = par0->Pt();
476 Double_t pt1V = par1->Pt();
477 Int_t ncl0 = seed0->GetNumberOfClusters();
478 Int_t ncl1 = seed1->GetNumberOfClusters();
479 Int_t nclmin=TMath::Min(ncl0,ncl1);
480 Double_t sign = (param0Combined->GetSigned1Pt()>0) ? 1:-1.;
481 //
482 TTreeSRedirector * pcstream = GetDebugStreamer();
483 if (pcstream){
484 (*pcstream)<<"material"<<
485 "run="<<fRun<< // run number
486 "event="<<fEvent<< // event number
487 "time="<<fTime<< // time stamp of event
488 "trigger="<<fTrigger<< // trigger
489 "triggerClass="<<&fTriggerClass<< // trigger
490 "mag="<<fMagF<< // magnetic field
491 "sign="<<sign<< // sign of the track
492 //
493 "ncl0="<<ncl0<<
494 "ncl1="<<ncl1<<
495 "nclmin="<<nclmin<<
496 //
497 "p0In="<<p0In<<
498 "p1In="<<p1In<<
499 "p0V="<<p0V<<
500 "p1V="<<p1V<<
501 "pt0In="<<pt0In<<
502 "pt1In="<<pt1In<<
503 "pt0V="<<pt0V<<
504 "pt1V="<<pt1V<<
505 "p0.="<<par0<<
506 "p1.="<<par1<<
507 "up0.="<<param0Combined<<
508 "up1.="<<param1Combined<<
509 "\n";
510 }
511
512}
f7f33dec 513
54b76c13 514void AliTPCcalibCosmic::Analyze() {
515
516 fMIPvalue = CalculateMIPvalue(fDeDxMIP);
517
518 return;
519
520}
521
f7f33dec 522
f7f33dec 523
9b27d39b 524void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
525 //
526 // Find cosmic pairs
04087794 527 //
528 // Track0 is choosen in upper TPC part
529 // Track1 is choosen in lower TPC part
9b27d39b 530 //
2acad464 531 if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
76273318 532 AliESDfriend *esdFriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
9b27d39b 533 Int_t ntracks=event->GetNumberOfTracks();
9b27d39b 534 TObjArray tpcSeeds(ntracks);
535 if (ntracks==0) return;
536 Double_t vtxx[3]={0,0,0};
537 Double_t svtxx[3]={0.000001,0.000001,100.};
538 AliESDVertex vtx(vtxx,svtxx);
539 //
540 //track loop
541 //
54b76c13 542 for (Int_t i=0;i<ntracks;++i) {
543 AliESDtrack *track = event->GetTrack(i);
544 fClusters->Fill(track->GetTPCNcls());
545
546 const AliExternalTrackParam * trackIn = track->GetInnerParam();
547 const AliExternalTrackParam * trackOut = track->GetOuterParam();
548 if (!trackIn) continue;
549 if (!trackOut) continue;
860b3d93 550 if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // filter laser
551
552
76273318 553 AliESDfriendTrack *friendTrack = esdFriend->GetTrack(i);
9b27d39b 554 TObject *calibObject;
555 AliTPCseed *seed = 0;
556 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
557 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
558 }
559 if (seed) tpcSeeds.AddAt(seed,i);
54b76c13 560
561 Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
562 if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
15e48021 563 fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159));
54b76c13 564 //
15e48021 565 if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159));
54b76c13 566 //
15e48021 567 if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) {
54b76c13 568 TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
b322e06a 569 //if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
570 // if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
54b76c13 571 }
572
573 }
574
9b27d39b 575 }
54b76c13 576
9b27d39b 577 if (ntracks<2) return;
578 //
579 // Find pairs
580 //
581 for (Int_t i=0;i<ntracks;++i) {
582 AliESDtrack *track0 = event->GetTrack(i);
04087794 583 // track0 - choosen upper part
584 if (!track0) continue;
585 if (!track0->GetOuterParam()) continue;
586 if (track0->GetOuterParam()->GetAlpha()<0) continue;
e9362f9d 587 Double_t dir0[3];
588 track0->GetDirection(dir0);
04087794 589 for (Int_t j=0;j<ntracks;++j) {
590 if (i==j) continue;
591 AliESDtrack *track1 = event->GetTrack(j);
592 //track 1 lower part
593 if (!track1) continue;
594 if (!track1->GetOuterParam()) continue;
595 if (track1->GetOuterParam()->GetAlpha()>0) continue;
596 //
e9362f9d 597 Double_t dir1[3];
598 track1->GetDirection(dir1);
54b76c13 599
9b27d39b 600 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
601 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
602 if (! seed0) continue;
603 if (! seed1) continue;
15e48021 604 Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159);
605 Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159);
bca20570 606 //
15e48021 607 Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
608 Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
bca20570 609 //
15e48021 610 Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159);
611 Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159);
bca20570 612 //
e9362f9d 613 Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
9b27d39b 614 Float_t d0 = track0->GetLinearD(0,0);
615 Float_t d1 = track1->GetLinearD(0,0);
616 //
617 // conservative cuts - convergence to be guarantied
618 // applying before track propagation
619 if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
620 if (dir>fCutMinDir) continue; // direction vector product
621 Float_t bz = AliTracker::GetBz();
622 Float_t dvertex0[2]; //distance to 0,0
623 Float_t dvertex1[2]; //distance to 0,0
624 track0->GetDZ(0,0,0,bz,dvertex0);
625 track1->GetDZ(0,0,0,bz,dvertex1);
626 if (TMath::Abs(dvertex0[1])>250) continue;
627 if (TMath::Abs(dvertex1[1])>250) continue;
628 //
629 //
630 //
631 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
632 AliExternalTrackParam param0(*track0);
633 AliExternalTrackParam param1(*track1);
634 //
635 // Propagate using Magnetic field and correct fo material budget
636 //
a390f11f 637 Double_t sign0=-1;
638 Double_t sign1=1;
639 Double_t maxsnp=0.90;
9963b5e2 640 AliTracker::PropagateTrackToBxByBz(&param0,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign0);
641 AliTracker::PropagateTrackToBxByBz(&param1,dmax+1,TDatabasePDG::Instance()->GetParticle("e-")->Mass(),3,kTRUE,maxsnp,sign1);
9b27d39b 642 //
643 // Propagate rest to the 0,0 DCA - z should be ignored
644 //
645 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
646 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
647 //
648 param0.GetDZ(0,0,0,bz,dvertex0);
649 param1.GetDZ(0,0,0,bz,dvertex1);
a6dc0cf6 650 if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
9b27d39b 651 //
652 Double_t xyz0[3];//,pxyz0[3];
653 Double_t xyz1[3];//,pxyz1[3];
654 param0.GetXYZ(xyz0);
655 param1.GetXYZ(xyz1);
656 Bool_t isPair = IsPair(&param0,&param1);
657 //
54b76c13 658 if (isPair) FillAcordeHist(track0);
15e48021 659 //
660 // combined track params
661 //
662 AliExternalTrackParam *par0U=MakeCombinedTrack(&param0,&param1);
663 AliExternalTrackParam *par1U=MakeCombinedTrack(&param1,&param0);
664
665
54b76c13 666 //
9b27d39b 667 if (fStreamLevel>0){
668 TTreeSRedirector * cstream = GetDebugStreamer();
108953e9 669 //printf("My stream=%p\n",(void*)cstream);
e9362f9d 670 AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
671 AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
672 AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
673 AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
674 Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
675 Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
676 Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
677 Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
91fd44c9 678 //
679 //
680 //
a390f11f 681 FillHistoPerformance(&param0, &param1, ip0, ip1, seed0, seed1,par0U);
682 MaterialBudgetDump(&param0, &param1, ip0, ip1, seed0, seed1,par0U,par1U);
9b27d39b 683 if (cstream) {
684 (*cstream) << "Track0" <<
108953e9 685 "run="<<fRun<< // run number
686 "event="<<fEvent<< // event number
687 "time="<<fTime<< // time stamp of event
688 "trigger="<<fTrigger<< // trigger
15e48021 689 "triggerClass="<<&fTriggerClass<< // trigger
108953e9 690 "mag="<<fMagF<< // magnetic field
9b27d39b 691 "dir="<<dir<< // direction
54b76c13 692 "OK="<<isPair<< // will be accepted
9b27d39b 693 "b0="<<b0<< // propagate status
694 "b1="<<b1<< // propagate status
e9362f9d 695 "crossI="<<isCrossI<< // cross inner
696 "crossO="<<isCrossO<< // cross outer
697 //
9b27d39b 698 "Orig0.=" << track0 << // original track 0
699 "Orig1.=" << track1 << // original track 1
700 "Tr0.="<<&param0<< // track propagated to the DCA 0,0
701 "Tr1.="<<&param1<< // track propagated to the DCA 0,0
e9362f9d 702 "Ip0.="<<ip0<< // inner param - upper
703 "Ip1.="<<ip1<< // inner param - lower
704 "Op0.="<<op0<< // outer param - upper
705 "Op1.="<<op1<< // outer param - lower
15e48021 706 "Up0.="<<par0U<< // combined track 0
707 "Up1.="<<par1U<< // combined track 1
e9362f9d 708 //
9b27d39b 709 "v00="<<dvertex0[0]<< // distance using kalman
710 "v01="<<dvertex0[1]<< //
711 "v10="<<dvertex1[0]<< //
712 "v11="<<dvertex1[1]<< //
713 "d0="<<d0<< // linear distance to 0,0
714 "d1="<<d1<< // linear distance to 0,0
715 //
e9362f9d 716 //
717 //
718 "x00="<<xyz0[0]<< // global position close to vertex
9b27d39b 719 "x01="<<xyz0[1]<<
720 "x02="<<xyz0[2]<<
721 //
e9362f9d 722 "x10="<<xyz1[0]<< // global position close to vertex
9b27d39b 723 "x11="<<xyz1[1]<<
724 "x12="<<xyz1[2]<<
725 //
e9362f9d 726 "alpha0="<<alpha0<<
727 "alpha1="<<alpha1<<
728 "dir00="<<dir0[0]<< // direction upper
729 "dir01="<<dir0[1]<<
730 "dir02="<<dir0[2]<<
731 //
732 "dir10="<<dir1[0]<< // direction lower
733 "dir11="<<dir1[1]<<
734 "dir12="<<dir1[2]<<
735 //
736 //
54b76c13 737 "Seed0.=" << seed0 << // original seed 0
738 "Seed1.=" << seed1 << // original seed 1
bca20570 739 //
740 "dedx0="<<dedx0<< // dedx0 - all
741 "dedx1="<<dedx1<< // dedx1 - all
742 //
54b76c13 743 "dedx0I="<<dedx0I<< // dedx0 - inner ROC
744 "dedx1I="<<dedx1I<< // dedx1 - inner ROC
bca20570 745 //
54b76c13 746 "dedx0O="<<dedx0O<< // dedx0 - outer ROC
747 "dedx1O="<<dedx1O<< // dedx1 - outer ROC
9b27d39b 748 "\n";
749 }
15e48021 750 }
751 delete par0U;
752 delete par1U;
9b27d39b 753 }
754 }
755}
756
9b27d39b 757
bca20570 758
759
54b76c13 760void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
761
762 // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
763 if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
764
76273318 765 const Double_t acordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
54b76c13 766 const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
767
768 Double_t r[3];
769 upperTrack->GetXYZ(r);
770 Double_t d[3];
771 upperTrack->GetDirection(d);
772 Double_t x,z;
76273318 773 z = r[2] + (d[2]/d[1])*(acordePlane - r[1]);
774 x = r[0] + (d[0]/d[1])*(acordePlane - r[1]);
54b76c13 775
776 if (x > roof) {
76273318 777 x = r[0] + (d[0]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
778 z = r[2] + (d[2]/(d[0]+d[1]))*(acordePlane+roof-r[0]-r[1]);
54b76c13 779 }
780 if (x < -roof) {
76273318 781 x = r[0] + (d[0]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
782 z = r[2] + (d[2]/(d[1]-d[0]))*(acordePlane+roof+r[0]-r[1]);
54b76c13 783 }
784
785 fModules->Fill(z, x);
786
787}
788
789
790
3326b323 791Long64_t AliTPCcalibCosmic::Merge(TCollection *const li) {
bca20570 792
793 TIterator* iter = li->MakeIterator();
794 AliTPCcalibCosmic* cal = 0;
795
796 while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
797 if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
860b3d93 798 //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
bca20570 799 return -1;
800 }
801
54b76c13 802 fHistNTracks->Add(cal->GetHistNTracks());
803 fClusters->Add(cal-> GetHistClusters());
804 fModules->Add(cal->GetHistAcorde());
805 fHistPt->Add(cal->GetHistPt());
806 fDeDx->Add(cal->GetHistDeDx());
807 fDeDxMIP->Add(cal->GetHistMIP());
91fd44c9 808 Add(cal);
bca20570 809 }
bca20570 810 return 0;
9b27d39b 811
f7f33dec 812}
813
54b76c13 814
9b27d39b 815Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
816 //
817 //
818 /*
819 // 0. Same direction - OPOSITE - cutDir +cutT
820 TCut cutDir("cutDir","dir<-0.99")
821 // 1.
822 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
823 //
824 // 2. The same rphi
825 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
826 //
827 //
828 //
829 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
830 // 1/Pt diff cut
831 */
832 const Double_t *p0 = tr0->GetParameter();
833 const Double_t *p1 = tr1->GetParameter();
834 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
a6dc0cf6 835 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
9b27d39b 836 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
a6dc0cf6 837
9b27d39b 838 Double_t d0[3], d1[3];
839 tr0->GetDirection(d0);
840 tr1->GetDirection(d1);
841 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
842 //
843 return kTRUE;
844}
54b76c13 845
846
847
848Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
849
850 TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
851 funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
852 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
853 hist->Fit(funcDoubleGaus);
3326b323 854 Double_t aMIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
54b76c13 855
856 delete funcDoubleGaus;
857
3326b323 858 return aMIPvalue;
54b76c13 859
860}
9b27d39b 861
862
f7f33dec 863
54b76c13 864
57dc06f2 865void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
866 //
867 // Not implemented yet
868 //
54b76c13 869 return;
870
871}
872
873
3326b323 874void AliTPCcalibCosmic::BinLogX(THnSparse *const h, Int_t axisDim) {
91fd44c9 875
876 // Method for the correct logarithmic binning of histograms
877
878 TAxis *axis = h->GetAxis(axisDim);
879 int bins = axis->GetNbins();
880
881 Double_t from = axis->GetXmin();
882 Double_t to = axis->GetXmax();
3326b323 883 Double_t *newBins = new Double_t[bins + 1];
91fd44c9 884
3326b323 885 newBins[0] = from;
91fd44c9 886 Double_t factor = pow(to/from, 1./bins);
887
888 for (int i = 1; i <= bins; i++) {
3326b323 889 newBins[i] = factor * newBins[i-1];
91fd44c9 890 }
3326b323 891 axis->Set(bins, newBins);
892 delete newBins;
91fd44c9 893
894}
895
896
3326b323 897void AliTPCcalibCosmic::BinLogX(TH1 *const h) {
f7f33dec 898
899 // Method for the correct logarithmic binning of histograms
900
901 TAxis *axis = h->GetXaxis();
902 int bins = axis->GetNbins();
903
904 Double_t from = axis->GetXmin();
905 Double_t to = axis->GetXmax();
3326b323 906 Double_t *newBins = new Double_t[bins + 1];
f7f33dec 907
3326b323 908 newBins[0] = from;
f7f33dec 909 Double_t factor = pow(to/from, 1./bins);
910
911 for (int i = 1; i <= bins; i++) {
3326b323 912 newBins[i] = factor * newBins[i-1];
f7f33dec 913 }
3326b323 914 axis->Set(bins, newBins);
915 delete newBins;
f7f33dec 916
917}
918
7b18d067 919
860b3d93 920AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
921 //
922 //
923 //
924 AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
925 par1R->Rotate(track0->GetAlpha());
15e48021 926 par1R->PropagateTo(track0->GetX(),AliTracker::GetBz());
860b3d93 927 //
928 //
929 Double_t * param = (Double_t*)par1R->GetParameter();
930 Double_t * covar = (Double_t*)par1R->GetCovariance();
15e48021 931
860b3d93 932 param[0]*=1; //OK
933 param[1]*=1; //OK
934 param[2]*=1; //?
935 param[3]*=-1; //OK
936 param[4]*=-1; //OK
937 //
938 covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
939 //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
940 covar[13]*=-1.;
860b3d93 941 return par1R;
942}
943
15e48021 944AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
945 //
946 // Make combined track
947 //
948 //
949 AliExternalTrackParam * par1T = MakeTrack(track0,track1);
950 AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
951 //
952 UpdateTrack(*par0U,*par1T);
953 delete par1T;
954 return par0U;
955}
956
957
860b3d93 958void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
959 //
960 // Update track 1 with track 2
961 //
962 //
963 //
964 TMatrixD vecXk(5,1); // X vector
965 TMatrixD covXk(5,5); // X covariance
966 TMatrixD matHk(5,5); // vector to mesurement
967 TMatrixD measR(5,5); // measurement error
968 TMatrixD vecZk(5,1); // measurement
969 //
970 TMatrixD vecYk(5,1); // Innovation or measurement residual
971 TMatrixD matHkT(5,5);
972 TMatrixD matSk(5,5); // Innovation (or residual) covariance
973 TMatrixD matKk(5,5); // Optimal Kalman gain
974 TMatrixD mat1(5,5); // update covariance matrix
975 TMatrixD covXk2(5,5); //
976 TMatrixD covOut(5,5);
977 //
978 Double_t *param1=(Double_t*) track1.GetParameter();
979 Double_t *covar1=(Double_t*) track1.GetCovariance();
980 Double_t *param2=(Double_t*) track2.GetParameter();
981 Double_t *covar2=(Double_t*) track2.GetCovariance();
982 //
983 // copy data to the matrix
984 for (Int_t ipar=0; ipar<5; ipar++){
860b3d93 985 for (Int_t jpar=0; jpar<5; jpar++){
986 covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
987 measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
15e48021 988 matHk(ipar,jpar)=0;
989 mat1(ipar,jpar)=0;
860b3d93 990 }
15e48021 991 vecXk(ipar,0) = param1[ipar];
992 vecZk(ipar,0) = param2[ipar];
993 matHk(ipar,ipar)=1;
994 mat1(ipar,ipar)=0;
860b3d93 995 }
996 //
997 //
998 //
999 //
860b3d93 1000 //
1001 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
1002 matHkT=matHk.T(); matHk.T();
1003 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
1004 matSk.Invert();
1005 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
1006 vecXk += matKk*vecYk; // updated vector
860b3d93 1007 covXk2 = (mat1-(matKk*matHk));
1008 covOut = covXk2*covXk;
1009 //
1010 //
1011 //
1012 // copy from matrix to parameters
1013 if (0) {
1014 vecXk.Print();
1015 vecZk.Print();
1016 //
1017 measR.Print();
1018 covXk.Print();
1019 covOut.Print();
1020 //
1021 track1.Print();
1022 track2.Print();
1023 }
1024
1025 for (Int_t ipar=0; ipar<5; ipar++){
1026 param1[ipar]= vecXk(ipar,0) ;
1027 for (Int_t jpar=0; jpar<5; jpar++){
1028 covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
1029 }
1030 }
1031}
1032
860b3d93 1033
1034
1035