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 {
  59   //
  60   // Default constructor
  61   //
  62
  63 }
  64
  65 //______________________________________________________________________________________
  66 AliTRDCalibraVector::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 //____________________________________________________________________________________
  85 AliTRDCalibraVector::~AliTRDCalibraVector()
  86 {
  87   //
  88   // AliTRDCalibraVector destructor
  89   //
  90
  91 }
  92
  93 //_____________________________________________________________________________
  94 Int_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 //_____________________________________________________________________________
 134 Int_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 //_____________________________________________________________________________
 173 Int_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 //_____________________________________________________________________________
 190 Bool_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 //_____________________________________________________________________________
 246 Bool_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 //_____________________________________________________________________________
 333 Bool_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 //_____________________________________________________________________________
 422 TGraphErrors *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 //_____________________________________________________________________________
 440 TGraphErrors *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 //_____________________________________________________________________________
 512 TH1F *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 //_____________________________________________________________________________
 530 TH1F *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 //_____________________________________________________________________________
 556 TTree *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 //_____________________________________________________________________________
 605 TTree *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 //_____________________________________________________________________________
 656 TObjArray *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 //_____________________________________________________________________________
 687 Bool_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 //_____________________________________________________________________________
 767 Bool_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 //_____________________________________________________________________________
 963 TTree *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 //_____________________________________________________________________________
1156 TGraphErrors *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 }