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