]> 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
Adding possibility to get run number (Marian)
[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(),
9b27d39b 47 fHistNTracks(0),
48 fClusters(0),
49 fModules(0),
50 fHistPt(0),
51 fPtResolution(0),
52 fDeDx(0),
53 fCutMaxD(5), // maximal distance in rfi ditection
54 fCutTheta(0.03), // maximal distan theta
55 fCutMinDir(-0.99) // direction vector products
f7f33dec 56{
57 AliInfo("Defualt Constructor");
58}
59
60
61AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
62 :AliTPCcalibBase(),
9b27d39b 63 fHistNTracks(0),
64 fClusters(0),
65 fModules(0),
66 fHistPt(0),
67 fPtResolution(0),
68 fDeDx(0),
69 fCutMaxD(5), // maximal distance in rfi ditection
70 fCutTheta(0.03), // maximal distan theta
71 fCutMinDir(-0.99) // direction vector products
f7f33dec 72{
73 SetName(name);
74 SetTitle(title);
75 AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
76 AliTracker::SetFieldMap(field, kTRUE);
77 fHistNTracks = new TH1F("ntracks","Number of Tracks per Event",501,-0.5,500.5);
78 fClusters = new TH1F("signal","Number of Clusters per track",160,0,160);
79 fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-1200,1200,600,-1000,1000);
80 fHistPt = new TH1F("Pt","Pt distribution",2000,0,50);
81 fPtResolution = new TH1F("PtResolution","Pt resolution",100,-50,50);
82 fDeDx = new TH2F("DeDx","dEdx",500,0.01,20.,500,0.,500);
83 BinLogX(fDeDx);
9b27d39b 84 AliInfo("Non Default Constructor");
f7f33dec 85}
86
87AliTPCcalibCosmic::~AliTPCcalibCosmic(){
88 //
89 //
90 //
91}
92
93
9b27d39b 94
95
96
f7f33dec 97void AliTPCcalibCosmic::Process(AliESDEvent *event) {
9b27d39b 98 //
99 //
100 //
f7f33dec 101 if (!event) {
102 Printf("ERROR: ESD not available");
103 return;
9b27d39b 104 }
f7f33dec 105 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
106 if (!ESDfriend) {
107 Printf("ERROR: ESDfriend not available");
108 return;
109 }
9b27d39b 110 FindPairs(event);
f7f33dec 111
9b27d39b 112 if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
113 Int_t ntracks=event->GetNumberOfTracks();
114 fHistNTracks->Fill(ntracks);
115 TObjArray tpcSeeds(ntracks);
116 if (ntracks==0) return;
117 //
f7f33dec 118 //track loop
9b27d39b 119 //
120 for (Int_t i=0;i<ntracks;++i) {
f7f33dec 121 AliESDtrack *track = event->GetTrack(i);
9b27d39b 122 fClusters->Fill(track->GetTPCNcls());
f7f33dec 123 AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
124
125 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
126 TObject *calibObject;
127 AliTPCseed *seed = 0;
9b27d39b 128 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
129 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
f7f33dec 130 }
9b27d39b 131 if (seed) tpcSeeds.AddAt(seed,i);
132 if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
f7f33dec 133 }
9b27d39b 134 if (ntracks<2) return;
135
136
137 // dE/dx,pt and ACORDE study --> studies which need the pair selection
138 for (Int_t i=0;i<ntracks;++i) {
f7f33dec 139 AliESDtrack *track1 = event->GetTrack(i);
140
141 Double_t d1[3];
142 track1->GetDirection(d1);
143
9b27d39b 144 for (Int_t j=i+1;j<ntracks;++j) {
f7f33dec 145 AliESDtrack *track2 = event->GetTrack(j);
146 Double_t d2[3];
147 track2->GetDirection(d2);
148
149 if (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2] < -0.999) {
150
151 /*___________________________________ Pt resolution ________________________________________*/
152 if (track1->Pt() != 0 && track1->GetTPCNcls() > 80 && track2->GetTPCNcls() > 80) {
153 Double_t res = (track1->Pt() - track2->Pt());
154 res = res/(2*(track1->Pt() + track2->Pt()));
155 fPtResolution->Fill(100*res);
156 }
157
158 /*_______________________________ Propagation to ACORDE ___________________________________*/
159 const Double_t AcordePlane = 850.; //distance of the central Acorde detectors to the beam line at y =0
160 const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
161
162 if (d1[1] > 0 && d2[1] < 0 && track1->GetTPCNcls() > 50) {
163 Double_t r[3];
164 track1->GetXYZ(r);
165 Double_t x,z;
166 z = r[2] + (d1[2]/d1[1])*(AcordePlane - r[1]);
167 x = r[0] + (d1[0]/d1[1])*(AcordePlane - r[1]);
168
169 if (x > roof) {
170 x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
171 z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
172 }
173 if (x < -roof) {
174 x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
175 z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
176 }
177
178 fModules->Fill(z, x);
179 }
180
181 if (d2[1] > 0 && d1[1] < 0 && track2->GetTPCNcls() > 50) {
182 Double_t r[3];
183 track2->GetXYZ(r);
184 Double_t x,z;
185 z = r[2] + (d2[2]/d2[1])*(AcordePlane - r[1]);
186 x = r[0] + (d2[0]/d2[1])*(AcordePlane - r[1]);
187
188 if (x > roof) {
189 x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
190 z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
191 }
192 if (x < -roof) {
193 x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
194 z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
195 }
196
197 fModules->Fill(z, x);
198 }
199
f7f33dec 200 // AliExternalTrackParam * trackOut = new AliExternalTrackParam(*track2->GetOuterParam());
201// AliTracker::PropagateTrackTo(trackOut,850.,105.658,30);
202// delete trackOut;
203
204
205
206
207
208 break;
209 }
210 }
211 }
212
213
214
215
216}
217
9b27d39b 218void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
219 //
220 // Find cosmic pairs
04087794 221 //
222 // Track0 is choosen in upper TPC part
223 // Track1 is choosen in lower TPC part
9b27d39b 224 //
225 if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
226 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
227 Int_t ntracks=event->GetNumberOfTracks();
9b27d39b 228 TObjArray tpcSeeds(ntracks);
229 if (ntracks==0) return;
230 Double_t vtxx[3]={0,0,0};
231 Double_t svtxx[3]={0.000001,0.000001,100.};
232 AliESDVertex vtx(vtxx,svtxx);
233 //
234 //track loop
235 //
236 for (Int_t i=0;i<ntracks;++i) {
237 AliESDtrack *track = event->GetTrack(i);
238 fClusters->Fill(track->GetTPCNcls());
239 AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
240
241 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
242 TObject *calibObject;
243 AliTPCseed *seed = 0;
244 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
245 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
246 }
247 if (seed) tpcSeeds.AddAt(seed,i);
248 if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
249 }
250 if (ntracks<2) return;
251 //
252 // Find pairs
253 //
254 for (Int_t i=0;i<ntracks;++i) {
255 AliESDtrack *track0 = event->GetTrack(i);
04087794 256 // track0 - choosen upper part
257 if (!track0) continue;
258 if (!track0->GetOuterParam()) continue;
259 if (track0->GetOuterParam()->GetAlpha()<0) continue;
9b27d39b 260 Double_t d1[3];
261 track0->GetDirection(d1);
04087794 262 for (Int_t j=0;j<ntracks;++j) {
263 if (i==j) continue;
264 AliESDtrack *track1 = event->GetTrack(j);
265 //track 1 lower part
266 if (!track1) continue;
267 if (!track1->GetOuterParam()) continue;
268 if (track1->GetOuterParam()->GetAlpha()>0) continue;
269 //
270 Double_t d2[3];
271 track1->GetDirection(d2);
9b27d39b 272 printf("My stream level=%d\n",fStreamLevel);
273 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
274 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
275 if (! seed0) continue;
276 if (! seed1) continue;
277 Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0);
278 Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0);
279 Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
280 Float_t d0 = track0->GetLinearD(0,0);
281 Float_t d1 = track1->GetLinearD(0,0);
282 //
283 // conservative cuts - convergence to be guarantied
284 // applying before track propagation
285 if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
286 if (dir>fCutMinDir) continue; // direction vector product
287 Float_t bz = AliTracker::GetBz();
288 Float_t dvertex0[2]; //distance to 0,0
289 Float_t dvertex1[2]; //distance to 0,0
290 track0->GetDZ(0,0,0,bz,dvertex0);
291 track1->GetDZ(0,0,0,bz,dvertex1);
292 if (TMath::Abs(dvertex0[1])>250) continue;
293 if (TMath::Abs(dvertex1[1])>250) continue;
294 //
295 //
296 //
297 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
298 AliExternalTrackParam param0(*track0);
299 AliExternalTrackParam param1(*track1);
300 //
301 // Propagate using Magnetic field and correct fo material budget
302 //
303 AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
304 AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
305 //
306 // Propagate rest to the 0,0 DCA - z should be ignored
307 //
308 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
309 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
310 //
311 param0.GetDZ(0,0,0,bz,dvertex0);
312 param1.GetDZ(0,0,0,bz,dvertex1);
313 //
314 Double_t xyz0[3];//,pxyz0[3];
315 Double_t xyz1[3];//,pxyz1[3];
316 param0.GetXYZ(xyz0);
317 param1.GetXYZ(xyz1);
318 Bool_t isPair = IsPair(&param0,&param1);
319 //
320 if (fStreamLevel>0){
321 TTreeSRedirector * cstream = GetDebugStreamer();
322 printf("My stream=%p\n",(void*)cstream);
323 if (cstream) {
324 (*cstream) << "Track0" <<
325 "dir="<<dir<< // direction
326 "OK="<<isPair<< // will be accepted
327 "b0="<<b0<< // propagate status
328 "b1="<<b1<< // propagate status
329 "Orig0.=" << track0 << // original track 0
330 "Orig1.=" << track1 << // original track 1
331 "Tr0.="<<&param0<< // track propagated to the DCA 0,0
332 "Tr1.="<<&param1<< // track propagated to the DCA 0,0
333 "v00="<<dvertex0[0]<< // distance using kalman
334 "v01="<<dvertex0[1]<< //
335 "v10="<<dvertex1[0]<< //
336 "v11="<<dvertex1[1]<< //
337 "d0="<<d0<< // linear distance to 0,0
338 "d1="<<d1<< // linear distance to 0,0
339 //
340 "x00="<<xyz0[0]<<
341 "x01="<<xyz0[1]<<
342 "x02="<<xyz0[2]<<
343 //
344 "x10="<<xyz1[0]<<
345 "x11="<<xyz1[1]<<
346 "x12="<<xyz1[2]<<
347 //
348 "Seed0.=" << track0 << // original seed 0
349 "Seed1.=" << track1 << // original seed 1
350 "dedx0="<<dedx0<< // dedx0
351 "dedx1="<<dedx1<< // dedx1
352 "\n";
353 }
354 }
355 }
356 }
357}
358
04087794 359/*
360
361
362void AliTPCcalibCosmic::dEdxCorrection(){
363 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
364 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
365 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
47328338 366 TCut cutN("cutN","min(Orig0.fTPCncls,Orig1.fTPCncls)>70");
04087794 367}
9b27d39b 368
04087794 369*/
f7f33dec 370
f7f33dec 371
9b27d39b 372
373Long64_t AliTPCcalibCosmic::Merge(TCollection */*li*/) {
374
f7f33dec 375}
376
9b27d39b 377Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
378 //
379 //
380 /*
381 // 0. Same direction - OPOSITE - cutDir +cutT
382 TCut cutDir("cutDir","dir<-0.99")
383 // 1.
384 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
385 //
386 // 2. The same rphi
387 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
388 //
389 //
390 //
391 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
392 // 1/Pt diff cut
393 */
394 const Double_t *p0 = tr0->GetParameter();
395 const Double_t *p1 = tr1->GetParameter();
396 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
397 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
398 Double_t d0[3], d1[3];
399 tr0->GetDirection(d0);
400 tr1->GetDirection(d1);
401 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
402 //
403 return kTRUE;
404}
405
406
f7f33dec 407
408void AliTPCcalibCosmic::BinLogX(TH1 *h) {
409
410 // Method for the correct logarithmic binning of histograms
411
412 TAxis *axis = h->GetXaxis();
413 int bins = axis->GetNbins();
414
415 Double_t from = axis->GetXmin();
416 Double_t to = axis->GetXmax();
417 Double_t *new_bins = new Double_t[bins + 1];
418
419 new_bins[0] = from;
420 Double_t factor = pow(to/from, 1./bins);
421
422 for (int i = 1; i <= bins; i++) {
423 new_bins[i] = factor * new_bins[i-1];
424 }
425 axis->Set(bins, new_bins);
426 delete new_bins;
427
428}
429