Flexible pt range for the efficiency histogramming
[u/mrichter/AliRoot.git] / T0 / MakeTrendT0.C
CommitLineData
b1a2430c 1#define NPMTs 24
2
3int MakeTrendT0( char *infile, int run) {
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();
105 man->SetDefaultStorage("raw://");
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//_____________________________________________________________________________
236double 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