1. Adding function to find cosmic track pairs
[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
221 //
222 //
223 if (GetDebugLevel()>1) printf("Hallo world: Im here\n");
224 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
225 Int_t ntracks=event->GetNumberOfTracks();
226 fHistNTracks->Fill(ntracks);
227 TObjArray tpcSeeds(ntracks);
228 if (ntracks==0) return;
229 Double_t vtxx[3]={0,0,0};
230 Double_t svtxx[3]={0.000001,0.000001,100.};
231 AliESDVertex vtx(vtxx,svtxx);
232 //
233 //track loop
234 //
235 for (Int_t i=0;i<ntracks;++i) {
236 AliESDtrack *track = event->GetTrack(i);
237 fClusters->Fill(track->GetTPCNcls());
238 AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
239
240 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
241 TObject *calibObject;
242 AliTPCseed *seed = 0;
243 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
244 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
245 }
246 if (seed) tpcSeeds.AddAt(seed,i);
247 if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
248 }
249 if (ntracks<2) return;
250 //
251 // Find pairs
252 //
253 for (Int_t i=0;i<ntracks;++i) {
254 AliESDtrack *track0 = event->GetTrack(i);
255 Double_t d1[3];
256 track0->GetDirection(d1);
257 for (Int_t j=i+1;j<ntracks;++j) {
258 AliESDtrack *track1 = event->GetTrack(j);
259 Double_t d2[3];
260 track1->GetDirection(d2);
261 printf("My stream level=%d\n",fStreamLevel);
262 AliTPCseed * seed0 = (AliTPCseed*) tpcSeeds.At(i);
263 AliTPCseed * seed1 = (AliTPCseed*) tpcSeeds.At(j);
264 if (! seed0) continue;
265 if (! seed1) continue;
266 Float_t dedx0 = seed0->CookdEdxNorm(0.05,0.55,0);
267 Float_t dedx1 = seed1->CookdEdxNorm(0.05,0.55,0);
268 Float_t dir = (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2]);
269 Float_t d0 = track0->GetLinearD(0,0);
270 Float_t d1 = track1->GetLinearD(0,0);
271 //
272 // conservative cuts - convergence to be guarantied
273 // applying before track propagation
274 if (TMath::Abs(d0+d1)>fCutMaxD) continue; // distance to the 0,0
275 if (dir>fCutMinDir) continue; // direction vector product
276 Float_t bz = AliTracker::GetBz();
277 Float_t dvertex0[2]; //distance to 0,0
278 Float_t dvertex1[2]; //distance to 0,0
279 track0->GetDZ(0,0,0,bz,dvertex0);
280 track1->GetDZ(0,0,0,bz,dvertex1);
281 if (TMath::Abs(dvertex0[1])>250) continue;
282 if (TMath::Abs(dvertex1[1])>250) continue;
283 //
284 //
285 //
286 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
287 AliExternalTrackParam param0(*track0);
288 AliExternalTrackParam param1(*track1);
289 //
290 // Propagate using Magnetic field and correct fo material budget
291 //
292 AliTracker::PropagateTrackTo(&param0,dmax+1,0.0005,3,kTRUE);
293 AliTracker::PropagateTrackTo(&param1,dmax+1,0.0005,3,kTRUE);
294 //
295 // Propagate rest to the 0,0 DCA - z should be ignored
296 //
297 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
298 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
299 //
300 param0.GetDZ(0,0,0,bz,dvertex0);
301 param1.GetDZ(0,0,0,bz,dvertex1);
302 //
303 Double_t xyz0[3];//,pxyz0[3];
304 Double_t xyz1[3];//,pxyz1[3];
305 param0.GetXYZ(xyz0);
306 param1.GetXYZ(xyz1);
307 Bool_t isPair = IsPair(&param0,&param1);
308 //
309 if (fStreamLevel>0){
310 TTreeSRedirector * cstream = GetDebugStreamer();
311 printf("My stream=%p\n",(void*)cstream);
312 if (cstream) {
313 (*cstream) << "Track0" <<
314 "dir="<<dir<< // direction
315 "OK="<<isPair<< // will be accepted
316 "b0="<<b0<< // propagate status
317 "b1="<<b1<< // propagate status
318 "Orig0.=" << track0 << // original track 0
319 "Orig1.=" << track1 << // original track 1
320 "Tr0.="<<&param0<< // track propagated to the DCA 0,0
321 "Tr1.="<<&param1<< // track propagated to the DCA 0,0
322 "v00="<<dvertex0[0]<< // distance using kalman
323 "v01="<<dvertex0[1]<< //
324 "v10="<<dvertex1[0]<< //
325 "v11="<<dvertex1[1]<< //
326 "d0="<<d0<< // linear distance to 0,0
327 "d1="<<d1<< // linear distance to 0,0
328 //
329 "x00="<<xyz0[0]<<
330 "x01="<<xyz0[1]<<
331 "x02="<<xyz0[2]<<
332 //
333 "x10="<<xyz1[0]<<
334 "x11="<<xyz1[1]<<
335 "x12="<<xyz1[2]<<
336 //
337 "Seed0.=" << track0 << // original seed 0
338 "Seed1.=" << track1 << // original seed 1
339 "dedx0="<<dedx0<< // dedx0
340 "dedx1="<<dedx1<< // dedx1
341 "\n";
342 }
343 }
344 }
345 }
346}
347
348
f7f33dec 349
f7f33dec 350
9b27d39b 351
352Long64_t AliTPCcalibCosmic::Merge(TCollection */*li*/) {
353
f7f33dec 354}
355
9b27d39b 356Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
357 //
358 //
359 /*
360 // 0. Same direction - OPOSITE - cutDir +cutT
361 TCut cutDir("cutDir","dir<-0.99")
362 // 1.
363 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
364 //
365 // 2. The same rphi
366 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
367 //
368 //
369 //
370 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
371 // 1/Pt diff cut
372 */
373 const Double_t *p0 = tr0->GetParameter();
374 const Double_t *p1 = tr1->GetParameter();
375 if (TMath::Abs(p0[3]+p1[3])>fCutTheta) return kFALSE;
376 if (TMath::Abs(p0[0]+p1[0])>fCutMaxD) return kFALSE;
377 Double_t d0[3], d1[3];
378 tr0->GetDirection(d0);
379 tr1->GetDirection(d1);
380 if (d0[0]*d1[0] + d0[1]*d1[1] + d0[2]*d1[2] >fCutMinDir) return kFALSE;
381 //
382 return kTRUE;
383}
384
385
f7f33dec 386
387void AliTPCcalibCosmic::BinLogX(TH1 *h) {
388
389 // Method for the correct logarithmic binning of histograms
390
391 TAxis *axis = h->GetXaxis();
392 int bins = axis->GetNbins();
393
394 Double_t from = axis->GetXmin();
395 Double_t to = axis->GetXmax();
396 Double_t *new_bins = new Double_t[bins + 1];
397
398 new_bins[0] = from;
399 Double_t factor = pow(to/from, 1./bins);
400
401 for (int i = 1; i <= bins; i++) {
402 new_bins[i] = factor * new_bins[i-1];
403 }
404 axis->Set(bins, new_bins);
405 delete new_bins;
406
407}
408