Non-implemented method are commented out or moved to the private part of the class
[u/mrichter/AliRoot.git] / T0 / MakeTrendT0.C
1 #define NPMTs 24
2
3 int MakeTrendT0( char *infile, int run, char* ocdbStorage="raw://") {
4
5   gSystem->Load("libANALYSIS");
6   gSystem->Load("libANALYSISalice");
7   gSystem->Load("libCORRFW");
8   gSystem->Load("libTENDER");
9   gSystem->Load("libPWGPP.so");
10
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(ocdbStorage);
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