]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - TRD/AliTRDCalibraVector.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
New DAQ calibration DAs by Raphaelle
[u/mrichter/AliRoot.git] / TRD / AliTRDCalibraVector.cxx
... / ...
CommitLineData
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 ,fPRFRange(1.5)
59{
60 //
61 // Default constructor
62 //
63
64}
65
66//______________________________________________________________________________________
67AliTRDCalibraVector::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//____________________________________________________________________________________
140AliTRDCalibraVector::~AliTRDCalibraVector()
141{
142 //
143 // AliTRDCalibraVector destructor
144 //
145
146}
147
148//_____________________________________________________________________________
149Int_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//_____________________________________________________________________________
189Int_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//_____________________________________________________________________________
228Int_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//_____________________________________________________________________________
245Bool_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//_____________________________________________________________________________
301Bool_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//_____________________________________________________________________________
388Bool_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//_____________________________________________________________________________
477TGraphErrors *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//_____________________________________________________________________________
495TGraphErrors *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//_____________________________________________________________________________
567TH1F *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//_____________________________________________________________________________
585TH1F *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//_____________________________________________________________________________
611TTree *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//_____________________________________________________________________________
660TTree *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//_____________________________________________________________________________
711TObjArray *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//_____________________________________________________________________________
742Bool_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//_____________________________________________________________________________
822Bool_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//_____________________________________________________________________________
1018TTree *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//_____________________________________________________________________________
1211TGraphErrors *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}