1. Making TPC calibration task working on PROOF
[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 {  
54   AliInfo("Defualt Constructor");  
55 }
56
57
58 AliTPCcalibCosmic::AliTPCcalibCosmic(const Text_t *name, const Text_t *title) 
59   :AliTPCcalibBase(),
60   fHistNTracks(0),
61   fClusters(0),
62   fModules(0),
63   fHistPt(0),
64   fPtResolution(0),
65   fDeDx(0)
66 {  
67   SetName(name);
68   SetTitle(title);
69   AliMagFMaps * field = new AliMagFMaps("dummy1", "dummy2",0,5,0);
70   AliTracker::SetFieldMap(field, kTRUE);  
71   fHistNTracks = new TH1F("ntracks","Number of Tracks per Event",501,-0.5,500.5);
72   fClusters = new TH1F("signal","Number of Clusters per track",160,0,160);
73   fModules = new TH2F("sector","Acorde hits; z (cm); x(cm)",1200,-1200,1200,600,-1000,1000);
74   fHistPt = new TH1F("Pt","Pt distribution",2000,0,50);  
75   fPtResolution = new TH1F("PtResolution","Pt resolution",100,-50,50);
76   fDeDx = new TH2F("DeDx","dEdx",500,0.01,20.,500,0.,500);
77   BinLogX(fDeDx);
78   AliInfo("Non Defualt Constructor");  
79 }
80
81 AliTPCcalibCosmic::~AliTPCcalibCosmic(){
82   //
83   //
84   //
85 }
86
87
88 void AliTPCcalibCosmic::Process(AliESDEvent *event) {
89     
90   if (!event) {
91     Printf("ERROR: ESD not available");
92     return;
93   }
94   
95   AliESDfriend *ESDfriend=static_cast<AliESDfriend*>(event->FindListObject("AliESDfriend"));
96   if (!ESDfriend) {
97    Printf("ERROR: ESDfriend not available");
98    return;
99   }
100
101   printf("Hallo world: Im here\n");
102   
103   Int_t n=event->GetNumberOfTracks(); 
104   fHistNTracks->Fill(n);
105   
106   //track loop
107   for (Int_t i=0;i<n;++i) { 
108    AliESDtrack *track = event->GetTrack(i); 
109    fClusters->Fill(track->GetTPCNcls());
110    
111    AliExternalTrackParam * trackIn = new AliExternalTrackParam(*track->GetInnerParam());
112    
113    AliESDfriendTrack *friendTrack = ESDfriend->GetTrack(i);
114    TObject *calibObject;
115    AliTPCseed *seed = 0;
116    for (Int_t l=0;calibObject=friendTrack->GetCalibObject(l);++l) {
117     if (seed=dynamic_cast<AliTPCseed*>(calibObject)) break;
118    }
119  
120    if (seed && track->GetTPCNcls() > 80) fDeDx->Fill(trackIn->GetP(), seed->CookdEdxNorm(0.05,0.45,0));
121
122     
123   }
124   
125   // dE/dx,pt and ACORDE study --> studies which need the pair selection  
126   if (n > 2) return;
127   
128   for (Int_t i=0;i<n;++i) {
129     AliESDtrack *track1 = event->GetTrack(i);
130      
131     Double_t d1[3];
132     track1->GetDirection(d1);
133     
134     for (Int_t j=i+1;j<n;++j) {
135      AliESDtrack *track2 = event->GetTrack(j);   
136      Double_t d2[3];
137      track2->GetDirection(d2);
138        
139      if (d1[0]*d2[0] + d1[1]*d2[1] + d1[2]*d2[2] < -0.999) {
140      
141       /*___________________________________ Pt resolution ________________________________________*/
142       if (track1->Pt() != 0 && track1->GetTPCNcls() > 80 && track2->GetTPCNcls() > 80) {
143        Double_t res = (track1->Pt() - track2->Pt());
144        res = res/(2*(track1->Pt() + track2->Pt()));
145        fPtResolution->Fill(100*res);
146       }
147       
148       /*_______________________________ Propagation to ACORDE ___________________________________*/
149       const Double_t AcordePlane = 850.; //distance of the central Acorde detectors to the beam line at y =0
150       const Double_t roof = 210.5; // distance from x =0 to end of magnet roof
151      
152       if (d1[1] > 0 && d2[1] < 0 && track1->GetTPCNcls() > 50) {        
153        Double_t r[3];
154        track1->GetXYZ(r);
155        Double_t x,z;
156        z = r[2] + (d1[2]/d1[1])*(AcordePlane - r[1]);
157        x = r[0] + (d1[0]/d1[1])*(AcordePlane - r[1]);
158        
159        if (x > roof) {
160         x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
161         z = z - TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
162        }
163        if (x < -roof) {
164         x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
165         z = z -  TMath::Abs(TMath::Tan(track1->Phi()))/TMath::Abs(TMath::Tan(track1->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track1->Phi())));
166        }
167        
168        fModules->Fill(z, x);
169       }
170       
171       if (d2[1] > 0 && d1[1] < 0 && track2->GetTPCNcls() > 50) {
172        Double_t r[3];
173        track2->GetXYZ(r);
174        Double_t x,z;
175        z = r[2] + (d2[2]/d2[1])*(AcordePlane - r[1]);
176        x = r[0] + (d2[0]/d2[1])*(AcordePlane - r[1]);
177        
178        if (x > roof) {
179         x = x - (x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
180         z = z - TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x-roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));  
181        }
182        if (x < -roof) {
183         x = x - (x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
184         z = z -  TMath::Abs(TMath::Tan(track2->Phi()))/TMath::Abs(TMath::Tan(track2->Theta()))*(x+roof)/(1 + TMath::Abs(TMath::Tan(track2->Phi())));
185        }       
186        
187        fModules->Fill(z, x);
188       }
189       
190       printf("My stream level=%d\n",fStreamLevel);
191
192       if (fStreamLevel>0){
193         TTreeSRedirector * cstream =  GetDebugStreamer();
194         printf("My stream=%p\n",cstream);
195         if (cstream) {
196           (*cstream) << "Track" <<
197             "Track1.=" << track1 <<      // original track 1
198             "Track2.=" << track2 <<      // original track2
199             "\n";
200         }
201       }
202   //     AliExternalTrackParam * trackOut = new AliExternalTrackParam(*track2->GetOuterParam());
203 //       AliTracker::PropagateTrackTo(trackOut,850.,105.658,30);
204 //       delete trackOut;
205       
206
207
208       
209
210       break;            
211      }     
212     }
213   }
214   
215   
216   
217   
218 }    
219
220
221 Long64_t AliTPCcalibCosmic::Merge(TCollection *li) {
222
223 }
224
225
226 void AliTPCcalibCosmic::BinLogX(TH1 *h) {
227
228   // Method for the correct logarithmic binning of histograms
229
230   TAxis *axis = h->GetXaxis();
231   int bins = axis->GetNbins();
232
233   Double_t from = axis->GetXmin();
234   Double_t to = axis->GetXmax();
235   Double_t *new_bins = new Double_t[bins + 1];
236    
237   new_bins[0] = from;
238   Double_t factor = pow(to/from, 1./bins);
239   
240   for (int i = 1; i <= bins; i++) {
241    new_bins[i] = factor * new_bins[i-1];
242   }
243   axis->Set(bins, new_bins);
244   delete new_bins;
245   
246 }
247