]> git.uio.no Git - u/mrichter/AliRoot.git/blame - T0/readESDlocal.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
bugfix
[u/mrichter/AliRoot.git] / T0 / readESDlocal.C
CommitLineData
fd154cff 1#include "TH1F.h"
2#include "TProfile.h"
3#include "TH2F.h"
4#include "TMap.h"
5#include "TFile.h"
6#include "TTree.h"
7#include "TBranch.h"
8
9#include "TGLabel.h"
10#include "TGrid.h"
11#include "TFitResult.h"
12
13#include "TControlBar.h"
14#include <Riostream.h>
15#include "AliCDBManager.h"
16#include "AliRawReader.h"
17#include "AliRawReaderRoot.h"
18#include "AliT0LookUpValue.h"
19#include "AliT0LookUpKey.h"
20#include "AliT0Parameters.h"
21#include "AliT0RawReader.h"
22#include "AliT0CalibLaserData.h"
23#include "AliESDTZERO.h"
24#include "AliESDVertex.h"
25#include "AliESDVZERO.h"
26#include "AliESDEvent.h"
27
28void readESDlocal(Int_t run)
29{
30 //read T0 ESD from reconstructed T0 data
31 TH1F *hMean= new TH1F("hMean"," T0 ", 1000,2000,3000);
32 // TH1F *hMeanTOF= new TH1F("hMeanTOF"," T0 ", 1000,-10,10);
33 // TH1F *hOrA= new TH1F("hOrA"," OrA T0 ",1000,-10,10);
34 // TH1F *hOrC= new TH1F("hOrC"," OrC T0 ",1000,-10,10);
35 TH1F *hVertex = new TH1F("hVertex","Z position of vertex", 180,-100,100);
36 TH1F *hVertexcalc = new TH1F("hVertexcalc","Z position of vertex", 180,-100,100);
37 TH1F *hVertexT0 = new TH1F("hVertexT0","Z position of vertex", 1000,-10,10);
38 TH1F *hVertexT0calc = new TH1F("hVertexT0calc","Z position of vertex", 1000,-10,10);
39 TH1F *hVertexComp = new TH1F("hVertexComparison ","Z position of vertex", 180,-30,30);
40 TH1F *hVertexT0only = new TH1F("hVertexT0only","Z position of vertex without SPD vertex", 180,-30,30);
41 TH1F *hVertexSPD = new TH1F("hVertexSPD","Z position of vertex SPD", 180,-30,30);
42
43 TH2F *hVertexSPDT0 = new TH2F("hVertexSPDT0","Z position of vertex SPD-T0", 180,-30,30, 180,-30,30);
44 //
45 TH1F *hMeanTOF= new TH1F("hMeanTOF"," T0 ", 8000, -100,100);
46 TH1F *hOrA= new TH1F("hOrA"," OrA T0 ",8000, -100,100);
47 TH1F *hOrC= new TH1F("hOrC"," OrC T0 ",8000, -100,100);
48 // TH1F *hMeanTOF= new TH1F("hMeanTOF"," T0 ", 2000, -8950,-8900);
49 // TH1F *hOrA= new TH1F("hOrA"," OrA T0 ",2000, -8950,-8900);
50 // TH1F *hOrC= new TH1F("hOrC"," OrC T0 ",2000, -8950,-8900);
51
52
53 TH1F *hMeanTOFcalc= new TH1F("hMeanTOFcalc"," T0 ", 2000, -8950,-8900);
54 TH1F *hOrAcalc= new TH1F("hOrAcalc"," OrA T0 ",2000, -8950,-8900);
55 TH1F *hOrCcalc= new TH1F("hOrCcalc"," OrC T0 ",2000, -8950,-8900);
56 // TH1F *hMeanTOFcalc= new TH1F("hMeanTOFcalc"," T0 ", 8000, -100,100);
57 // TH1F *hOrAcalc= new TH1F("hOrAcalc"," OrA T0 ",8000, -100,100);
58 // TH1F *hOrCcalc= new TH1F("hOrCcalc"," OrC T0 ",8000, -100,100);
59
60 TH1F *hAmp[24]; TH1F * hTime[24]; TH1F * hTimeDiff[24]; TH1F * hTimecorr[24];
61 // TProfile * hAmpTime[24];
62 TH2F * hAmpTime[24];
63
64 for(Int_t ic=0; ic<24; ic++)
65 {
66 hAmp[ic] = new TH1F(Form("hAmp%i",ic+1),"Amp",100, -10, 20 );
67 // hAmp[ic] = new TH1F(Form("hAmp%i",ic),"Amp",400, 300, 700 );
68 hTime[ic] = new TH1F(Form("hTime%i",ic+1)," time", 5000,0,10000);
69 hTimecorr[ic] = new TH1F(Form("hTimecorr%i",ic+1)," time", 5000,0,10000);
70 hTimeDiff[ic] = new TH1F(Form("hTimeDiff%i",ic+1)," time", 300,-300,300);
71 // hAmpTime[ic] = new TProfile(Form("hAmpTime%i",ic)," time",200, 300,700, 49000, 53000);
72 hAmpTime[ic] = new TH2F(Form("hAmpTime%i",ic+1)," time",100, 0,10, 100,3000,4000);
73 // hAmpTime[ic] = new TH2F(Form("hAmpTime%i",ic)," time",250, 350,600, 400, 49400, 49800);
74 }
75
76 Float_t channelWidth=24.4;
77 Float_t c = 0.0299792458; // cm/ps
78 Float_t shift=0;
79 Float_t fLatencyHPTDC = 9000;
80
81
82 //run 124702
83 // Float_t fLatencyL1 = 8.91358e+03;
84 // Float_t fLatencyL1A = 8.91352e+03;
85 // Float_t fLatencyL1C =8.91361e+03;
86
87 //run 125097
88 // Float_t fLatencyL1 = 8.914520e+03 ;
89 // Float_t fLatencyL1A = 8.914860e+03 ;
90 // Float_t fLatencyL1C = 8.914180e+03;
91
92//run 125295
93 Float_t fLatencyL1 = 8.91406e+03 ;
94 Float_t fLatencyL1A = 8.91401e+03;
95 Float_t fLatencyL1C = 8.91412e+03 ;
96
97
98 // Float_t fLatencyL1 = 0;
99 // Float_t fLatencyL1A = 0;
100 // Float_t fLatencyL1C =0;
101 //125842
102 // Float_t fLatencyL1 = 8.91306e+03;
103 // Float_t fLatencyL1A =8.91338e+03;
104 // Float_t fLatencyL1C =8.91274e+03;
105 //
106 //126407
107 // Float_t fLatencyL1 = 8.91345e+03;
108 // Float_t fLatencyL1A =8.91378e+03;
109 // Float_t fLatencyL1C =8.91311e+03;
110 //
111 const Double32_t *amp, *time, *mean;
112 Double32_t time1[24];
113 TString filenam[100];
114
115
116
117 // Float_t equal[24]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0 };
118 //124358
119 //Float_t equal[24]={22, 0 ,0 ,1 ,3 ,1 ,0 ,0 ,-4 ,1 ,1 ,-1 ,0 ,-3 ,0 ,-6 ,0 ,-3 ,2 ,0 ,-1 ,-1 ,-3 ,0 };
120
121 //124702
122 // Float_t equal[24]={20, 0 ,15, 50, 37, 36, 23, -20, 77, 54, 44, -12,
123 // 0, 19, -25, -13, 35, 1, 23, 46, 20, 77, 50, 122};
124
125
126 //125097
127 Float_t equal[24]={8,0 ,15 ,50 ,36 ,34 ,24 ,-21 ,81 ,53 ,42 ,-12 ,0 ,19 ,-20 ,-7 ,37 ,5 ,25 ,49 ,19 ,78 ,52 ,124};
128
129 // 125295
130
131 //Float_t equal[24]={16, 0, 2, 2, 30, 3, -5, -8, 5, -4, 16, -5,
132 // 0, 10, 10, 12, 13, -24, 15, 26, -2, 17, 10, -30};
133
134
135
136 //125842
137 // Float_t equal[24]={16, 0, 1, 3, 32, 4, -2, -8, 5, -3, 15, -5,
138 // 0, 10, 9, 12, 14, -23, 16, 26, 0, 19, 9, -31 };
139
140 //126407
141 // Float_t equal[24]={16, 0, 1, 4, 33, 4, 0, -7, 6, -3, 15, -4,
142 // 0 , 10, 10, 13, 13, -23, 16, 28, 0, 19, 11, -29};
143
144 //for ( Int_t indexfile=filestart; indexfile < filestop+1;indexfile++ )
145 for ( Int_t indexfile=0; indexfile < 1; indexfile++ )
146 {
147 TString fFileName=Form("%i/noeq/AliESDs.root", run);
148 // cout<<" filenam "<<fFileName<<endl;
149 // fFileName=filenam[indexfile];
150
151 TFile *file = TFile::Open(fFileName);
152 // cout<<file<<endl;
153
154 // TFile *file = new TFile("AliESDs.root");
155
156 TTree *esdTree = (TTree*)file->Get("esdTree");
157
158 AliESDEvent *esdevent = new AliESDEvent;
159 esdevent->ReadFromTree(esdTree);
160 cout<<"esdevent "<<esdevent<<endl;
161
162 TBranch *brSPD=esdTree->GetBranch("SPDVertex.");
163 AliESDVertex *fverSPD = new AliESDVertex();
164 if (brSPD) {
165 brSPD->SetAddress(&fverSPD);
166 }else{
167 cerr<<"EXEC Branch SPDvertex not found"<<endl;
168 return;
169 }
170
171 TBranch *brRec=esdTree->GetBranch("AliESDTZERO.");
172 AliESDTZERO *fesdT0 = new AliESDTZERO();
173 if (brRec) {
174 brRec->SetAddress(&fesdT0);
175 }else{
176 cerr<<"EXEC Branch T0 not found"<<endl;
177 return;
178 }
179 cout<<" branches "<<brRec<<endl;
180 // Event ------------------------- LOOP
181 for (Int_t ievent=0; ievent<brRec->GetEntries(); ievent++){
182
183 Double32_t fOrA, besttimeA=9999999;
184 Double32_t fOrC, besttimeC=9999999;
185 for (Int_t ip=0; ip<24; ip++) time1[ip]=0;
186 Float_t timecorr=0;
187 Int_t ncont=0;
188 Float_t shift=0;
189 brRec->GetEntry(ievent);
190 brSPD->GetEntry(ievent);
191 esdTree->GetEvent(ievent);
192
193 TString triggers = esdevent->GetFiredTriggerClasses();
194 // printf("Event %d: trigger classes: \"%s\"\n",
195 // esdevent->GetEventNumberInFile(),
196 // esdevent->GetFiredTriggerClasses().Data());
197 if ( !triggers.Contains("CINT1B-ABCE-NOPF-ALL") ) {
198 // Printf("Skip event with trigger class \"%s\"",triggers.Data());
199 continue;
200 }
201 Double_t spdver = fverSPD->GetZ();
202 ncont = fverSPD->GetNContributors();
203 // cout<<" spdver "<<spdver<<" ncont "<<ncont<<endl;
204 if(ncont>2) {
205 hVertexSPD->Fill(spdver);
206 shift = spdver/(c*channelWidth );
207 // cout<<ievent<<" vertex shif "<<shift<<" vertex "<<spdver<<" IsFromVertexer3D "<<fverSPD->IsFromVertexer3D()<<endl;
208 }
209
210 Int_t trig=fesdT0->GetT0Trig();
211 Float_t meanTOF = fesdT0->GetT0();
212 mean = fesdT0->GetT0TOF();
213 Double32_t vertex= fesdT0->GetT0zVertex();
214 Double32_t orA=0.001* mean[1] ;
215 Double32_t orC= 0.001*mean[2] ;
216 if (orA<99){
217 hOrA->Fill(orA);
218 // cout<<ievent<<" ORA "<< orA<<endl;
219 }
220 if (orC<99) {
221 hOrC->Fill(orC);
222 // cout<<ievent<<" ORC "<< orC<<endl;
223 }
224
225 if(orA<99 && orC<99) hVertexT0->Fill((orA-orC)/2.);
226
227 if (vertex<99990) {
228 hMeanTOF->Fill(0.001*mean[0] );
229 hMean->Fill(meanTOF);
230 hVertex->Fill(vertex);
231 if(ncont>1) {
232 hVertexComp->Fill(vertex-spdver);
233 hVertexSPDT0->Fill(vertex,spdver);
234 }
235 if(ncont<0) hVertexT0only->Fill(vertex);
236
237 }
238 amp=fesdT0->GetT0amplitude();
239 time=fesdT0->GetT0time();
240 for (Int_t i=0; i<24; i++){
241 if( time[i]>0){
242 hAmp[i]->Fill(amp[i]);
243 hTime[i]->Fill(time[i]);
244 time1[i]=time[i] -equal[i];
245
246 if(i<12){
247 if(time[1] >0) hTimeDiff[i]->Fill((time1[i] - time1[1]));
248 if(ncont>2) { timecorr = time[i] - shift;
249 hTimecorr[i]->Fill(timecorr);
250 hAmpTime[i]->Fill(amp[i],timecorr);
251 }
252 }
253
254
255 if(i>11) {
256 if(time[12] >0) hTimeDiff[i]->Fill((time1[i]- time1[12]));
257 if(ncont>2) { timecorr = time[i] + shift;
258 hTimecorr[i]->Fill(timecorr);
259 hAmpTime[i]->Fill(amp[i],timecorr);
260 }
261 }
262 // cout<<ievent<<" "<<i<<" "<<<endl;
263 }
264 }
265
266 for (Int_t ipmt=0; ipmt<12; ipmt++){
267
268 if( time1[ipmt]>1 && time1[ipmt] < besttimeC ){
269 besttimeC=time1[ipmt]; //timeC
270
271 }
272 }
273 for ( Int_t ipmt=12; ipmt<23; ipmt++){
274 if( time1[ipmt]>1 && time1[ipmt]<besttimeA ) {
275 besttimeA=time1[ipmt]; //timeA
276 }
277 }
278 if(besttimeA<9999999 ) {
279 fOrA=0.001*besttimeA * channelWidth - fLatencyHPTDC + fLatencyL1A +0.001*shift*channelWidth;
280 hOrAcalc->Fill(fOrA);
281 // cout<<ievent<<" calc OrA "<<fOrA<<endl;
282 }
283 if( besttimeC < 9999999 ) {
284 fOrC=0.001*besttimeC * channelWidth - fLatencyHPTDC +fLatencyL1C -0.001*shift*channelWidth;
285 hOrCcalc->Fill(fOrC);
286 // cout<<ievent<<" calc OrC "<<fOrC<<endl;
287 }
288 if(besttimeA <9999999 && besttimeC < 9999999 )
289 {
290 Double_t meanCalc= 0.001*channelWidth * Float_t( besttimeA+besttimeC)/2.- fLatencyHPTDC + fLatencyL1;
291 hMeanTOFcalc->Fill(meanCalc );
292 hVertexT0calc->Fill((fOrA-fOrC)/2.);
293 Float_t timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth + fLatencyL1A - fLatencyL1C;
294 Float_t vertexcal = - 1000*c*(timeDiff)/2. ; //+ (fdZonA - fdZonC)/2;
295
296 hVertexcalc->Fill(vertexcal);
297
298
299 // cout<<ievent<<" calc meanTOF "<<meanCalc<<endl;
300 }
301
302 }
303
304 }
305 TFile* Hfile = new TFile(Form("FigESD_%i.noeq.root",run),"RECREATE","Histograms for T0 digits");
306 printf("Writting histograms to root file \n");
307 Hfile->cd();
308 //Create a canvas, set the view range, show histograms
309 //
310 Float_t differ;
311 hMean->Write();
312 hMeanTOF->Write();
313 hMeanTOFcalc->Write();
314 hVertex->Write();
315 hVertexT0->Write();
316 hVertexT0calc->Write();
317 hVertexSPD->Write();
318 hVertexComp->Write();
319 hVertexT0only->Write();
320 hOrC->Write();
321 hOrA->Write();
322 hOrCcalc->Write();
323 hOrAcalc->Write();
324 hVertexSPDT0->Write();
325 hVertexcalc->Write();
326 for (Int_t i=0; i<24; i++){
327 hAmp[i]->Write();
328 hTime[i]->Write();
329 TFitResultPtr r = hTimeDiff[i]->Fit("gaus","SQ");
330 Double_t par0;
331 if ((Int_t)r == 0) par0 = r->Parameters()[1];
332 else par0 = 99999;
333 hTimeDiff[i]->Write();
334 cout<<Int_t(par0)<<", ";
335 // cout<<Int_t(hTimeDiff[i]->GetMean())<<" fit "<<par0<<endl;
336 // cout<<Int_t(hTimeDiff[i]->GetMean())<<" , ";
337 hAmpTime[i]->Write();
338 hTimecorr[i]->Write();
339 }
340 cout<<endl;
341 for (Int_t i=1; i<24; i++) cout<<i<<" ,";
342 cout<<endl;
343 for (Int_t i=1; i<24; i++) cout<<Int_t(hTimeDiff[i]->GetMean())<<" ,";
344 cout<<endl;
345
346}