]> git.uio.no Git - u/mrichter/AliRoot.git/blame - T0/MakeTrendT0.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
files for central QA
[u/mrichter/AliRoot.git] / T0 / MakeTrendT0.C
CommitLineData
b1a2430c 1#include "TChain.h"
2#include "TTree.h"
3#include "TH1F.h"
4#include "TF1.h"
5#include "TH2F.h"
6#include "TCanvas.h"
7#include "TObjArray.h"
8#include "/home/alla/alice/AliRoot/AliRoot/STEER/STEERBase/TTreeStream.h"
9
10#define NPMTs 24
11
12int MakeTrendT0( char *infile, int run) {
13
14 char *outfile = "trending.root";
15
16 if(!infile) return -1;
17 if(!outfile) return -1;
18
19 TFile *f = TFile::Open(infile,"read");
20 if (!f) {
21 printf("File %s not available\n", infile);
22 return -1;
23 }
24
25 // LOAD HISTOGRAMS FROM QAresults.root
26 TObjArray *fTzeroObject = (TObjArray*) f->Get("T0_Performance/QAT0chists");
27
28 TH1F* fTzeroORAplusORC =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORAplusORC"))->Clone("A");
29 TH1F* fResolution =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fResolution"))->Clone("B");
30 TH1F* fTzeroORA =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORA"))->Clone("C");
31 TH1F* fTzeroORC =(TH1F*) ((TH1F*) fTzeroObject->FindObject("fTzeroORC"))->Clone("D");
32
33 TH2F *fTimeVSAmplitude[NPMTs];//counting PMTs from 0
34 TH1D *fAmplitude[NPMTs];
35 TH1D *fTime[NPMTs];
36 //-----> add new histogram here
37
38
39
40 // sigma fit of resolution, mean T0A, T0C and T0C&A, mean amplitude for each PMT
41 double resolutionSigma = -9999; //dummy init
42 double tzeroOrAPlusOrC = -9999; //dummy init
43 double tzeroOrA = -9999; //dummy init
44 double tzeroOrC = -9999; //dummy init
45 double meanAmplitude[NPMTs]; //dummy init
46 double meanTime[NPMTs]; //dummy init
47 double timeDelayOCDB[NPMTs]; //dummy init
48
49 //.................................................................................
50
51 for(int ipmt=1; ipmt<=NPMTs; ipmt++){ //loop over all PMTs
52 fTimeVSAmplitude[ipmt-1] =(TH2F*) ((TH2F*) fTzeroObject->FindObject(Form("fTimeVSAmplitude%d",ipmt)))->Clone(Form("E%d",ipmt));
53 int nbinsX = fTimeVSAmplitude[ipmt-1]->GetNbinsX();
54 int nbinsY = fTimeVSAmplitude[ipmt-1]->GetNbinsY();
55 //Amplitude
56 fAmplitude[ipmt-1] = (TH1D*) fTimeVSAmplitude[ipmt-1]->ProjectionX(Form("fAmplitude%d",ipmt), 1, nbinsY);
57 meanAmplitude[ipmt-1] = -9999; //dummy init
58 if(fAmplitude[ipmt-1]->GetEntries()>0){
59 meanAmplitude[ipmt-1] = fAmplitude[ipmt-1]->GetMean();
60 }
61
62 //Time
63 fTime[ipmt-1] = (TH1D*) fTimeVSAmplitude[ipmt-1]->ProjectionY(Form("fTime%d",ipmt), 1, nbinsX);
64 meanTime[ipmt-1] = -9999; //dummy init
65 if(fTime[ipmt-1]->GetEntries()>0){
66 if(fTime[ipmt-1]->GetEntries()>20){
67 meanTime[ipmt-1] = GetParameterGaus((TH1F*) fTime[ipmt-1], 1); // Mean Time
68 }else if(fTime[ipmt-1]->GetEntries()>0){
69 meanTime[ipmt-1] = fTime[ipmt-1]->GetMean();
70 }
71 }
72 }
73
74
75
76 if(fResolution->GetEntries()>20){
77 resolutionSigma = GetParameterGaus(fResolution, 2); //gaussian sigma
78 }else if(fResolution->GetEntries()>0){
79 resolutionSigma = fResolution->GetRMS(); //gaussian sigma
80 }
81
82 if(fTzeroORAplusORC->GetEntries()>20){
83 tzeroOrAPlusOrC = GetParameterGaus(fTzeroORAplusORC, 1); //gaussian mean
84 }else if(fTzeroORAplusORC->GetEntries()>0){
85 tzeroOrAPlusOrC = fTzeroORAplusORC->GetMean();
86 }
87
88 if(fTzeroORA->GetEntries()>20){
89 tzeroOrA = GetParameterGaus(fTzeroORA, 1); //gaussian mean
90 }else if(fTzeroORA->GetEntries()>0){
91 tzeroOrA = fTzeroORA->GetMean();
92 }
93
94 if(fTzeroORC->GetEntries()>20){
95 tzeroOrC = GetParameterGaus(fTzeroORC, 1); //gaussian mean
96 }else if(fTzeroORC->GetEntries()>0){
97 tzeroOrC = fTzeroORC->GetMean(); //gaussian mean
98 }
99
100 //-----> analyze the new histogram here and set mean/sigma
101
102
103 f->Close();
104
105 //-------------------- READ OCDB TIME DELAYS ---------------------------
106 // Arguments:
107 AliCDBManager* man = AliCDBManager::Instance();
108 man->SetDefaultStorage("raw://");
109 man->SetRun(run);
110 AliCDBEntry *entry = AliCDBManager::Instance()->Get("T0/Calib/TimeDelay");
111 AliT0CalibTimeEq *clb = (AliT0CalibTimeEq*)entry->GetObject();
112 for (Int_t i=0; i<NPMTs; i++){
113 timeDelayOCDB[i] = 0;
114 if(clb)
115 timeDelayOCDB[i] = clb->GetCFDvalue(i,0);
116 }
117
118
119
120 //--------------- write walues to the output ------------
121 TTreeSRedirector* pcstream = NULL;
122 pcstream = new TTreeSRedirector(outfile);
123 if (!pcstream) return;
124
125
126 TFile *x = pcstream->GetFile();
127 x->cd();
128
129 TObjString runType;
130 Int_t startTimeGRP=0;
131 Int_t stopTimeGRP=0;
132 Int_t time=0;
133 Int_t duration=0;
134
135 time = (startTimeGRP+stopTimeGRP)/2;
136 duration = (stopTimeGRP-startTimeGRP);
137
138 (*pcstream)<<"t0QA"<<
139 "run="<<run;//<<
140 //"time="<<time<<
141 // "startTimeGRP="<<startTimeGRP<<
142 // "stopTimeGRP="<<stopTimeGRP<<
143 /// "duration="<<duration<<
144 // "runType.="<<&runType;
145
146
147 (*pcstream)<<"t0QA"<<
148 "resolution="<< resolutionSigma<<
149 "tzeroOrAPlusOrC="<< tzeroOrAPlusOrC<<
150 "tzeroOrA="<< tzeroOrA<<
151 "tzeroOrC="<< tzeroOrC<<
152 "amplPMT1="<<meanAmplitude[0]<<
153 "amplPMT2="<<meanAmplitude[1]<<
154 "amplPMT3="<<meanAmplitude[2]<<
155 "amplPMT4="<<meanAmplitude[3]<<
156 "amplPMT5="<<meanAmplitude[4]<<
157 "amplPMT6="<<meanAmplitude[5]<<
158 "amplPMT7="<<meanAmplitude[6]<<
159 "amplPMT8="<<meanAmplitude[7]<<
160 "amplPMT9="<<meanAmplitude[8]<<
161 "amplPMT10="<<meanAmplitude[9]<<
162 "amplPMT11="<<meanAmplitude[10]<<
163 "amplPMT12="<<meanAmplitude[11]<<
164 "amplPMT13="<<meanAmplitude[12]<<
165 "amplPMT14="<<meanAmplitude[13]<<
166 "amplPMT15="<<meanAmplitude[14]<<
167 "amplPMT16="<<meanAmplitude[15]<<
168 "amplPMT17="<<meanAmplitude[16]<<
169 "amplPMT18="<<meanAmplitude[17]<<
170 "amplPMT19="<<meanAmplitude[18]<<
171 "amplPMT20="<<meanAmplitude[19]<<
172 "amplPMT21="<<meanAmplitude[20]<<
173 "amplPMT22="<<meanAmplitude[21]<<
174 "amplPMT23="<<meanAmplitude[22]<<
175 "amplPMT24="<<meanAmplitude[23]<<
176 "timePMT1="<<meanTime[0]<<
177 "timePMT2="<<meanTime[1]<<
178 "timePMT3="<<meanTime[2]<<
179 "timePMT4="<<meanTime[3]<<
180 "timePMT5="<<meanTime[4]<<
181 "timePMT6="<<meanTime[5]<<
182 "timePMT7="<<meanTime[6]<<
183 "timePMT8="<<meanTime[7]<<
184 "timePMT9="<<meanTime[8]<<
185 "timePMT10="<<meanTime[9]<<
186 "timePMT11="<<meanTime[10]<<
187 "timePMT12="<<meanTime[11]<<
188 "timePMT13="<<meanTime[12]<<
189 "timePMT14="<<meanTime[13]<<
190 "timePMT15="<<meanTime[14]<<
191 "timePMT16="<<meanTime[15]<<
192 "timePMT17="<<meanTime[16]<<
193 "timePMT18="<<meanTime[17]<<
194 "timePMT19="<<meanTime[18]<<
195 "timePMT20="<<meanTime[19]<<
196 "timePMT21="<<meanTime[20]<<
197 "timePMT22="<<meanTime[21]<<
198 "timePMT23="<<meanTime[22]<<
199 "timePMT24="<<meanTime[23];
200(*pcstream)<<"t0QA"<<
201 "timeDelayPMT1="<<timeDelayOCDB[0]<<
202 "timeDelayPMT2="<<timeDelayOCDB[1]<<
203 "timeDelayPMT3="<<timeDelayOCDB[2]<<
204 "timeDelayPMT4="<<timeDelayOCDB[3]<<
205 "timeDelayPMT5="<<timeDelayOCDB[4]<<
206 "timeDelayPMT6="<<timeDelayOCDB[5]<<
207 "timeDelayPMT7="<<timeDelayOCDB[6]<<
208 "timeDelayPMT8="<<timeDelayOCDB[7]<<
209 "timeDelayPMT9="<<timeDelayOCDB[8]<<
210 "timeDelayPMT10="<<timeDelayOCDB[9]<<
211 "timeDelayPMT11="<<timeDelayOCDB[10]<<
212 "timeDelayPMT12="<<timeDelayOCDB[11]<<
213 "timeDelayPMT13="<<timeDelayOCDB[12]<<
214 "timeDelayPMT14="<<timeDelayOCDB[13]<<
215 "timeDelayPMT15="<<timeDelayOCDB[14]<<
216 "timeDelayPMT16="<<timeDelayOCDB[15]<<
217 "timeDelayPMT17="<<timeDelayOCDB[16]<<
218 "timeDelayPMT18="<<timeDelayOCDB[17]<<
219 "timeDelayPMT19="<<timeDelayOCDB[18]<<
220 "timeDelayPMT20="<<timeDelayOCDB[19]<<
221 "timeDelayPMT21="<<timeDelayOCDB[20]<<
222 "timeDelayPMT22="<<timeDelayOCDB[21]<<
223 "timeDelayPMT23="<<timeDelayOCDB[22]<<
224 "timeDelayPMT24="<<timeDelayOCDB[23];
225
226 //-----> add the mean/sigma of the new histogram here
227
228 (*pcstream)<<"t0QA"<<"\n";
229
230 pcstream->Close();
231
232 delete pcstream;
233
234 return;
235
236}
237
238//_____________________________________________________________________________
239double GetParameterGaus(TH1F *histo, int whichParameter){
240
241 int maxBin = histo->GetMaximumBin();
242 double max = (histo->GetBinContent(maxBin-1) + histo->GetBinContent(maxBin) + histo->GetBinContent(maxBin+1))/3;
243 double mean = histo->GetBinCenter(maxBin); //mean
244 double lowfwhm = histo->GetBinCenter(histo->FindFirstBinAbove(max/2));
245 double highfwhm = histo->GetBinCenter(histo->FindLastBinAbove(max/2));
246 double sigma = (highfwhm - lowfwhm)/2.35482; //estimate fwhm FWHM = 2.35482*sigma
247
248 TF1 *gaussfit = new TF1("gaussfit","gaus", mean - 4*sigma, mean + 4*sigma); // fit in +- 4 sigma window
249 gaussfit->SetParameters(max, mean, sigma);
250
251 if(whichParameter==2)
252 histo->Fit(gaussfit,"RQNI");
253 else
254 histo->Fit(gaussfit,"RQN");
255
256 double parValue = gaussfit->GetParameter(whichParameter);
257
258 delete gaussfit;
259
260 return parValue;
261}
262
263