1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 #include "Riostream.h"
26 #include "AliTPCseed.h"
27 #include "AliESDVertex.h"
28 #include "AliESDEvent.h"
29 #include "AliESDfriend.h"
30 #include "AliESDInputHandler.h"
32 #include "AliTracker.h"
33 #include "AliMagFMaps.h"
37 #include "AliTPCcalibCosmic.h"
39 #include "TTreeStream.h"
40 #include "AliTPCTracklet.h"
42 ClassImp(AliTPCcalibCosmic)
45 AliTPCcalibCosmic::AliTPCcalibCosmic()
53 fCutMaxD(5), // maximal distance in rfi ditection
54 fCutTheta(0.03), // maximal distan theta
55 fCutMinDir(-0.99) // direction vector products
57 AliInfo("Defualt Constructor");
61 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title)
69 fCutMaxD(5), // maximal distance in rfi ditection
70 fCutTheta(0.03), // maximal distan theta
71 fCutMinDir(-0.99) // direction vector products
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);
84 AliInfo("Non Default Constructor");
87 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
97 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
102 Printf("ERROR: ESD not available");
105 AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
107 Printf("ERROR: ESDfriend not available");
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;
120 for (Int_t i=0;i<ntracks;++i) {
121 AliESDtrack *track = event->GetTrack(i);
122 fClusters->Fill(track->GetTPCNcls());
123 AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
125 AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
126 TObject *calibObject;
127 AliTPCseed *seed = 0;
128 for (Int_t l=0;(calibObject=friendTrack->GetCalibObject(l));++l) {
129 if ((seed=dynamic_cast<AliTPCseed*>(calibObject))) break;
131 if (seed) tpcSeeds.AddAt(seed,i);
132 if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
134 if (ntracks<2) return;
137 // dE/dx,pt and ACORDE study --> studies which need the pair selection
138 for (Int_t i=0;i<ntracks;++i) {
139 AliESDtrack *track1 = event->GetTrack(i);
142 track1->GetDirection(d1);
144 for (Int_t j=i+1;j<ntracks;++j) {
145 AliESDtrack *track2 = event->GetTrack(j);
147 track2->GetDirection(d2);
149 if (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2] < -0.999) {
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);
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
162 if (d1[1] > 0 && d2[1] < 0 && track1->GetTPCNcls() > 50) {
166 z = r[2] + (d1[2]/d1[1])*(AcordePlane - r[1]);
167 x = r[0] + (d1[0]/d1[1])*(AcordePlane - r[1]);
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())));
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())));
178 fModules->Fill(z, x);
181 if (d2[1] > 0 && d1[1] < 0 && track2->GetTPCNcls() > 50) {
185 z = r[2] + (d2[2]/d2[1])*(AcordePlane - r[1]);
186 x = r[0] + (d2[0]/d2[1])*(AcordePlane - r[1]);
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())));
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())));
197 fModules->Fill(z, x);
200 // AliExternalTrackParam * trackOut = new AliExternalTrackParam(*track2->GetOuterParam());
201 // AliTracker::PropagateTrackTo(trackOut,850.,105.658,30);
218 void AliTPCcalibCosmic::FindPairs(AliESDEvent *event) {
222 // Track0 is choosen in upper TPC part
223 // Track1 is choosen in lower TPC part
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();
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);
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());
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;
247 if (seed) tpcSeeds.AddAt(seed,i);
248 if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
250 if (ntracks<2) return;
254 for (Int_t i=0;i<ntracks;++i) {
255 AliESDtrack *track0 = event->GetTrack(i);
256 // track0 - choosen upper part
257 if (!track0) continue;
258 if (!track0->GetOuterParam()) continue;
259 if (track0->GetOuterParam()->GetAlpha()<0) continue;
261 track0->GetDirection(d1);
262 for (Int_t j=0;j<ntracks;++j) {
264 AliESDtrack *track1 = event->GetTrack(j);
266 if (!track1) continue;
267 if (!track1->GetOuterParam()) continue;
268 if (track1->GetOuterParam()->GetAlpha()>0) continue;
271 track1->GetDirection(d2);
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);
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;
297 Float_t dmax = TMath::Max(TMath::Abs(d0),TMath::Abs(d1));
298 AliExternalTrackParam param0(*track0);
299 AliExternalTrackParam param1(*track1);
301 // Propagate using Magnetic field and correct fo material budget
303 AliTracker::PropagateTrackTo(¶m0,dmax+1,0.0005,3,kTRUE);
304 AliTracker::PropagateTrackTo(¶m1,dmax+1,0.0005,3,kTRUE);
306 // Propagate rest to the 0,0 DCA - z should be ignored
308 Bool_t b0 = param0.PropagateToDCA(&vtx,bz,1000);
309 Bool_t b1 = param1.PropagateToDCA(&vtx,bz,1000);
311 param0.GetDZ(0,0,0,bz,dvertex0);
312 param1.GetDZ(0,0,0,bz,dvertex1);
314 Double_t xyz0[3];//,pxyz0[3];
315 Double_t xyz1[3];//,pxyz1[3];
318 Bool_t isPair = IsPair(¶m0,¶m1);
321 TTreeSRedirector * cstream = GetDebugStreamer();
322 printf("My stream=%p\n",(void*)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.="<<¶m0<< // track propagated to the DCA 0,0
332 "Tr1.="<<¶m1<< // 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
348 "Seed0.=" << track0 << // original seed 0
349 "Seed1.=" << track1 << // original seed 1
350 "dedx0="<<dedx0<< // dedx0
351 "dedx1="<<dedx1<< // dedx1
362 void 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");
373 Long64_t AliTPCcalibCosmic::Merge(TCollection */*li*/) {
377 Bool_t AliTPCcalibCosmic::IsPair(AliExternalTrackParam *tr0, AliExternalTrackParam *tr1){
381 // 0. Same direction - OPOSITE - cutDir +cutT
382 TCut cutDir("cutDir","dir<-0.99")
384 TCut cutT("cutT","abs(Tr1.fP[3]+Tr0.fP[3])<0.03")
387 TCut cutD("cutD","abs(Tr0.fP[0]+Tr1.fP[0])<5")
391 TCut cutPt("cutPt","abs(Tr1.fP[4]+Tr0.fP[4])<1&&abs(Tr0.fP[4])+abs(Tr1.fP[4])<10");
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;
408 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
410 // Method for the correct logarithmic binning of histograms
412 TAxis *axis = h->GetXaxis();
413 int bins = axis->GetNbins();
415 Double_t from = axis->GetXmin();
416 Double_t to = axis->GetXmax();
417 Double_t *new_bins = new Double_t[bins + 1];
420 Double_t factor = pow(to/from, 1./bins);
422 for (int i = 1; i <= bins; i++) {
423 new_bins[i] = factor * new_bins[i-1];
425 axis->Set(bins, new_bins);