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