New DAQ calibration DAs by Raphaelle
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraVector.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 /////////////////////////////////////////////////////////////////////////////////
19 //                                                                             
20 // AliTRDCalibraVector                                                             
21 //                                                                             
22 // This class is for the vector methode of the TRD calibration.        
23 //                            
24 // Author:
25 //   R. Bailhache (R.Bailhache@gsi.de)
26 //                            
27 //////////////////////////////////////////////////////////////////////////////////////
28
29 #include <TGraphErrors.h>
30 #include <TH1F.h>
31 #include <TTree.h>
32 #include <TObjArray.h>
33 #include <TMath.h>
34 #include <TChain.h>
35 #include <TDirectory.h>
36 #include <TROOT.h>
37 #include <TFile.h>
38
39 #include "AliLog.h"
40
41 #include "AliTRDCalibraVector.h"
42 #include "AliTRDCommonParam.h"
43
44 ClassImp(AliTRDCalibraVector)
45
46 //______________________________________________________________________________________
47 AliTRDCalibraVector::AliTRDCalibraVector()
48   :TObject()
49   ,fVectorPH(new TObjArray())
50   ,fPlaPH(new TObjArray())
51   ,fVectorCH(new TObjArray())
52   ,fPlaCH(new TObjArray())
53   ,fVectorPRF(new TObjArray())
54   ,fPlaPRF(new TObjArray())
55   ,fNumberBinCharge(0)
56   ,fNumberBinPRF(0)
57   ,fTimeMax(0)
58   ,fPRFRange(1.5)
59 {
60   //
61   // Default constructor
62   //
63    
64 }
65
66 //______________________________________________________________________________________
67 AliTRDCalibraVector::AliTRDCalibraVector(const AliTRDCalibraVector &c)
68   :TObject(c)
69   ,fVectorPH(new TObjArray())
70   ,fPlaPH(new TObjArray())
71   ,fVectorCH(new TObjArray())
72   ,fPlaCH(new TObjArray())
73   ,fVectorPRF(new TObjArray())
74   ,fPlaPRF(new TObjArray())
75   ,fNumberBinCharge(c.fNumberBinCharge)
76   ,fNumberBinPRF(c.fNumberBinPRF)
77   ,fTimeMax(c.fTimeMax)
78   ,fPRFRange(c.fPRFRange)
79 {
80   //
81   // Copy constructor
82   //
83   for(Int_t k = 0; k < c.fVectorPH->GetEntriesFast(); k++){
84     AliTRDPInfo *obj = new AliTRDPInfo();
85     UShort_t  *entries    = new UShort_t[fTimeMax];
86     Float_t   *sum        = new Float_t[fTimeMax];
87     Float_t   *sumsquare  = new Float_t[fTimeMax];
88     for(Int_t i = 0; i < fTimeMax; i++){
89       entries[i]    = ((AliTRDPInfo *)c.fVectorPH->UncheckedAt(k))->GetEntries()[i];
90       sum[i]        = ((AliTRDPInfo *)c.fVectorPH->UncheckedAt(k))->GetSum()[i];
91       sumsquare[i]  = ((AliTRDPInfo *)c.fVectorPH->UncheckedAt(k))->GetSumSquare()[i];
92     }
93     obj->SetEntries(entries);
94     obj->SetSum(sum);
95     obj->SetSumSquare(sumsquare);
96     fVectorPH->AddAt((TObject *)obj,k);
97   } 
98   for(Int_t k = 0; k < c.fVectorPRF->GetEntriesFast(); k++){
99     AliTRDPInfo *obj = new AliTRDPInfo();
100     UShort_t  *entries    = new UShort_t[fNumberBinPRF];
101     Float_t   *sum        = new Float_t[fNumberBinPRF];
102     Float_t   *sumsquare  = new Float_t[fNumberBinPRF];
103     for(Int_t i = 0; i < fNumberBinPRF; i++){
104       entries[i]    = ((AliTRDPInfo *)c.fVectorPRF->UncheckedAt(k))->GetEntries()[i];
105       sum[i]        = ((AliTRDPInfo *)c.fVectorPRF->UncheckedAt(k))->GetSum()[i];
106       sumsquare[i]  = ((AliTRDPInfo *)c.fVectorPRF->UncheckedAt(k))->GetSumSquare()[i];
107     }
108     obj->SetEntries(entries);
109     obj->SetSum(sum);
110     obj->SetSumSquare(sumsquare);   
111     fVectorPRF->AddAt((TObject *)obj,k);
112   } 
113   for(Int_t k = 0; k < c.fVectorCH->GetEntriesFast(); k++){
114     AliTRDCTInfo *obj = new AliTRDCTInfo();
115     UShort_t  *entries = new UShort_t[fNumberBinCharge];
116     for(Int_t i = 0; i < fNumberBinCharge; i++){
117       entries[i] = ((AliTRDCTInfo *)c.fVectorCH->UncheckedAt(k))->GetEntries()[i];
118     }
119     obj->SetEntries(entries);
120     fVectorCH->AddAt((TObject *)obj,k);
121   } 
122   for(Int_t k = 0; k < c.fPlaPH->GetEntriesFast(); k++){
123     AliTRDPlace *obj = new AliTRDPlace();
124     obj->SetPlace(((AliTRDPlace *)c.fPlaPH->UncheckedAt(k))->GetPlace());
125     fPlaPH->AddAt((TObject *)obj,k);
126   } 
127   for(Int_t k = 0; k < c.fPlaCH->GetEntriesFast(); k++){
128     AliTRDPlace *obj = new AliTRDPlace();
129     obj->SetPlace(((AliTRDPlace *)c.fPlaCH->UncheckedAt(k))->GetPlace());
130     fPlaCH->AddAt((TObject *)obj,k);
131   } 
132   for(Int_t k = 0; k < c.fPlaPRF->GetEntriesFast(); k++){
133     AliTRDPlace *obj = new AliTRDPlace();
134     obj->SetPlace(((AliTRDPlace *)c.fPlaPRF->UncheckedAt(k))->GetPlace());
135     fPlaPRF->AddAt((TObject *)obj,k);
136   } 
137
138 }
139 //____________________________________________________________________________________
140 AliTRDCalibraVector::~AliTRDCalibraVector()
141 {
142   //
143   // AliTRDCalibraVector destructor
144   //
145
146 }
147
148 //_____________________________________________________________________________
149 Int_t AliTRDCalibraVector::SearchBin(Float_t value, Int_t i) const
150 {
151   //
152   // Search the bin
153   //
154
155   Int_t reponse        = 0;
156   Float_t fbinmin      = 0;
157   Float_t fbinmax      = value;
158   Int_t fNumberOfBin   = -1;
159
160   // Charge
161   if (i == 0) {
162     fbinmax      = 300.0;
163     fbinmin      = 0.0;
164     fNumberOfBin = fNumberBinCharge;
165   }
166
167   // PRF
168   if (i == 2) {
169     fbinmax      =   TMath::Abs(fPRFRange);
170     fbinmin      =  -TMath::Abs(fPRFRange);
171     fNumberOfBin = fNumberBinPRF;
172   }
173
174   // Return -1 if out
175   if ((value >= fbinmax) || 
176       (value <  fbinmin)) {
177     return -1;
178   }
179   // Sinon
180   else {
181     reponse = (Int_t) ((fNumberOfBin*(value-fbinmin)) / (fbinmax-fbinmin));
182   }
183
184   return reponse;
185
186 }
187
188 //_____________________________________________________________________________
189 Int_t AliTRDCalibraVector::SearchInVector(Int_t group, Int_t i) const
190 {
191   //
192   // Search if the calibration group "group" has already been
193   // initialised by a previous track in the vector
194   //
195
196   if (i == 0) {
197     for (Int_t k = 0; k < (Int_t) fPlaCH->GetEntriesFast(); k++) {
198       if (((AliTRDPlace *) fPlaCH->At(k))->GetPlace() == group) {
199         return k;
200       }
201     }
202     return -1;
203   }
204
205   if (i == 1) {
206     for (Int_t k = 0; k < (Int_t) fPlaPH->GetEntriesFast(); k++) {
207       if (((AliTRDPlace *) fPlaPH->At(k))->GetPlace() == group) {
208         return k;
209       }
210     }
211     return -1;
212   }
213
214   if (i == 2) {
215     for (Int_t k = 0; k < (Int_t) fPlaPRF->GetEntriesFast(); k++) {
216       if (((AliTRDPlace *) fPlaPRF->At(k))->GetPlace() == group) {
217         return k;
218       }
219     }
220     return -1;
221   }
222
223   return -1;
224
225 }
226
227 //_____________________________________________________________________________
228 Int_t AliTRDCalibraVector::SearchInTreeVector(TObjArray *vectorplace, Int_t group) const
229 {
230   //
231   // Search if the calibration group "group" is present in the tree
232   //
233
234   for (Int_t k = 0; k < (Int_t) vectorplace->GetEntriesFast(); k++) {
235     if (((AliTRDPlace *) vectorplace->At(k))->GetPlace() == group) {
236       return k;
237     }
238   }
239
240   return -1;
241
242 }
243
244 //_____________________________________________________________________________
245 Bool_t AliTRDCalibraVector::UpdateVectorCH(Int_t group, Float_t value)
246 {
247   //
248   // Fill the vector if a new calibration group "group" or update the
249   // values of the calibration group "group" if already here  
250   //
251
252   // Search bin
253   Int_t bin = SearchBin(value,0);
254   // Out
255   if ((bin < 0) || (bin >= fNumberBinCharge)) {
256     return kFALSE; 
257   }
258
259   // Search place
260   Int_t place = SearchInVector(group,0);
261
262   // New group
263   if (place == -1) {
264     AliTRDPlace *placegroup = new AliTRDPlace();
265     placegroup->SetPlace(group);
266     fPlaCH->Add((TObject *) placegroup);
267     // Variable
268     AliTRDCTInfo *fCHInfo = new AliTRDCTInfo();
269     UShort_t *entries = new UShort_t[fNumberBinCharge];
270     // Initialise first
271     for(Int_t k = 0; k < fNumberBinCharge; k++) {
272       entries[k] = 0;
273     }
274     // Add the value
275     entries[bin]= 1;
276     // Set
277     fCHInfo->SetEntries(entries);
278     // Set in the vector
279     fVectorCH->Add((TObject *) fCHInfo);
280   }
281   // Group already exits
282   else {
283     // Variable
284     AliTRDCTInfo *fCHInfo = new AliTRDCTInfo();
285     // Retrieve
286     fCHInfo = ((AliTRDCTInfo *) fVectorCH->At(place));
287     UShort_t *entries = fCHInfo->GetEntries();
288     // Add
289     entries[bin]++;
290     // Set
291     fCHInfo->SetEntries(entries);
292     // Update the vector
293     fVectorCH->AddAt((TObject *) fCHInfo,place);
294   }
295
296   return kTRUE;
297
298 }
299
300 //_____________________________________________________________________________
301 Bool_t AliTRDCalibraVector::UpdateVectorPRF(Int_t group, Float_t x, Float_t y)
302 {
303   //
304   // Fill the vector if a new calibration group "group" or update the
305   // values of the calibration group "group" if already here  
306   //
307
308   // Search bin
309   Int_t bin = SearchBin(x,2);
310   // Out
311   if ((bin < 0) || (bin >= fNumberBinPRF)) {
312     return kFALSE; 
313   }
314
315   // Search place
316   Int_t place = SearchInVector(group,2);
317
318   // New group
319   if (place == -1) {
320
321     AliTRDPlace *placegroup = new AliTRDPlace();
322     placegroup->SetPlace(group);
323     fPlaPRF->Add((TObject *) placegroup);
324     AliTRDPInfo *fPRFInfo = new AliTRDPInfo();
325
326     Float_t  *sum       = new Float_t[fNumberBinPRF];
327     Float_t  *sumsquare = new Float_t[fNumberBinPRF];
328     UShort_t *entries   = new UShort_t[fNumberBinPRF];
329
330     // Initialise first
331     for (Int_t k = 0; k < fNumberBinPRF; k++) {
332       sum[k]       = 0.0;
333       sumsquare[k] = 0.0;
334       entries[k]   = 0;
335     }
336
337     // Add the value
338     sum[bin]       += y;
339     sumsquare[bin] += y*y;
340     entries[bin]++;
341
342     // Set
343     fPRFInfo->SetSum(sum);
344     fPRFInfo->SetSumSquare(sumsquare);
345     fPRFInfo->SetEntries(entries);
346
347     // Set in the vector
348     fVectorPRF->Add((TObject *) fPRFInfo);
349         
350   }
351   // Group already exits
352   else {
353
354     AliTRDPInfo *fPRFInfo = new AliTRDPInfo();
355     // Retrieve
356     fPRFInfo = (AliTRDPInfo *) fVectorPRF->At(place);
357
358     Float_t  *sum       = fPRFInfo->GetSum();
359     Float_t  *sumsquare = fPRFInfo->GetSumSquare();
360     UShort_t *entries   = fPRFInfo->GetEntries();
361
362     // Add
363     Double_t calcul       = (((Double_t) fPRFInfo->GetEntries()[bin])
364                            * ((Double_t) fPRFInfo->GetSum()[bin]) + (Double_t) y)
365                           / (((Double_t) fPRFInfo->GetEntries()[bin]) + 1);
366     sum[bin]       = (Float_t) calcul;
367     Double_t calculsquare = (((Double_t) fPRFInfo->GetSumSquare()[bin])
368                            * ((Double_t) fPRFInfo->GetEntries()[bin]) + ((Double_t) y)*((Double_t) y))
369                           / (((Double_t) fPRFInfo->GetEntries()[bin]) + 1);
370     sumsquare[bin] = (Float_t) calculsquare;
371     entries[bin]++;
372
373     // Set
374     fPRFInfo->SetSum(sum);
375     fPRFInfo->SetSumSquare(sumsquare);
376     fPRFInfo->SetEntries(entries);
377  
378     // Update the vector
379     fVectorPRF->AddAt((TObject *) fPRFInfo,place);
380
381   }
382
383   return kTRUE;
384
385 }
386
387 //_____________________________________________________________________________
388 Bool_t AliTRDCalibraVector::UpdateVectorPH(Int_t group, Int_t time, Float_t value)
389 {
390   //
391   // Fill the vector if a new calibration group "group" or update
392   // the values of the calibration group "group" if already here  
393   //
394
395   // Search bin
396   Int_t bin = time;
397   // Out
398   if ((bin <         0) || 
399       (bin >= fTimeMax)) {
400     return kFALSE; 
401   }
402
403   // Search place
404   Int_t place = SearchInVector(group,1);
405
406   // New group
407   if(place == -1){
408
409     AliTRDPlace *placegroup = new AliTRDPlace();
410     placegroup->SetPlace(group);
411     fPlaPH->Add((TObject *) placegroup);
412     AliTRDPInfo *fPHInfo = new AliTRDPInfo();
413
414     Float_t  *sum       = new Float_t[fTimeMax];
415     Float_t  *sumsquare = new Float_t[fTimeMax];
416     UShort_t *entries   = new UShort_t[fTimeMax];
417
418     // Initialise first
419     for (Int_t k = 0; k < fTimeMax; k++) {
420       sum[k]       = 0.0;
421       sumsquare[k] = 0.0;
422       entries[k]   = 0;
423     }
424
425     // Add the value
426     sum[bin]       += value;
427     sumsquare[bin] += value*value;
428     entries[bin]++;
429
430     // Set
431     fPHInfo->SetSum(sum);
432     fPHInfo->SetSumSquare(sumsquare);
433     fPHInfo->SetEntries(entries);
434
435     // Set in the vector
436     fVectorPH->Add((TObject *) fPHInfo);
437
438   }
439   // Group already exits
440   else {
441
442     AliTRDPInfo *fPHInfo = new AliTRDPInfo();
443     // Retrieve
444     fPHInfo = (AliTRDPInfo *) fVectorPH->At(place);
445
446     Float_t  *sum       = fPHInfo->GetSum();
447     Float_t  *sumsquare = fPHInfo->GetSumSquare();
448     UShort_t *entries   = fPHInfo->GetEntries();
449
450     // Add
451     Double_t calcul       = (((Double_t) fPHInfo->GetEntries()[bin])
452                            * ((Double_t) fPHInfo->GetSum()[bin]) + (Double_t) value)
453                           / (((Double_t) fPHInfo->GetEntries()[bin]) + 1);
454     sum[bin]       = (Float_t) calcul;
455     Double_t calculsquare = ((((Double_t) fPHInfo->GetSumSquare()[bin])
456                             * ((Double_t) fPHInfo->GetEntries()[bin])) 
457                           + (((Double_t) value) * ((Double_t)value))) 
458                           / (((Double_t) fPHInfo->GetEntries()[bin]) + 1);
459     sumsquare[bin] = (Float_t) calculsquare;
460     entries[bin]++;
461
462     // Set
463     fPHInfo->SetSum(sum);
464     fPHInfo->SetSumSquare(sumsquare);
465     fPHInfo->SetEntries(entries);
466
467     // Update the vector
468     fVectorPH->AddAt((TObject *) fPHInfo,place);
469
470   }
471
472   return kTRUE;
473
474 }  
475
476 //_____________________________________________________________________________
477 TGraphErrors *AliTRDCalibraVector::ConvertVectorPHisto(Int_t place
478                                         , const Char_t * name) const
479 {
480   //
481   // Convert the PInfo in a TGraphErrors
482   //
483
484   AliTRDPInfo *pInfo = new AliTRDPInfo();
485   // Retrieve
486   pInfo = ((AliTRDPInfo *) fVectorPH->At(place));
487   TGraphErrors *histo;
488   histo = ConvertVectorPHistoI((AliTRDPInfo *)pInfo, name);
489   
490   return histo;
491
492
493
494 //_____________________________________________________________________________
495 TGraphErrors *AliTRDCalibraVector::ConvertVectorPHistoI(AliTRDPInfo *pInfo
496                                                , const Char_t *name) const
497 {
498   //
499   // Convert the PInfo in a 1D grapherror, name must contains "PRF"
500   // if PRF calibration and not "PRF" for Vdrift calibration
501   //
502  
503   TGraphErrors *histo;
504   const Char_t *pattern1 = "PRF";
505
506   Float_t sf = 10.0;
507   AliTRDCommonParam *parCom = AliTRDCommonParam::Instance();
508   if (!parCom) {
509     AliInfo("Could not get CommonParam, take the default 10MHz");
510   }
511   
512   sf = parCom->GetSamplingFrequency();
513
514   // Axis
515   Double_t *x;
516   Double_t *y;
517   Double_t *ex;
518   Double_t *ey;
519   Double_t step = 0.0;
520   Double_t min  = 0.0;
521
522   // Ntimes
523   Int_t ntimes = 0;
524   if (strstr(name,pattern1)) {
525     ntimes = fNumberBinPRF;
526   }
527   else {
528     ntimes = fTimeMax;
529   }
530   x  = new Double_t[ntimes]; // Xaxis
531   y  = new Double_t[ntimes]; // Sum
532   ex = new Double_t[ntimes]; // Nentries
533   ey = new Double_t[ntimes]; // Sum of square/nentries
534
535   // Init histo
536   if (!strstr(name,pattern1)) {
537     step = 1.0 / sf;
538     min  = 0.0;
539   }
540   else {
541     step = (2*TMath::Abs(fPRFRange)) / fNumberBinPRF;
542     min  = -TMath::Abs(fPRFRange) + step / 2.0;
543   }
544
545   // Fill histo
546   for (Int_t k = 0; k < ntimes; k++) {
547     x[k]  = min + k*step;
548     y[k]  = 0.0;
549     ex[k] = 0.0;
550     ey[k] = 0.0;
551     // Fill only if there is more than 0 something
552     if (pInfo->GetEntries()[k] > 0) {
553       ex[k] = pInfo->GetEntries()[k];
554       y[k]  = pInfo->GetSum()[k];
555       ey[k] =  pInfo->GetSumSquare()[k];
556     }
557   }
558
559   // Define the TGraphErrors
560   histo = new TGraphErrors(ntimes,x,y,ex,ey);
561   histo->SetTitle(name); 
562   return histo;
563
564 }
565
566 //_____________________________________________________________________________
567 TH1F *AliTRDCalibraVector::ConvertVectorCTHisto(Int_t place
568                                               , const Char_t * name) const
569 {
570   //
571   // Convert the CTInfo in a 1D histo
572   //
573
574   AliTRDCTInfo *cHInfo = new AliTRDCTInfo();
575   // Retrieve
576   cHInfo = ((AliTRDCTInfo *) fVectorCH->At(place));
577   TH1F *histo;
578   histo = ConvertVectorCTHistoI((AliTRDCTInfo *)cHInfo,(const Char_t *) name);
579   
580   return histo;
581
582
583
584 //_____________________________________________________________________________
585 TH1F *AliTRDCalibraVector::ConvertVectorCTHistoI(AliTRDCTInfo *cTInfo
586                                                , const Char_t * name) const
587 {
588   //
589   // Convert the CTInfo in a 1D histo
590   //
591
592   TH1F *histo;
593   
594   Int_t     ntimes  = fNumberBinCharge;
595   UShort_t *entries = cTInfo->GetEntries();
596   
597   // Init histo
598   histo = new TH1F(name,name,fNumberBinCharge,0,300);
599   histo->Sumw2();
600   // Fill histo
601   for (Int_t k = 0; k < ntimes; k++) {
602     histo->SetBinContent(k+1,entries[k]);
603     histo->SetBinError(k+1,TMath::Sqrt(TMath::Abs(entries[k])));
604   }
605   
606   return histo;
607
608 }
609
610 //_____________________________________________________________________________
611 TTree *AliTRDCalibraVector::ConvertVectorCTTreeHisto(TObjArray *vVectorCT
612                                                    , TObjArray *pPlaCT
613                                                    , const Char_t *name
614                                                    , const Char_t *nametitle) const
615 {
616   //
617   // Convert the vector in a tree with two branchs: the group number
618   // and the TH1F histo reconstructed from the vector
619   //
620
621   // Size of the things
622   Int_t ntotal = (Int_t) pPlaCT->GetEntriesFast();
623   if (ntotal == 0) {
624     AliInfo("nothing to write!");
625     TTree *treeCT = new TTree(name,nametitle);
626     return treeCT;
627   }
628   
629   // Variable of the tree
630   Int_t groupnumber = -1; // Group calibration
631   TH1F      *histo = 0x0;
632   TObjArray  vectorCT = *vVectorCT;
633   TObjArray  plaCT    = *pPlaCT;
634
635   // Init the tree
636   TTree *treeCT = new TTree(name,nametitle);
637   treeCT->Branch("groupnumber",&groupnumber,"groupnumber/I");
638   treeCT->Branch("histo","TH1F",&histo,32000,0);
639
640   // Fill
641   Int_t k = 0;
642   while (k < ntotal) {
643     TString nome(name);
644     groupnumber  = ((AliTRDPlace *) plaCT.At(0))->GetPlace();
645     nome        += groupnumber;
646     histo        = ConvertVectorCTHistoI(((AliTRDCTInfo *) vectorCT.At(0)),nome);
647     treeCT->Fill();
648     vectorCT.RemoveAt(0);
649     vectorCT.Compress();
650     plaCT.RemoveAt(0);
651     plaCT.Compress();
652     k++;
653   } 
654
655   return treeCT;
656
657 }
658
659 //_____________________________________________________________________________
660 TTree *AliTRDCalibraVector::ConvertVectorPTreeHisto(TObjArray *vVectorP
661                                                   , TObjArray *pPlaP
662                                                   , const Char_t *name
663                                                   , const Char_t *nametitle) const
664 {
665   //
666   // Convert the vector in a tree with two branchs: the group number
667   // and the TGraphErrors histo reconstructed from the vector.
668   // The name must contain "PRF" for PRF calibration and not "PRF"
669   // for Vdrift calibration
670   //
671
672   // Size of the things
673   Int_t ntotal = (Int_t) pPlaP->GetEntriesFast();
674   if (ntotal == 0) {
675     AliInfo("nothing to write!");
676     TTree *treeP = new TTree(name,nametitle);
677     return treeP;
678   }
679
680   // Variable of the tree
681   Int_t groupnumber = -1; // Group calibration
682   TGraphErrors *histo   = 0x0;
683   TObjArray     vectorP = *vVectorP;
684   TObjArray     plaP    = *pPlaP;
685
686   // Init the tree
687   TTree *treeP = new TTree(name,nametitle);
688   treeP->Branch("groupnumber",&groupnumber,"groupnumber/I");
689   treeP->Branch("histo","TGraphErrors",&histo,32000,0);
690
691   // Fill
692   Int_t k = 0;
693   while (k < ntotal) {
694     TString nome(name);
695     groupnumber = ((AliTRDPlace *) plaP.At(0))->GetPlace();
696     nome       += groupnumber;
697     histo       = ConvertVectorPHistoI((AliTRDPInfo *) vectorP.At(0),nome);
698     treeP->Fill();
699     vectorP.RemoveAt(0);
700     vectorP.Compress();
701     plaP.RemoveAt(0);
702     plaP.Compress();
703     k++;
704   } 
705
706   return treeP;
707
708 }
709
710 //_____________________________________________________________________________
711 TObjArray *AliTRDCalibraVector::ConvertTreeVector(TTree *tree) const
712 {
713   //
714   // Convert the branch groupnumber of the tree taken from
715   // TRD.calibration.root in case of vector method in a std::vector 
716   // to be faster
717   //
718
719   // Initialise
720   TObjArray *vectorplace = new TObjArray();
721   
722   // Variable of the tree
723   Int_t groupnumber = -1; // Group calibration
724
725   // Set the branch
726   tree->SetBranchAddress("groupnumber",&groupnumber);
727     
728   // Fill
729   Int_t ntotal = tree->GetEntries();
730   for (Int_t k = 0; k < ntotal; k++) {
731     tree->GetEntry(k);
732     AliTRDPlace *placegroupnumber = new AliTRDPlace();
733     placegroupnumber->SetPlace(groupnumber);
734     vectorplace->Add((TObject *) placegroupnumber);
735   }
736   
737   return vectorplace;
738
739 }
740
741 //_____________________________________________________________________________
742 Bool_t AliTRDCalibraVector::MergeVectorCT(TObjArray *vVectorCT2, TObjArray *pPlaCT2)
743 {
744   //
745   // Add the two vectors and place the result in the first
746   //
747
748   if (((Int_t) pPlaCT2->GetEntriesFast()) != ((Int_t) vVectorCT2->GetEntriesFast())){
749     AliInfo("VectorCT2 doesn't correspond to PlaCT2!");
750     return kFALSE;
751   }
752  
753   // CH case
754   for (Int_t k = 0; k < (Int_t) fPlaCH->GetEntriesFast(); k++) {
755     
756     // Look if PlaCT1[k] it is also in the second vector
757     Int_t place = -1;
758     for (Int_t j = 0; j < (Int_t) pPlaCT2->GetEntriesFast(); j++) {
759       if (((AliTRDPlace *) pPlaCT2->At(j))->GetPlace() == 
760             ((AliTRDPlace *) fPlaCH->At(k))->GetPlace()) {
761         place = j;
762         break;
763       }
764     }
765     
766     // If not in the second vector nothing to do
767
768     // If in the second vector
769     if (place != -1) {
770       
771       AliTRDCTInfo *fCTInfo = new AliTRDCTInfo();
772       UShort_t *entries = new UShort_t[fNumberBinCharge];
773       
774       for (Int_t nu = 0; nu < fNumberBinCharge; nu++) {
775         entries[nu] = ((AliTRDCTInfo *)  fVectorCH->At(((AliTRDPlace *) fPlaCH->At(k))->GetPlace()))->GetEntries()[nu]
776                     + ((AliTRDCTInfo *) vVectorCT2->At(((AliTRDPlace *) fPlaCH->At(k))->GetPlace()))->GetEntries()[nu];
777       }
778       
779       // Set
780       fCTInfo->SetEntries(entries);
781
782       // Nothing to do on PlaCT1
783       
784       // Update the vector 
785       fVectorCH->AddAt((TObject *) fCTInfo,((AliTRDPlace *) fPlaCH->At(k))->GetPlace());
786
787     }
788     
789   } 
790  
791   // And at the end the vector in CT2 but not in CH1
792   for (Int_t k = 0; k < (Int_t) pPlaCT2->GetEntriesFast(); k++) {
793     
794     // Look if pPlaCT2[k] it is also in the second vector
795     Int_t place = -1;
796     for (Int_t j = 0; j < (Int_t) fPlaCH->GetEntriesFast(); j++) {
797       if (((AliTRDPlace *) fPlaCH->At(j))->GetPlace() == ((AliTRDPlace *) pPlaCT2->At(k))->GetPlace()) {
798         place = j;
799         break;
800       }
801     }
802
803     // If not in the first vector
804     if (place == -1) {
805       
806       AliTRDCTInfo *fCTInfo = new AliTRDCTInfo();     
807       fCTInfo = ((AliTRDCTInfo *) vVectorCT2->At(((AliTRDPlace *) pPlaCT2->At(k))->GetPlace()));
808       
809       // Add at the end 
810       fPlaCH->Add((TObject *) (pPlaCT2->At(k)));
811       fVectorCH->Add((TObject *) fCTInfo);
812
813     }
814     
815   }
816   
817   return kTRUE;
818   
819 }
820
821 //_____________________________________________________________________________
822 Bool_t AliTRDCalibraVector::MergeVectorP(TObjArray *vVectorP2
823                                        , TObjArray *pPlaP2
824                                        , Int_t i)
825 {
826   //
827   // Add the two vectors and place the result in the first
828   //
829
830   if (((Int_t) pPlaP2->GetEntriesFast()) != ((Int_t) vVectorP2->GetEntriesFast())) {
831     AliInfo("VectorP2 doesn't correspond to PlaP2!");
832     return kFALSE;
833   }
834
835   // PH case
836   if (i == 1) {
837
838      for (Int_t k = 0; k < (Int_t) fPlaPH->GetEntriesFast(); k++) {
839        
840        // Look if fPlaPH[k] it is also in the second vector
841        Int_t place = -1;
842        for (Int_t j = 0; j < (Int_t) pPlaP2->GetEntriesFast(); j++) {
843          if (((AliTRDPlace *) pPlaP2->At(j))->GetPlace() == ((AliTRDPlace *) fPlaPH->At(k))->GetPlace()) {
844            place = j;
845            break;
846          }
847        }
848        
849        // If not in the second vector nothing to do
850
851        // If in the second vector
852        if (place != -1) {
853
854          AliTRDPInfo *fPInfo = new AliTRDPInfo();
855          UShort_t *entries   = new UShort_t[fTimeMax];
856          Float_t  *sum       = new Float_t[fTimeMax];
857          Float_t  *sumsquare = new Float_t[fTimeMax];
858
859          for (Int_t nu = 0; nu < fTimeMax; nu++) {
860            
861            entries[nu]   = ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu]
862                          + ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu];
863            
864            Double_t calcul       = ((((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSum()[nu])
865                                    * ((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu]))
866                                   + (((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSum()[nu])
867                                    * ((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu])))
868                                  / ((Double_t) fPInfo->GetEntries()[nu]);
869
870            sum[nu]       = (Float_t) calcul;
871            
872            Double_t calculsquare = ((((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSumSquare()[nu])
873                                    * ((Double_t) ((AliTRDPInfo *) fVectorPH->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu]))
874                                   + (((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetSumSquare()[nu])
875                                    * ((Double_t) ((AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) fPlaPH->At(k))->GetPlace()))->GetEntries()[nu])))
876                                  / ((Double_t) fPInfo->GetEntries()[nu]);
877
878            sumsquare[nu] = calculsquare;
879
880          }
881
882          // Set
883          fPInfo->SetSum(sum);
884          fPInfo->SetSumSquare(sumsquare);
885          fPInfo->SetEntries(entries);
886          
887          // Nothing to do on PlaCT1
888          
889          // Update the vector VectorCT1
890          fVectorPH->AddAt((TObject *) fPInfo,((AliTRDPlace *) fPlaPH->At(k))->GetPlace());
891          
892        }
893
894      }
895
896      // And at the end the vector in P2 but not in CH1
897      for (Int_t k = 0; k < (Int_t) pPlaP2->GetEntriesFast(); k++) {
898        
899        // Look if PlaCT2[k] it is also in the second vector
900        Int_t place = -1;
901        for (Int_t j = 0; j < (Int_t) fPlaPH->GetEntriesFast(); j++) {
902          if (((AliTRDPlace *) fPlaPH->At(j))->GetPlace() == ((AliTRDPlace *) pPlaP2->At(k))->GetPlace()) {
903            place = j;
904            break;
905          }
906        }
907        
908        // If not in the first vector
909        if (place == -1) {
910                  
911          AliTRDPInfo *fPInfo = new AliTRDPInfo();
912          fPInfo = (AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) pPlaP2->At(k))->GetPlace());
913          
914          // Add at the end of CH1
915          fPlaPH->Add(((TObject *) pPlaP2->At(k)));
916          fVectorPH->Add((TObject *) fPInfo);
917
918        }
919
920      }
921
922    }
923
924    // PRF case
925    if (i == 1) {
926
927      for (Int_t k = 0; k < (Int_t) fPlaPRF->GetEntriesFast(); k++) {
928
929        // Look if fPlaPRF[k] it is also in the second vector
930        Int_t place = -1;
931        for (Int_t j = 0; j < (Int_t) pPlaP2->GetEntriesFast(); j++) {
932          if (((AliTRDPlace *) pPlaP2->At(j))->GetPlace() == ((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()) {
933            place = j;
934            break;
935          }
936        }
937
938        // If not in the second vector nothing to do
939
940        // If in the second vector
941        if (place != -1) {
942         
943          AliTRDPInfo *fPInfo = new AliTRDPInfo();
944          UShort_t *entries   = new UShort_t[fNumberBinPRF];
945          Float_t  *sum       = new Float_t[fNumberBinPRF];
946          Float_t  *sumsquare = new Float_t[fNumberBinPRF];
947
948          for (Int_t nu = 0; nu < fNumberBinPRF; nu++) {
949            
950            entries[nu]           = ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu]
951                                  + ((AliTRDPInfo *)  vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu];
952            
953            Double_t calcul       = ((((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSum()[nu])
954                                    * ((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu]))
955                                   + (((Double_t) ((AliTRDPInfo *)  vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSum()[nu])
956                                    * ((Double_t) ((AliTRDPInfo *)  vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu])))
957                                  / ((Double_t) fPInfo->GetEntries()[nu]);
958
959            sum[nu]               = (Float_t) calcul;
960            
961            Double_t calculsquare = ((((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSumSquare()[nu])
962                                    * ((Double_t) ((AliTRDPInfo *) fVectorPRF->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu]))
963                                   + (((Double_t) ((AliTRDPInfo *)  vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetSumSquare()[nu])
964                                    * ((Double_t) ((AliTRDPInfo *)  vVectorP2->At(((AliTRDPlace *) fPlaPRF->At(k))->GetPlace()))->GetEntries()[nu])))
965                                  / ((Double_t) fPInfo->GetEntries()[nu]);
966
967            sumsquare[nu]         = calculsquare;
968
969          }
970
971          // Set
972          fPInfo->SetSum(sum);
973          fPInfo->SetSumSquare(sumsquare);
974          fPInfo->SetEntries(entries);
975
976          // Nothing to do on PlaCT1
977          
978          // Update the vector VectorCT1
979          fVectorPRF->AddAt((TObject *) fPInfo,((AliTRDPlace *) fPlaPRF->At(k))->GetPlace());
980          
981        }
982
983      }
984
985      // And at the end the vector in P2 but not in CH1
986      for (Int_t k = 0; k < (Int_t) pPlaP2->GetEntriesFast(); k++) {
987        
988        // Look if PlaCT2[k] it is also in the second vector
989        Int_t place = -1;
990        for (Int_t j = 0; j < (Int_t) fPlaPRF->GetEntriesFast(); j++) {
991          if (((AliTRDPlace *) fPlaPRF->At(j))->GetPlace() == ((AliTRDPlace *) pPlaP2->At(k))->GetPlace()) {
992            place = j;
993            break;
994          }
995        }
996
997        // If not in the first vector
998        if (place == -1) {
999
1000          AliTRDPInfo *fPInfo = new AliTRDPInfo();
1001          fPInfo = (AliTRDPInfo *) vVectorP2->At(((AliTRDPlace *) pPlaP2->At(k))->GetPlace());
1002
1003          // Add at the end of CH1
1004          fPlaPRF->Add(((TObject *) pPlaP2->At(k)));
1005          fVectorPRF->Add((TObject *) fPInfo);
1006
1007        }
1008        
1009      }
1010
1011    } 
1012
1013    return kTRUE;
1014
1015 }   
1016
1017 //_____________________________________________________________________________
1018 TTree *AliTRDCalibraVector::Sum2Trees(const Char_t *filename1
1019                                     , const Char_t *filename2
1020                                     , const Char_t *variablecali)
1021 {
1022   //
1023   // It returns the sum of two trees with the name variablecali
1024   // in the files filenam1 and filename2 equivalent of merging two 2D histos
1025   // The name of the resulting tree is the same as the two input trees
1026   // variablecali can be treeCH2d, treePH2d or treePRF2d 
1027   //
1028
1029   // Variables
1030   TChain    *treeChain   = new TChain(variablecali);
1031   TObjArray *vectorplace = new TObjArray();
1032   TObjArray *where       = new TObjArray();
1033   
1034   // First tree
1035   // Take the tree
1036   TFile *file1 = new TFile(filename1,"READ");
1037   TTree *tree1 = (TTree *) file1->Get(variablecali);
1038
1039   gDirectory = gROOT;
1040
1041   // Take the places
1042   vectorplace = ConvertTreeVector(tree1);
1043
1044   // Say where it is in tree 1
1045   for (Int_t jui = 0; jui < (Int_t) vectorplace->GetEntriesFast(); jui++) {
1046     AliTRDPlace *placejui = new AliTRDPlace();
1047     placejui->SetPlace(jui);
1048     TObjArray *chainplace = new TObjArray();
1049     chainplace->Add((TObject *) placejui);
1050     where->Add((TObject *) chainplace);
1051   }
1052
1053   // Add to the chain
1054   treeChain->Add(filename1);
1055   delete file1;
1056
1057   // Second tree
1058   // Take the tree
1059   TFile *file2 = new TFile(filename2,"READ");
1060   TTree *tree2 = (TTree *) file2->Get(variablecali);
1061
1062   gDirectory = gROOT;
1063
1064   // Take the places
1065   TObjArray *vector2 = ConvertTreeVector(tree2);
1066   Int_t j = treeChain->GetEntries();
1067
1068   for (Int_t jui = 0; jui < (Int_t) vector2->GetEntriesFast(); jui++) {
1069     // Search if already found
1070     Int_t place = SearchInTreeVector(vectorplace,((AliTRDPlace *) vector2->At(jui))->GetPlace());
1071     // Create a new element in the two std vectors
1072     if (place == -1) {
1073       AliTRDPlace *placejjui  = new AliTRDPlace();
1074       placejjui->SetPlace((j+jui));
1075       TObjArray   *chainplace = new TObjArray();
1076       chainplace->Add((TObject *) placejjui);
1077       vectorplace->Add((TObject *) (vector2->At(jui)));
1078       where->Add((TObject *) chainplace);
1079     }
1080     // Update the element at the place "place" in the std vector whereinthechain
1081     else {
1082       AliTRDPlace *placejjui  = new AliTRDPlace();
1083       placejjui->SetPlace((j+jui));
1084       TObjArray   *chainplace = ((TObjArray *) where->At(place));
1085       chainplace->Add((TObject *) placejjui);
1086       where->AddAt((TObject *) chainplace,place);
1087     }
1088   }
1089
1090   // Add to the Chain
1091   treeChain->Add(filename2);
1092   delete file2; 
1093
1094   // Take care of the profile
1095   const Char_t *pattern = "P";
1096   TTree *tree = 0x0;
1097
1098   if (!strstr(variablecali,pattern)) {
1099
1100     // Ready to read the chain
1101     TH1F *his = 0x0;
1102     treeChain->SetBranchAddress("histo",&his);
1103
1104     // Initialise the final tree
1105     Int_t group   = -1;
1106     TH1F *histsum = 0x0;
1107    
1108     tree = new TTree(variablecali,variablecali);
1109     tree->Branch("groupnumber",&group,"groupnumber/I");
1110     tree->Branch("histo","TH1F",&histsum,32000,0);
1111
1112     // Init histsum
1113     if (treeChain->GetEntries() < 1) {
1114       return tree1; 
1115     }
1116     
1117     for (Int_t h = 0; h < (Int_t) vectorplace->GetEntriesFast(); h++) {
1118       group = ((AliTRDPlace *) vectorplace->At(h))->GetPlace();
1119       TObjArray *chainplace = ((TObjArray *) where->At(h));
1120       treeChain->GetEntry(((AliTRDPlace *) chainplace->At(0))->GetPlace());
1121       //Init for the first time
1122       if (h == 0)  {
1123         histsum = new TH1F("","",his->GetXaxis()->GetNbins()
1124                                 ,his->GetXaxis()->GetBinLowEdge(1)
1125                                 ,his->GetXaxis()->GetBinUpEdge(his->GetXaxis()->GetNbins()));
1126         histsum->Sumw2();
1127       }
1128       // Reset for each new group
1129       histsum->SetEntries(0.0);
1130       for (Int_t l = 0; l <= histsum->GetXaxis()->GetNbins(); l++) {
1131         histsum->SetBinContent(l,0.0);
1132         histsum->SetBinError(l,0.0);
1133       }
1134       histsum->Add(his,1);
1135       if ((Int_t) chainplace->GetEntriesFast() > 1) {
1136         for (Int_t s = 1; s < (Int_t) chainplace->GetEntriesFast(); s++) {
1137           treeChain->GetEntry(((AliTRDPlace *) chainplace->At(s))->GetPlace());
1138           histsum->Add(his,1);
1139         }
1140       }
1141       tree->Fill();
1142     }
1143
1144   }
1145   else {
1146
1147     // Ready to read the chain
1148     TGraphErrors *his = 0x0;
1149     treeChain->SetBranchAddress("histo",&his);
1150     
1151     // Initialise the final tree
1152     Int_t         group   = -1;
1153     TGraphErrors *histsum = 0x0;
1154     Double_t     *xref    = 0x0;
1155   
1156     tree = new TTree(variablecali,variablecali);
1157     tree->Branch("groupnumber",&group,"groupnumber/I");
1158     tree->Branch("histo","TGraphErrors",&histsum,32000,0);
1159
1160     // Init histsum
1161     if (treeChain->GetEntries() < 1) {
1162       return tree1; 
1163     }
1164
1165     for (Int_t h = 0; h < (Int_t) vectorplace->GetEntriesFast(); h++) {
1166
1167       group = ((AliTRDPlace *) vectorplace->At(h))->GetPlace();
1168       TObjArray *chainplace = ((TObjArray *) where->At(h));
1169       treeChain->GetEntry(((AliTRDPlace *) chainplace->At(0))->GetPlace());
1170       //Init or reset for a new group
1171       Int_t nbins = his->GetN();
1172       Double_t *x;
1173       x    = new Double_t[nbins];
1174       xref = his->GetX();
1175       Double_t *ex;
1176       ex   = new Double_t[nbins];
1177       Double_t *y;
1178       y    = new Double_t[nbins];
1179       Double_t *ey;
1180       ey   = new Double_t[nbins];
1181      
1182       for (Int_t lo = 0; lo < nbins; lo++) {
1183         x[lo]  = xref[lo];
1184         ex[lo] = 0.0;
1185         y[lo]  = 0.0;
1186         ey[lo] = 0.0;
1187       }
1188       delete histsum;
1189       histsum = new TGraphErrors(nbins,x,y,ex,ey);
1190
1191       // Add the first
1192       histsum = AddProfiles(his,histsum);
1193       if ((Int_t) chainplace->GetEntriesFast() > 1) {
1194         for (Int_t s = 1; s < (Int_t) chainplace->GetEntriesFast(); s++) {
1195           treeChain->GetEntry(((AliTRDPlace *) chainplace->At(s))->GetPlace());
1196           histsum = AddProfiles(his,histsum);
1197         }
1198       }
1199
1200       tree->Fill();
1201
1202     }
1203
1204   }
1205     
1206   return tree;
1207
1208 }
1209
1210 //_____________________________________________________________________________
1211 TGraphErrors *AliTRDCalibraVector::AddProfiles(TGraphErrors *hist1
1212                                              , TGraphErrors *hist2) const
1213 {
1214   //
1215   // In the case of the vectors method we use TGraphErrors for PH and PRF
1216   // to be able to add the them after
1217   // Here we add the TGraphErrors  
1218   //
1219
1220   // First TGraphErrors
1221   Int_t     nbins1 = hist1->GetN();
1222   Double_t *x1     = hist1->GetX();
1223   Double_t *ex1    = hist1->GetEX();
1224   Double_t *y1     = hist1->GetY();
1225   Double_t *ey1    = hist1->GetEY();
1226
1227   TGraphErrors *rehist = new TGraphErrors(nbins1);
1228
1229   // Second TGraphErrors
1230   Double_t *ex2    = hist2->GetEX();
1231   Double_t *y2     = hist2->GetY();
1232   Double_t *ey2    = hist2->GetEY();
1233
1234   // Define the Variables for the new TGraphErrors
1235   Double_t x;
1236   Double_t ex;
1237   Double_t y;
1238   Double_t ey;
1239   
1240   for (Int_t k = 0; k < nbins1; k++) {
1241     Double_t nentries = 0.0;
1242     x  = x1[k];
1243     y  = 0.0;
1244     ey = 0.0;
1245     ex = 0.0;
1246     if ((ex2[k] == 0.0) && 
1247         (ex1[k] == 0.0)) {
1248       nentries = 0.0;
1249     }
1250     if ((ex2[k] == 0.0) && 
1251         (ex1[k]  > 0.0)) {
1252       nentries = ex1[k];
1253       y  = y1[k];
1254       ey = ey1[k];
1255       ex = ex1[k];
1256     }
1257     if ((ex2[k]  > 0.0) && 
1258         (ex1[k] == 0.0)) {
1259       nentries = ex2[k];
1260       y  = y2[k];
1261       ey = ey2[k];
1262       ex = ex2[k];
1263     }
1264     if ((ex2[k] > 0.0) && 
1265         (ex1[k] > 0.0)) { 
1266      nentries = ex1[k] + ex2[k];
1267      y  = ( y1[k]*ex1[k]+ y2[k]*ex2[k]) / nentries;
1268      ey = (ey1[k]*ex1[k]+ey2[k]*ex2[k]) / nentries;
1269      ex = nentries;
1270    }
1271    rehist->SetPoint(k,x,y);
1272    rehist->SetPointError(k,ex,ey);
1273  }
1274
1275  return rehist;
1276
1277 }