]> 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
Enable all components in Caibration example macro
[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
16#include "Riostream.h"
17#include "TChain.h"
18#include "TTree.h"
19#include "TH1F.h"
20#include "TH2F.h"
21#include "TList.h"
22#include "TMath.h"
23#include "TCanvas.h"
24#include "TFile.h"
25
26#include "AliTPCseed.h"
27#include "AliESDVertex.h"
28#include "AliESDEvent.h"
29#include "AliESDfriend.h"
30#include "AliESDInputHandler.h"
31
32#include "AliTracker.h"
33#include "AliMagFMaps.h"
34
35#include "AliLog.h"
36
37#include "AliTPCcalibCosmic.h"
38
39#include "TTreeStream.h"
40#include "AliTPCTracklet.h"
41
42ClassImp(AliTPCcalibCosmic)
43
44
45AliTPCcalibCosmic::AliTPCcalibCosmic()
46 :AliTPCcalibBase(),
bca20570 47 fGainMap(0),
9b27d39b 48 fHistNTracks(0),
49 fClusters(0),
50 fModules(0),
51 fHistPt(0),
52 fPtResolution(0),
53 fDeDx(0),
54 fCutMaxD(5), // maximal distance in rfi ditection
55 fCutTheta(0.03), // maximal distan theta
56 fCutMinDir(-0.99) // direction vector products
f7f33dec 57{
58 AliInfo("Defualt Constructor");
bca20570 59 TFile f("/u/miranov/calibKr.root");
60 AliTPCCalPad *gainMap = (AliTPCCalPad *)f.Get("spectrMean");
61 gainMap->Multiply(1/gainMap->GetMedian());
62 fGainMap =gainMap;
f7f33dec 63}
64
65
66AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
67 :AliTPCcalibBase(),
bca20570 68 fGainMap(0),
9b27d39b 69 fHistNTracks(0),
70 fClusters(0),
71 fModules(0),
72 fHistPt(0),
73 fPtResolution(0),
74 fDeDx(0),
75 fCutMaxD(5), // maximal distance in rfi ditection
76 fCutTheta(0.03), // maximal distan theta
77 fCutMinDir(-0.99) // direction vector products
f7f33dec 78{
79 SetName(name);
80 SetTitle(title);
81 AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
82 AliTracker::SetFieldMap(field, kTRUE);
83 fHistNTracks = new TH1F("ntracks","Number of Tracks per Event",501,-0.5,500.5);
84 fClusters = new TH1F("signal","Number of Clusters per track",160,0,160);
85 fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-1200,1200,600,-1000,1000);
86 fHistPt = new TH1F("Pt","Pt distribution",2000,0,50);
87 fPtResolution = new TH1F("PtResolution","Pt resolution",100,-50,50);
88 fDeDx = new TH2F("DeDx","dEdx",500,0.01,20.,500,0.,500);
89 BinLogX(fDeDx);
9b27d39b 90 AliInfo("Non Default Constructor");
bca20570 91 //
92 TFile f("/u/miranov/calibKr.root");
93 AliTPCCalPad *gainMap = (AliTPCCalPad *)f.Get("spectrMean");
94 gainMap->Multiply(1/gainMap->GetMedian());
95 fGainMap =gainMap;
f7f33dec 96}
97
98AliTPCcalibCosmic::~AliTPCcalibCosmic(){
99 //
100 //
101 //
102}
103
104
9b27d39b 105
106
107
f7f33dec 108void AliTPCcalibCosmic::Process(AliESDEvent *event) {
9b27d39b 109 //
110 //
111 //
f7f33dec 112 if (!event) {
113 Printf("ERROR: ESD not available");
114 return;
9b27d39b 115 }
f7f33dec 116 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
117 if (!ESDfriend) {
118 Printf("ERROR: ESDfriend not available");
119 return;
120 }
9b27d39b 121 FindPairs(event);
f7f33dec 122
9b27d39b 123 if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
124 Int_t ntracks=event->GetNumberOfTracks();
125 fHistNTracks->Fill(ntracks);
126 TObjArray tpcSeeds(ntracks);
127 if (ntracks==0) return;
128 //
f7f33dec 129 //track loop
9b27d39b 130 //
131 for (Int_t i=0;i<ntracks;++i) {
f7f33dec 132 AliESDtrack *track = event->GetTrack(i);
9b27d39b 133 fClusters->Fill(track->GetTPCNcls());
f7f33dec 134 AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
135
136 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
137 TObject *calibObject;
138 AliTPCseed *seed = 0;
9b27d39b 139 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
140 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
f7f33dec 141 }
9b27d39b 142 if (seed) tpcSeeds.AddAt(seed,i);
bca20570 143 if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0,0,159,fGainMap));
f7f33dec 144 }
9b27d39b 145 if (ntracks<2) return;
146
147
148 // dE/dx,pt and ACORDE study --> studies which need the pair selection
149 for (Int_t i=0;i<ntracks;++i) {
f7f33dec 150 AliESDtrack *track1 = event->GetTrack(i);
151
152 Double_t d1[3];
153 track1->GetDirection(d1);
154
9b27d39b 155 for (Int_t j=i+1;j<ntracks;++j) {
f7f33dec 156 AliESDtrack *track2 = event->GetTrack(j);
157 Double_t d2[3];
158 track2->GetDirection(d2);
159
160 if (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2] < -0.999) {
161
162 /*___________________________________ Pt resolution ________________________________________*/
163 if (track1->Pt() != 0 && track1->GetTPCNcls() > 80 && track2->GetTPCNcls() > 80) {
164 Double_t res = (track1->Pt() - track2->Pt());
165 res = res/(2*(track1->Pt() + track2->Pt()));
166 fPtResolution->Fill(100*res);
167 }
168
169 /*_______________________________ Propagation to ACORDE ___________________________________*/
170 const Double_t AcordePlane = 850.; //distance of the central Acorde detectors to the beam line at y =0
171 const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
172
173 if (d1[1] > 0 && d2[1] < 0 && track1->GetTPCNcls() > 50) {
174 Double_t r[3];
175 track1->GetXYZ(r);
176 Double_t x,z;
177 z = r[2] + (d1[2]/d1[1])*(AcordePlane - r[1]);
178 x = r[0] + (d1[0]/d1[1])*(AcordePlane - r[1]);
179
180 if (x > roof) {
181 x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
182 z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
183 }
184 if (x < -roof) {
185 x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
186 z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
187 }
188
189 fModules->Fill(z, x);
190 }
191
192 if (d2[1] > 0 && d1[1] < 0 && track2->GetTPCNcls() > 50) {
193 Double_t r[3];
194 track2->GetXYZ(r);
195 Double_t x,z;
196 z = r[2] + (d2[2]/d2[1])*(AcordePlane - r[1]);
197 x = r[0] + (d2[0]/d2[1])*(AcordePlane - r[1]);
198
199 if (x > roof) {
200 x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
201 z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
202 }
203 if (x < -roof) {
204 x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
205 z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
206 }
207
208 fModules->Fill(z, x);
209 }
210
f7f33dec 211 // AliExternalTrackParam * trackOut = new AliExternalTrackParam(*track2->GetOuterParam());
212// AliTracker::PropagateTrackTo(trackOut,850.,105.658,30);
213// delete trackOut;
214
215
216
217
218
219 break;
220 }
221 }
222 }
223
224
225
226
227}
228
9b27d39b 229void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
230 //
231 // Find cosmic pairs
04087794 232 //
233 // Track0 is choosen in upper TPC part
234 // Track1 is choosen in lower TPC part
9b27d39b 235 //
236 if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
237 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
238 Int_t ntracks=event->GetNumberOfTracks();
9b27d39b 239 TObjArray tpcSeeds(ntracks);
240 if (ntracks==0) return;
241 Double_t vtxx[3]={0,0,0};
242 Double_t svtxx[3]={0.000001,0.000001,100.};
243 AliESDVertex vtx(vtxx,svtxx);
244 //
245 //track loop
246 //
247 for (Int_t i=0;i<ntracks;++i) {
248 AliESDtrack *track = event->GetTrack(i);
249 fClusters->Fill(track->GetTPCNcls());
bca20570 250 if (!track->GetInnerParam()) continue;
9b27d39b 251 AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
252
253 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
254 TObject *calibObject;
255 AliTPCseed *seed = 0;
256 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
257 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
258 }
259 if (seed) tpcSeeds.AddAt(seed,i);
bca20570 260 if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0,0,159,fGainMap));
9b27d39b 261 }
262 if (ntracks<2) return;
263 //
264 // Find pairs
265 //
266 for (Int_t i=0;i<ntracks;++i) {
267 AliESDtrack *track0 = event->GetTrack(i);
04087794 268 // track0 - choosen upper part
269 if (!track0) continue;
270 if (!track0->GetOuterParam()) continue;
271 if (track0->GetOuterParam()->GetAlpha()<0) continue;
9b27d39b 272 Double_t d1[3];
273 track0->GetDirection(d1);
04087794 274 for (Int_t j=0;j<ntracks;++j) {
275 if (i==j) continue;
276 AliESDtrack *track1 = event->GetTrack(j);
277 //track 1 lower part
278 if (!track1) continue;
279 if (!track1->GetOuterParam()) continue;
280 if (track1->GetOuterParam()->GetAlpha()>0) continue;
281 //
282 Double_t d2[3];
283 track1->GetDirection(d2);
9b27d39b 284 printf("My stream level=%d\n",fStreamLevel);
285 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
286 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
287 if (! seed0) continue;
288 if (! seed1) continue;
bca20570 289 Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
290 Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0,0,159,fGainMap);
291 //
292 Float_t dedx0I = seed0->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
293 Float_t dedx1I = seed1->CookdEdxNorm(0.05,0.55,0,0,63,fGainMap);
294 //
295 Float_t dedx0O = seed0->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
296 Float_t dedx1O = seed1->CookdEdxNorm(0.05,0.55,0,64,159,fGainMap);
297 //
9b27d39b 298 Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
299 Float_t d0 = track0->GetLinearD(0,0);
300 Float_t d1 = track1->GetLinearD(0,0);
301 //
302 // conservative cuts - convergence to be guarantied
303 // applying before track propagation
304 if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
305 if (dir>fCutMinDir) continue; // direction vector product
306 Float_t bz = AliTracker::GetBz();
307 Float_t dvertex0[2]; //distance to 0,0
308 Float_t dvertex1[2]; //distance to 0,0
309 track0->GetDZ(0,0,0,bz,dvertex0);
310 track1->GetDZ(0,0,0,bz,dvertex1);
311 if (TMath::Abs(dvertex0[1])>250) continue;
312 if (TMath::Abs(dvertex1[1])>250) continue;
313 //
314 //
315 //
316 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
317 AliExternalTrackParam param0(*track0);
318 AliExternalTrackParam param1(*track1);
319 //
320 // Propagate using Magnetic field and correct fo material budget
321 //
322 AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
323 AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
324 //
325 // Propagate rest to the 0,0 DCA - z should be ignored
326 //
327 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
328 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
329 //
330 param0.GetDZ(0,0,0,bz,dvertex0);
331 param1.GetDZ(0,0,0,bz,dvertex1);
332 //
333 Double_t xyz0[3];//,pxyz0[3];
334 Double_t xyz1[3];//,pxyz1[3];
335 param0.GetXYZ(xyz0);
336 param1.GetXYZ(xyz1);
337 Bool_t isPair = IsPair(&param0,&param1);
338 //
339 if (fStreamLevel>0){
340 TTreeSRedirector * cstream = GetDebugStreamer();
341 printf("My stream=%p\n",(void*)cstream);
342 if (cstream) {
343 (*cstream) << "Track0" <<
344 "dir="<<dir<< // direction
345 "OK="<<isPair<< // will be accepted
346 "b0="<<b0<< // propagate status
347 "b1="<<b1<< // propagate status
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
352 "v00="<<dvertex0[0]<< // distance using kalman
353 "v01="<<dvertex0[1]<< //
354 "v10="<<dvertex1[0]<< //
355 "v11="<<dvertex1[1]<< //
356 "d0="<<d0<< // linear distance to 0,0
357 "d1="<<d1<< // linear distance to 0,0
358 //
359 "x00="<<xyz0[0]<<
360 "x01="<<xyz0[1]<<
361 "x02="<<xyz0[2]<<
362 //
363 "x10="<<xyz1[0]<<
364 "x11="<<xyz1[1]<<
365 "x12="<<xyz1[2]<<
366 //
367 "Seed0.=" << track0 << // original seed 0
368 "Seed1.=" << track1 << // original seed 1
bca20570 369 //
370 "dedx0="<<dedx0<< // dedx0 - all
371 "dedx1="<<dedx1<< // dedx1 - all
372 //
373 "dedx0I="<<dedx0I<< // dedx0 - inner
374 "dedx1I="<<dedx1I<< // dedx1 - inner
375 //
376 "dedx0O="<<dedx0O<< // dedx0 - outer
377 "dedx1O="<<dedx1O<< // dedx1 -outer
9b27d39b 378 "\n";
379 }
380 }
381 }
382 }
383}
384
9b27d39b 385
bca20570 386
387
388Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
389
390 TIterator* iter = li->MakeIterator();
391 AliTPCcalibCosmic* cal = 0;
392
393 while ((cal = (AliTPCcalibCosmic*)iter->Next())) {
394 if (!cal->InheritsFrom(AliTPCcalibCosmic::Class())) {
395 Error("Merge","Attempt to add object of class %s to a %s", cal->ClassName(), this->ClassName());
396 return -1;
397 }
398
399 fHistNTracks->Add(cal->fHistNTracks);
400 fClusters->Add(cal->fClusters);
401 fModules->Add(cal->fModules);
402 fHistPt->Add(cal->fHistPt);
403 fPtResolution->Add(cal->fPtResolution);
404 fDeDx->Add(cal->fDeDx);
405
406 }
407
408 return 0;
9b27d39b 409
f7f33dec 410}
411
9b27d39b 412Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
413 //
414 //
415 /*
416 // 0. Same direction - OPOSITE - cutDir +cutT
417 TCut cutDir("cutDir","dir<-0.99")
418 // 1.
419 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
420 //
421 // 2. The same rphi
422 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
423 //
424 //
425 //
426 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
427 // 1/Pt diff cut
428 */
429 const Double_t *p0 = tr0->GetParameter();
430 const Double_t *p1 = tr1->GetParameter();
431 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
432 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
433 Double_t d0[3], d1[3];
434 tr0->GetDirection(d0);
435 tr1->GetDirection(d1);
436 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
437 //
438 return kTRUE;
439}
440
441
f7f33dec 442
443void AliTPCcalibCosmic::BinLogX(TH1 *h) {
444
445 // Method for the correct logarithmic binning of histograms
446
447 TAxis *axis = h->GetXaxis();
448 int bins = axis->GetNbins();
449
450 Double_t from = axis->GetXmin();
451 Double_t to = axis->GetXmax();
452 Double_t *new_bins = new Double_t[bins + 1];
453
454 new_bins[0] = from;
455 Double_t factor = pow(to/from, 1./bins);
456
457 for (int i = 1; i <= bins; i++) {
458 new_bins[i] = factor * new_bins[i-1];
459 }
460 axis->Set(bins, new_bins);
461 delete new_bins;
462
463}
464
7b18d067 465
466
467/*
468
469
bca20570 470
471
7b18d067 472void AliTPCcalibCosmic::dEdxCorrection(){
473 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03");
474 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5");
475 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<0.2&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
bca20570 476 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>110");
7b18d067 477 TCut cutA=cutT+cutD+cutPt+cutN;
478
bca20570 479 TChain * chain = tool.MakeChain("cosmic.txt","Track0",0,1000000);
480 chain->Lookup();
7b18d067 481
482 .x ~/rootlogon.C
483 gSystem->Load("libSTAT.so");
484
485 Double_t chi2=0;
486 Int_t npoints=0;
487 TVectorD fitParam;
488 TMatrixD covMatrix;
489
490 chain->Draw("Tr0.fP[4]+Tr1.fP[4]","OK"+cutA);
491
492 TString strFit;
493 strFit+="(Tr0.fP[1]/250)++";
494 strFit+="(Tr0.fP[1]/250)^2++";
495 strFit+="(Tr0.fP[3])++";
496 strFit+="(Tr0.fP[3])^2++";
497
498 TString * ptParam = TStatToolkit::FitPlane(chain,"Tr0.fP[4]+Tr1.fP[4]", strFit.Data(),cutA, chi2,npoints,fitParam,covMatrix)
499
500strFit+="(Tr0.fP[1]/250)++";
501strFit+="(Tr0.fP[1]/250)^2++";
502strFit+="(Tr0.fP[3])++";
503strFit+="(Tr0.fP[3])^2++";
504strFit+="(Tr0.fP[1]/250)^2*Tr0.fP[3]++";
505strFit+="(Tr0.fP[1]/250)^2*Tr0.fP[3]^2++";
506//
507
508strFit+="sign(Tr0.fP[1])++"
509strFit+="sign(Tr0.fP[1])*(1-abs(Tr0.fP[1]/250))"
510
511TString * thetaParam = TStatToolkit::FitPlane(chain,"Tr0.fP[3]+Tr1.fP[3]", strFit.Data(),cutA, chi2,npoints,fitParam,covMatrix)
bca20570 512
513
514
515(-0.009263+(Tr0.fP[1]/250)*(-0.009693)+(Tr0.fP[1]/250)^2*(0.009062)+(Tr0.fP[3])*(0.002256)+(Tr0.fP[3])^2*(-0.052775)+(Tr0.fP[1]/250)^2*Tr0.fP[3]*(-0.020627)+(Tr0.fP[1]/250)^2*Tr0.fP[3]^2*(0.049030))
516
7b18d067 517*/
518
519