1. Adding function to find cosmic track pairs
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibCosmic.cxx
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
42 ClassImp(AliTPCcalibCosmic)
43
44
45 AliTPCcalibCosmic::AliTPCcalibCosmic() 
46   :AliTPCcalibBase(),
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
56 {  
57   AliInfo("Defualt Constructor");  
58 }
59
60
61 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title) 
62   :AliTPCcalibBase(),
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
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);
84   AliInfo("Non Default Constructor");  
85 }
86
87 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
88   //
89   //
90   //
91 }
92
93
94
95
96
97 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
98   //
99   //
100   //
101   if (!event) {
102     Printf("ERROR: ESD not available");
103     return;
104   }  
105   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
106   if (!ESDfriend) {
107    Printf("ERROR: ESDfriend not available");
108    return;
109   }
110   FindPairs(event);
111
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   //
118   //track loop
119   //
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());
124    
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;
130    }
131    if (seed) tpcSeeds.AddAt(seed,i);
132    if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0)); 
133   }
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) {
139     AliESDtrack *track1 = event->GetTrack(i);
140      
141     Double_t d1[3];
142     track1->GetDirection(d1);
143     
144     for (Int_t j=i+1;j<ntracks;++j) {
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       
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
218 void 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
349
350
351
352 Long64_t AliTPCcalibCosmic::Merge(TCollection */*li*/) {
353   
354 }
355
356 Bool_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
386
387 void 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