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