- Prepare covariance matrix manipulation in the tracklet for resolution
[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
36
37*/
38
39
40
f7f33dec 41#include "Riostream.h"
42#include "TChain.h"
43#include "TTree.h"
44#include "TH1F.h"
45#include "TH2F.h"
46#include "TList.h"
47#include "TMath.h"
48#include "TCanvas.h"
49#include "TFile.h"
54b76c13 50#include "TF1.h"
f7f33dec 51
54b76c13 52#include "AliTPCclusterMI.h"
f7f33dec 53#include "AliTPCseed.h"
54#include "AliESDVertex.h"
55#include "AliESDEvent.h"
56#include "AliESDfriend.h"
57#include "AliESDInputHandler.h"
54b76c13 58#include "AliAnalysisManager.h"
15e48021 59#include "AliExternalComparison.h"
f7f33dec 60
61#include "AliTracker.h"
f7a1cc68 62#include "AliMagF.h"
54b76c13 63#include "AliTPCCalROC.h"
f7f33dec 64
65#include "AliLog.h"
66
67#include "AliTPCcalibCosmic.h"
f7f33dec 68#include "TTreeStream.h"
69#include "AliTPCTracklet.h"
15e48021 70#include "AliESDcosmic.h"
71
f7f33dec 72
73ClassImp(AliTPCcalibCosmic)
74
75
76AliTPCcalibCosmic::AliTPCcalibCosmic()
77 :AliTPCcalibBase(),
15e48021 78 fComp(0),
9b27d39b 79 fHistNTracks(0),
80 fClusters(0),
81 fModules(0),
82 fHistPt(0),
9b27d39b 83 fDeDx(0),
54b76c13 84 fDeDxMIP(0),
85 fMIPvalue(1),
9b27d39b 86 fCutMaxD(5), // maximal distance in rfi ditection
a6dc0cf6 87 fCutMaxDz(40), // maximal distance in z ditection
9b27d39b 88 fCutTheta(0.03), // maximal distan theta
89 fCutMinDir(-0.99) // direction vector products
f7f33dec 90{
54b76c13 91 AliInfo("Default Constructor");
f7f33dec 92}
93
94
95AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
96 :AliTPCcalibBase(),
15e48021 97 fComp(0),
9b27d39b 98 fHistNTracks(0),
99 fClusters(0),
100 fModules(0),
101 fHistPt(0),
9b27d39b 102 fDeDx(0),
54b76c13 103 fDeDxMIP(0),
104 fMIPvalue(1),
a6dc0cf6 105 fCutMaxD(5), // maximal distance in rfi ditection
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{
110 SetName(name);
111 SetTitle(title);
54b76c13 112
113 fHistNTracks = new TH1F("ntracks","Number of Tracks per Event; number of tracks per event; number of tracks",501,-0.5,500.5);
114 fClusters = new TH1F("signal","Number of Clusters per track; number of clusters per track n_{cl}; counts",160,0,160);
115 fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-650,650,600,-700,700);
116 fHistPt = new TH1F("Pt","Pt distribution; p_{T} (GeV); counts",2000,0,50);
117 fDeDx = new TH2F("DeDx","dEdx; momentum p (GeV); TPC signal (a.u.)",500,0.01,100.,500,2.,1000);
f7f33dec 118 BinLogX(fDeDx);
54b76c13 119 fDeDxMIP = new TH1F("DeDxMIP","MIP region; TPC signal (a.u.);counts ",500,2.,1000);
120
9b27d39b 121 AliInfo("Non Default Constructor");
bca20570 122 //
f7f33dec 123}
124
125AliTPCcalibCosmic::~AliTPCcalibCosmic(){
126 //
127 //
128 //
15e48021 129 delete fComp;
f7f33dec 130}
131
132
9b27d39b 133
134
135
f7f33dec 136void AliTPCcalibCosmic::Process(AliESDEvent *event) {
9b27d39b 137 //
138 //
139 //
f7f33dec 140 if (!event) {
141 Printf("ERROR: ESD not available");
142 return;
9b27d39b 143 }
f7f33dec 144 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
145 if (!ESDfriend) {
146 Printf("ERROR: ESDfriend not available");
147 return;
148 }
2acad464 149
108953e9 150
54b76c13 151 FindPairs(event); // nearly everything takes place in find pairs...
152
2acad464 153 if (GetDebugLevel()>20) printf("Hallo world: Im here and processing an event\n");
9b27d39b 154 Int_t ntracks=event->GetNumberOfTracks();
155 fHistNTracks->Fill(ntracks);
9b27d39b 156 if (ntracks==0) return;
15e48021 157 AliESDcosmic cosmicESD;
158 TTreeSRedirector * cstream = GetDebugStreamer();
159 cosmicESD.SetDebugStreamer(cstream);
160 cosmicESD.ProcessEvent(event);
161 if (cstream) cosmicESD.DumpToTree();
162
163
54b76c13 164}
9b27d39b 165
f7f33dec 166
167
54b76c13 168void AliTPCcalibCosmic::Analyze() {
169
170 fMIPvalue = CalculateMIPvalue(fDeDxMIP);
171
172 return;
173
174}
175
f7f33dec 176
f7f33dec 177
9b27d39b 178void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
179 //
180 // Find cosmic pairs
04087794 181 //
182 // Track0 is choosen in upper TPC part
183 // Track1 is choosen in lower TPC part
9b27d39b 184 //
2acad464 185 if (GetDebugLevel()>20) printf("Hallo world: Im here\n");
9b27d39b 186 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
187 Int_t ntracks=event->GetNumberOfTracks();
9b27d39b 188 TObjArray tpcSeeds(ntracks);
189 if (ntracks==0) return;
190 Double_t vtxx[3]={0,0,0};
191 Double_t svtxx[3]={0.000001,0.000001,100.};
192 AliESDVertex vtx(vtxx,svtxx);
193 //
194 //track loop
195 //
54b76c13 196 for (Int_t i=0;i<ntracks;++i) {
197 AliESDtrack *track = event->GetTrack(i);
198 fClusters->Fill(track->GetTPCNcls());
199
200 const AliExternalTrackParam * trackIn = track->GetInnerParam();
201 const AliExternalTrackParam * trackOut = track->GetOuterParam();
202 if (!trackIn) continue;
203 if (!trackOut) continue;
860b3d93 204 if (ntracks>4 && TMath::Abs(trackIn->GetTgl())<0.0015) continue; // filter laser
205
206
9b27d39b 207 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
208 TObject *calibObject;
209 AliTPCseed *seed = 0;
210 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
211 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
212 }
213 if (seed) tpcSeeds.AddAt(seed,i);
54b76c13 214
215 Double_t meanP = 0.5*(trackIn->GetP() + trackOut->GetP());
216 if (seed && track->GetTPCNcls() > 80 + 60/(1+TMath::Exp(-meanP+5))) {
15e48021 217 fDeDx->Fill(meanP, seed->CookdEdxNorm(0.0,0.45,0,0,159));
54b76c13 218 //
15e48021 219 if (meanP > 0.4 && meanP < 0.45) fDeDxMIP->Fill(seed->CookdEdxNorm(0.0,0.45,0,0,159));
54b76c13 220 //
15e48021 221 if (GetDebugLevel()>0&&meanP>0.2&&seed->CookdEdxNorm(0.0,0.45,0,0,159)>300) {
54b76c13 222 TFile *curfile = AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile();
223 if (curfile) printf(">>> p+ in file: %s \t event: %i \t Number of ESD tracks: %i \n", curfile->GetName(), (int)event->GetEventNumberInFile(), (int)ntracks);
224 if (track->GetOuterParam()->GetAlpha()<0) cout << " Polartiy: " << track->GetSign() << endl;
225 }
226
227 }
228
9b27d39b 229 }
54b76c13 230
9b27d39b 231 if (ntracks<2) return;
232 //
233 // Find pairs
234 //
235 for (Int_t i=0;i<ntracks;++i) {
236 AliESDtrack *track0 = event->GetTrack(i);
04087794 237 // track0 - choosen upper part
238 if (!track0) continue;
239 if (!track0->GetOuterParam()) continue;
240 if (track0->GetOuterParam()->GetAlpha()<0) continue;
e9362f9d 241 Double_t dir0[3];
242 track0->GetDirection(dir0);
04087794 243 for (Int_t j=0;j<ntracks;++j) {
244 if (i==j) continue;
245 AliESDtrack *track1 = event->GetTrack(j);
246 //track 1 lower part
247 if (!track1) continue;
248 if (!track1->GetOuterParam()) continue;
249 if (track1->GetOuterParam()->GetAlpha()>0) continue;
250 //
e9362f9d 251 Double_t dir1[3];
252 track1->GetDirection(dir1);
54b76c13 253
9b27d39b 254 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
255 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
256 if (! seed0) continue;
257 if (! seed1) continue;
15e48021 258 Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159);
259 Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159);
bca20570 260 //
15e48021 261 Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63);
262 Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63);
bca20570 263 //
15e48021 264 Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159);
265 Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159);
bca20570 266 //
e9362f9d 267 Float_t dir = (dir0[0]*dir1[0] + dir0[1]*dir1[1] + dir0[2]*dir1[2]);
9b27d39b 268 Float_t d0 = track0->GetLinearD(0,0);
269 Float_t d1 = track1->GetLinearD(0,0);
270 //
271 // conservative cuts - convergence to be guarantied
272 // applying before track propagation
273 if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
274 if (dir>fCutMinDir) continue; // direction vector product
275 Float_t bz = AliTracker::GetBz();
276 Float_t dvertex0[2]; //distance to 0,0
277 Float_t dvertex1[2]; //distance to 0,0
278 track0->GetDZ(0,0,0,bz,dvertex0);
279 track1->GetDZ(0,0,0,bz,dvertex1);
280 if (TMath::Abs(dvertex0[1])>250) continue;
281 if (TMath::Abs(dvertex1[1])>250) continue;
282 //
283 //
284 //
285 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
286 AliExternalTrackParam param0(*track0);
287 AliExternalTrackParam param1(*track1);
288 //
289 // Propagate using Magnetic field and correct fo material budget
290 //
291 AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
292 AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
293 //
294 // Propagate rest to the 0,0 DCA - z should be ignored
295 //
296 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
297 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
298 //
299 param0.GetDZ(0,0,0,bz,dvertex0);
300 param1.GetDZ(0,0,0,bz,dvertex1);
a6dc0cf6 301 if (TMath::Abs(param0.GetZ()-param1.GetZ())>fCutMaxDz) continue;
9b27d39b 302 //
303 Double_t xyz0[3];//,pxyz0[3];
304 Double_t xyz1[3];//,pxyz1[3];
305 param0.GetXYZ(xyz0);
306 param1.GetXYZ(xyz1);
307 Bool_t isPair = IsPair(&param0,&param1);
308 //
54b76c13 309 if (isPair) FillAcordeHist(track0);
15e48021 310 if (fComp){
311 Bool_t acceptComp = fComp->AcceptPair(&param0,&param1);
312 if (acceptComp) fComp->Process(&param0,&param1);
313 }
314 //
315 // combined track params
316 //
317 AliExternalTrackParam *par0U=MakeCombinedTrack(&param0,&param1);
318 AliExternalTrackParam *par1U=MakeCombinedTrack(&param1,&param0);
319
320
54b76c13 321 //
9b27d39b 322 if (fStreamLevel>0){
323 TTreeSRedirector * cstream = GetDebugStreamer();
108953e9 324 //printf("My stream=%p\n",(void*)cstream);
e9362f9d 325 AliExternalTrackParam *ip0 = (AliExternalTrackParam *)track0->GetInnerParam();
326 AliExternalTrackParam *ip1 = (AliExternalTrackParam *)track1->GetInnerParam();
327 AliExternalTrackParam *op0 = (AliExternalTrackParam *)track0->GetOuterParam();
328 AliExternalTrackParam *op1 = (AliExternalTrackParam *)track1->GetOuterParam();
329 Bool_t isCrossI = ip0->GetZ()*ip1->GetZ()<0;
330 Bool_t isCrossO = op0->GetZ()*op1->GetZ()<0;
331 Double_t alpha0 = TMath::ATan2(dir0[1],dir0[0]);
332 Double_t alpha1 = TMath::ATan2(dir1[1],dir1[0]);
9b27d39b 333 if (cstream) {
334 (*cstream) << "Track0" <<
108953e9 335 "run="<<fRun<< // run number
336 "event="<<fEvent<< // event number
337 "time="<<fTime<< // time stamp of event
338 "trigger="<<fTrigger<< // trigger
15e48021 339 "triggerClass="<<&fTriggerClass<< // trigger
108953e9 340 "mag="<<fMagF<< // magnetic field
9b27d39b 341 "dir="<<dir<< // direction
54b76c13 342 "OK="<<isPair<< // will be accepted
9b27d39b 343 "b0="<<b0<< // propagate status
344 "b1="<<b1<< // propagate status
e9362f9d 345 "crossI="<<isCrossI<< // cross inner
346 "crossO="<<isCrossO<< // cross outer
347 //
9b27d39b 348 "Orig0.=" << track0 << // original track 0
349 "Orig1.=" << track1 << // original track 1
350 "Tr0.="<<&param0<< // track propagated to the DCA 0,0
351 "Tr1.="<<&param1<< // track propagated to the DCA 0,0
e9362f9d 352 "Ip0.="<<ip0<< // inner param - upper
353 "Ip1.="<<ip1<< // inner param - lower
354 "Op0.="<<op0<< // outer param - upper
355 "Op1.="<<op1<< // outer param - lower
15e48021 356 "Up0.="<<par0U<< // combined track 0
357 "Up1.="<<par1U<< // combined track 1
e9362f9d 358 //
9b27d39b 359 "v00="<<dvertex0[0]<< // distance using kalman
360 "v01="<<dvertex0[1]<< //
361 "v10="<<dvertex1[0]<< //
362 "v11="<<dvertex1[1]<< //
363 "d0="<<d0<< // linear distance to 0,0
364 "d1="<<d1<< // linear distance to 0,0
365 //
e9362f9d 366 //
367 //
368 "x00="<<xyz0[0]<< // global position close to vertex
9b27d39b 369 "x01="<<xyz0[1]<<
370 "x02="<<xyz0[2]<<
371 //
e9362f9d 372 "x10="<<xyz1[0]<< // global position close to vertex
9b27d39b 373 "x11="<<xyz1[1]<<
374 "x12="<<xyz1[2]<<
375 //
e9362f9d 376 "alpha0="<<alpha0<<
377 "alpha1="<<alpha1<<
378 "dir00="<<dir0[0]<< // direction upper
379 "dir01="<<dir0[1]<<
380 "dir02="<<dir0[2]<<
381 //
382 "dir10="<<dir1[0]<< // direction lower
383 "dir11="<<dir1[1]<<
384 "dir12="<<dir1[2]<<
385 //
386 //
54b76c13 387 "Seed0.=" << seed0 << // original seed 0
388 "Seed1.=" << seed1 << // original seed 1
bca20570 389 //
390 "dedx0="<<dedx0<< // dedx0 - all
391 "dedx1="<<dedx1<< // dedx1 - all
392 //
54b76c13 393 "dedx0I="<<dedx0I<< // dedx0 - inner ROC
394 "dedx1I="<<dedx1I<< // dedx1 - inner ROC
bca20570 395 //
54b76c13 396 "dedx0O="<<dedx0O<< // dedx0 - outer ROC
397 "dedx1O="<<dedx1O<< // dedx1 - outer ROC
9b27d39b 398 "\n";
399 }
15e48021 400 }
401 delete par0U;
402 delete par1U;
9b27d39b 403 }
404 }
405}
406
9b27d39b 407
bca20570 408
409
54b76c13 410void AliTPCcalibCosmic::FillAcordeHist(AliESDtrack *upperTrack) {
411
412 // Pt cut to select straight tracks which can be easily propagated to ACORDE which is outside the magnetic field
413 if (upperTrack->Pt() < 10 || upperTrack->GetTPCNcls() < 80) return;
414
415 const Double_t AcordePlane = 850.; // distance of the central Acorde detectors to the beam line at y =0
416 const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
417
418 Double_t r[3];
419 upperTrack->GetXYZ(r);
420 Double_t d[3];
421 upperTrack->GetDirection(d);
422 Double_t x,z;
423 z = r[2] + (d[2]/d[1])*(AcordePlane - r[1]);
424 x = r[0] + (d[0]/d[1])*(AcordePlane - r[1]);
425
426 if (x > roof) {
427 x = r[0] + (d[0]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
428 z = r[2] + (d[2]/(d[0]+d[1]))*(AcordePlane+roof-r[0]-r[1]);
429 }
430 if (x < -roof) {
431 x = r[0] + (d[0]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
432 z = r[2] + (d[2]/(d[1]-d[0]))*(AcordePlane+roof+r[0]-r[1]);
433 }
434
435 fModules->Fill(z, x);
436
437}
438
439
440
bca20570 441Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
442
443 TIterator* iter = li->MakeIterator();
444 AliTPCcalibCosmic* cal = 0;
445
446 while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
447 if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
860b3d93 448 //Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
bca20570 449 return -1;
450 }
451
54b76c13 452 fHistNTracks->Add(cal->GetHistNTracks());
453 fClusters->Add(cal-> GetHistClusters());
454 fModules->Add(cal->GetHistAcorde());
455 fHistPt->Add(cal->GetHistPt());
456 fDeDx->Add(cal->GetHistDeDx());
457 fDeDxMIP->Add(cal->GetHistMIP());
bca20570 458
459 }
15e48021 460 //if (fComp && cal->fComp) fComp->Add(cal->fComp);
bca20570 461 return 0;
9b27d39b 462
f7f33dec 463}
464
54b76c13 465
9b27d39b 466Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
467 //
468 //
469 /*
470 // 0. Same direction - OPOSITE - cutDir +cutT
471 TCut cutDir("cutDir","dir<-0.99")
472 // 1.
473 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
474 //
475 // 2. The same rphi
476 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
477 //
478 //
479 //
480 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
481 // 1/Pt diff cut
482 */
483 const Double_t *p0 = tr0->GetParameter();
484 const Double_t *p1 = tr1->GetParameter();
485 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
a6dc0cf6 486 if (TMath::Abs(p0[1]-p1[1])>fCutMaxDz) return kFALSE;
9b27d39b 487 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
a6dc0cf6 488
9b27d39b 489 Double_t d0[3], d1[3];
490 tr0->GetDirection(d0);
491 tr1->GetDirection(d1);
492 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
493 //
494 return kTRUE;
495}
54b76c13 496
497
498
499Double_t AliTPCcalibCosmic::CalculateMIPvalue(TH1F * hist) {
500
501 TF1 * funcDoubleGaus = new TF1("funcDoubleGaus", "gaus(0)+gaus(3)",0,1000);
502 funcDoubleGaus->SetParameters(hist->GetEntries()*0.75,hist->GetMean()/1.3,hist->GetMean()*0.10,
503 hist->GetEntries()*0.25,hist->GetMean()*1.3,hist->GetMean()*0.10);
504 hist->Fit(funcDoubleGaus);
505 Double_t MIPvalue = TMath::Min(funcDoubleGaus->GetParameter(1),funcDoubleGaus->GetParameter(4));
506
507 delete funcDoubleGaus;
508
509 return MIPvalue;
510
511}
9b27d39b 512
513
f7f33dec 514
54b76c13 515
57dc06f2 516void AliTPCcalibCosmic::CalculateBetheParams(TH2F */*hist*/, Double_t * /*initialParam*/) {
517 //
518 // Not implemented yet
519 //
54b76c13 520 return;
521
522}
523
524
f7f33dec 525void AliTPCcalibCosmic::BinLogX(TH1 *h) {
526
527 // Method for the correct logarithmic binning of histograms
528
529 TAxis *axis = h->GetXaxis();
530 int bins = axis->GetNbins();
531
532 Double_t from = axis->GetXmin();
533 Double_t to = axis->GetXmax();
534 Double_t *new_bins = new Double_t[bins + 1];
535
536 new_bins[0] = from;
537 Double_t factor = pow(to/from, 1./bins);
538
539 for (int i = 1; i <= bins; i++) {
540 new_bins[i] = factor * new_bins[i-1];
541 }
542 axis->Set(bins, new_bins);
543 delete new_bins;
544
545}
546
7b18d067 547
860b3d93 548AliExternalTrackParam *AliTPCcalibCosmic::MakeTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
549 //
550 //
551 //
552 AliExternalTrackParam *par1R= new AliExternalTrackParam(*track1);
553 par1R->Rotate(track0->GetAlpha());
15e48021 554 par1R->PropagateTo(track0->GetX(),AliTracker::GetBz());
860b3d93 555 //
556 //
557 Double_t * param = (Double_t*)par1R->GetParameter();
558 Double_t * covar = (Double_t*)par1R->GetCovariance();
15e48021 559
860b3d93 560 param[0]*=1; //OK
561 param[1]*=1; //OK
562 param[2]*=1; //?
563 param[3]*=-1; //OK
564 param[4]*=-1; //OK
565 //
566 covar[6] *=-1.; covar[7] *=-1.; covar[8] *=-1.;
567 //covar[10]*=-1.; covar[11]*=-1.; covar[12]*=-1.;
568 covar[13]*=-1.;
860b3d93 569 return par1R;
570}
571
15e48021 572AliExternalTrackParam *AliTPCcalibCosmic::MakeCombinedTrack(const AliExternalTrackParam *track0, const AliExternalTrackParam *track1){
573 //
574 // Make combined track
575 //
576 //
577 AliExternalTrackParam * par1T = MakeTrack(track0,track1);
578 AliExternalTrackParam * par0U = new AliExternalTrackParam(*track0);
579 //
580 UpdateTrack(*par0U,*par1T);
581 delete par1T;
582 return par0U;
583}
584
585
860b3d93 586void AliTPCcalibCosmic::UpdateTrack(AliExternalTrackParam &track1, const AliExternalTrackParam &track2){
587 //
588 // Update track 1 with track 2
589 //
590 //
591 //
592 TMatrixD vecXk(5,1); // X vector
593 TMatrixD covXk(5,5); // X covariance
594 TMatrixD matHk(5,5); // vector to mesurement
595 TMatrixD measR(5,5); // measurement error
596 TMatrixD vecZk(5,1); // measurement
597 //
598 TMatrixD vecYk(5,1); // Innovation or measurement residual
599 TMatrixD matHkT(5,5);
600 TMatrixD matSk(5,5); // Innovation (or residual) covariance
601 TMatrixD matKk(5,5); // Optimal Kalman gain
602 TMatrixD mat1(5,5); // update covariance matrix
603 TMatrixD covXk2(5,5); //
604 TMatrixD covOut(5,5);
605 //
606 Double_t *param1=(Double_t*) track1.GetParameter();
607 Double_t *covar1=(Double_t*) track1.GetCovariance();
608 Double_t *param2=(Double_t*) track2.GetParameter();
609 Double_t *covar2=(Double_t*) track2.GetCovariance();
610 //
611 // copy data to the matrix
612 for (Int_t ipar=0; ipar<5; ipar++){
860b3d93 613 for (Int_t jpar=0; jpar<5; jpar++){
614 covXk(ipar,jpar) = covar1[track1.GetIndex(ipar, jpar)];
615 measR(ipar,jpar) = covar2[track2.GetIndex(ipar, jpar)];
15e48021 616 matHk(ipar,jpar)=0;
617 mat1(ipar,jpar)=0;
860b3d93 618 }
15e48021 619 vecXk(ipar,0) = param1[ipar];
620 vecZk(ipar,0) = param2[ipar];
621 matHk(ipar,ipar)=1;
622 mat1(ipar,ipar)=0;
860b3d93 623 }
624 //
625 //
626 //
627 //
860b3d93 628 //
629 vecYk = vecZk-matHk*vecXk; // Innovation or measurement residual
630 matHkT=matHk.T(); matHk.T();
631 matSk = (matHk*(covXk*matHkT))+measR; // Innovation (or residual) covariance
632 matSk.Invert();
633 matKk = (covXk*matHkT)*matSk; // Optimal Kalman gain
634 vecXk += matKk*vecYk; // updated vector
860b3d93 635 covXk2 = (mat1-(matKk*matHk));
636 covOut = covXk2*covXk;
637 //
638 //
639 //
640 // copy from matrix to parameters
641 if (0) {
642 vecXk.Print();
643 vecZk.Print();
644 //
645 measR.Print();
646 covXk.Print();
647 covOut.Print();
648 //
649 track1.Print();
650 track2.Print();
651 }
652
653 for (Int_t ipar=0; ipar<5; ipar++){
654 param1[ipar]= vecXk(ipar,0) ;
655 for (Int_t jpar=0; jpar<5; jpar++){
656 covar1[track1.GetIndex(ipar, jpar)]=covOut(ipar,jpar);
657 }
658 }
659}
660
661void AliTPCcalibCosmic::ProcessTree(TTree * chainTracklet, AliExternalComparison *comp){
662 //
663 // Process the debug streamer tree
664 // Possible to modify selection criteria
665 //
666 TTreeSRedirector * cstream = new TTreeSRedirector("cosmicdump.root");
667 //AliTPCcalibCosmic *cosmic = this;
668 //
669 AliExternalTrackParam * tr0 = 0;
670 AliExternalTrackParam * tr1 = 0;
671 Int_t npoints =0;
672 {
673 Int_t entries=chainTracklet->GetEntries();
674 for (Int_t i=0; i< entries; i++){
675 chainTracklet->GetBranch("Tr0.")->SetAddress(&tr0);
676 chainTracklet->GetBranch("Tr1.")->SetAddress(&tr1);
677 chainTracklet->GetEntry(i);
678 if (!tr0) continue;
679 if (!tr1) continue;
680 if (tr0->GetY()==0) continue;
681 if (tr1->GetY()==0) continue;
682 // make a local copy
683 AliExternalTrackParam par0(*tr0);
684 AliExternalTrackParam par1(*tr1);
685 AliExternalTrackParam par1R(*tr1);
686 par1R.Rotate(par1.GetAlpha()+TMath::Pi());
687 AliExternalTrackParam *par1T = MakeTrack(tr0,tr1);
688 if (0) {
689 printf("%d\t%d\t\n",i, npoints);
690 par1R.Print();
691 par1T->Print();
692 }
693 AliExternalTrackParam par0U=par0;
694 AliExternalTrackParam par1U=*par1T;
695 //
696 UpdateTrack(par0U,*par1T);
697 UpdateTrack(par1U,par0);
698 //
699 //
700 if (i%100==0) printf("%d\t%d\tt\n",i, npoints);
701 Bool_t accept = comp->AcceptPair(&par0,par1T);
702
703 if (1||fStreamLevel>0){
704 (*cstream)<<"Tracklet"<<
705 "accept="<<accept<<
706 "tr0.="<<&par0<< //original track up
707 "tr1.="<<&par1<< //original track down
708 "tr1R.="<<&par1R<< //track1 rotated to 0 frame
709 "tr1T.="<<par1T<< //track1 transformed to the track 0 frame
710 //
711 "tr0U.="<<&par0U<< //track 0 updated with track 1
712 "tr1U.="<<&par1U<< //track 1 updated with track 0
713 "\n";
714 }
715 //
716 if (accept) {
717 npoints++;
718 if (comp) comp->Process(&par0,par1T);
719 }
720 delete par1T;
721 }
722 }
723 delete cstream;
724}
725
726
727
728