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