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