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