TRD/AliTRDCalibraVdriftLinearFit.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 // AliTRDCalibraVdriftLinearFit                                           //
  21 //                                                                        //
  22 // Does the Vdrift an ExB calibration by applying a linear fit            //
  23 //                                                                        //
  24 // Author:                                                                //
  25 //   R. Bailhache (R.Bailhache@gsi.de)                                    //
  26 //                                                                        //
  27 ////////////////////////////////////////////////////////////////////////////
  28
  29 //Root includes
  30 #include <TObjArray.h>
  31 #include <TH2F.h>
  32 #include <TString.h>
  33 #include <TVectorD.h>
  34 #include <TAxis.h>
  35 #include <TLinearFitter.h>
  36 #include <TMath.h>
  37 #include <TStyle.h>
  38 #include <TCanvas.h>
  39 #include <TTreeStream.h>
  40 #include <TGraphErrors.h>
  41 #include <TDirectory.h>
  42 #include <TTreeStream.h>
  43 #include <TF1.h>
  44
  45 //header file
  46 #include "AliTRDCalibraVdriftLinearFit.h"
  47
  48 ClassImp(AliTRDCalibraVdriftLinearFit) /*FOLD00*/
  49
  50 //_____________________________________________________________________
  51 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit() : /*FOLD00*/
  52   TObject(),
  53   fVersion(0),
  54   fNameCalibUsed(""),
  55   fLinearFitterHistoArray(540),
  56   fLinearFitterPArray(540),
  57   fLinearFitterEArray(540),
  58   fRobustFit(kTRUE),
  59   fMinNpointsFit(11),
  60   fNbBindx(32),
  61   fNbBindy(70),
  62   fRangedx(0.8),
  63   fRangedy(1.4),
  64   fDebugStreamer(0x0),
  65   fDebugLevel(0),
  66   fSeeDetector(0)
  67 {
  68   //
  69   // default constructor
  70   //
  71 }
  72 //_____________________________________________________________________
  73 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const AliTRDCalibraVdriftLinearFit &ped) : /*FOLD00*/
  74   TObject(ped),
  75   fVersion(ped.fVersion),
  76   fNameCalibUsed(ped.fNameCalibUsed),
  77   fLinearFitterHistoArray(540),
  78   fLinearFitterPArray(540),
  79   fLinearFitterEArray(540),
  80   fRobustFit(kTRUE),
  81   fMinNpointsFit(10),
  82   fNbBindx(32),
  83   fNbBindy(70),
  84   fRangedx(0.8),
  85   fRangedy(1.4),
  86   fDebugStreamer(0x0),
  87   fDebugLevel(0),
  88   fSeeDetector(0)
  89 {
  90     //
  91     // copy constructor
  92     //
  93   for (Int_t idet = 0; idet < 540; idet++){
  94
  95     const TVectorD     *vectorE     = (TVectorD*)ped.fLinearFitterEArray.UncheckedAt(idet);
  96     const TVectorD     *vectorP     = (TVectorD*)ped.fLinearFitterPArray.UncheckedAt(idet);
  97     const TH2S         *hped        = (TH2S*)ped.fLinearFitterHistoArray.UncheckedAt(idet);
  98
  99     if ( vectorE != 0x0 ) fLinearFitterEArray.AddAt(new TVectorD(*vectorE), idet);
 100     if ( vectorP != 0x0 ) fLinearFitterPArray.AddAt(new TVectorD(*vectorP), idet);
 101     if ( hped != 0x0 ){
 102       TH2S *hNew = (TH2S *)hped->Clone();
 103       //hNew->SetDirectory(0);
 104       fLinearFitterHistoArray.AddAt(hNew,idet);
 105     }
 106   }
 107 }
 108 //_____________________________________________________________________
 109 AliTRDCalibraVdriftLinearFit::AliTRDCalibraVdriftLinearFit(const TObjArray &obja) : /*FOLD00*/
 110   TObject(),
 111   fVersion(0),
 112   fNameCalibUsed(""),
 113   fLinearFitterHistoArray(540),
 114   fLinearFitterPArray(540),
 115   fLinearFitterEArray(540),
 116   fRobustFit(kTRUE),
 117   fMinNpointsFit(10),
 118   fNbBindx(32),
 119   fNbBindy(70),
 120   fRangedx(0.8),
 121   fRangedy(1.4),
 122   fDebugStreamer(0x0),
 123   fDebugLevel(0),
 124   fSeeDetector(0)
 125 {
 126   //
 127   // constructor from a TObjArray
 128   //
 129   for (Int_t idet = 0; idet < 540; idet++){
 130     const TH2S         *hped        = (TH2S*)obja.UncheckedAt(idet);
 131     if ( hped != 0x0 ){
 132       TH2S *hNew = (TH2S *)hped->Clone();
 133       //hNew->SetDirectory(0);
 134       fLinearFitterHistoArray.AddAt(hNew,idet);
 135     }
 136   }
 137 }
 138 //_____________________________________________________________________
 139 AliTRDCalibraVdriftLinearFit& AliTRDCalibraVdriftLinearFit::operator = (const  AliTRDCalibraVdriftLinearFit &source)
 140 {
 141   //
 142   // assignment operator
 143   //
 144   if (&source == this) return *this;
 145   new (this) AliTRDCalibraVdriftLinearFit(source);
 146
 147   return *this;
 148 }
 149 //_____________________________________________________________________
 150 AliTRDCalibraVdriftLinearFit::~AliTRDCalibraVdriftLinearFit() /*FOLD00*/
 151 {
 152   //
 153   // destructor
 154   //
 155   fLinearFitterHistoArray.SetOwner();
 156   fLinearFitterPArray.SetOwner();
 157   fLinearFitterEArray.SetOwner();
 158
 159   fLinearFitterHistoArray.Delete();
 160   fLinearFitterPArray.Delete();
 161   fLinearFitterEArray.Delete();
 162
 163   if ( fDebugStreamer ) delete fDebugStreamer;
 164
 165 }
 166 //_____________________________________________________________________________
 167 void AliTRDCalibraVdriftLinearFit::Copy(TObject &c) const
 168 {
 169   //
 170   // Copy function
 171   //
 172
 173   AliTRDCalibraVdriftLinearFit& target = (AliTRDCalibraVdriftLinearFit &) c;
 174
 175   // Copy only the histos
 176   for (Int_t idet = 0; idet < 540; idet++){
 177     if(fLinearFitterHistoArray.UncheckedAt(idet)){
 178       TH2S *hped1 = (TH2S *)target.GetLinearFitterHisto(idet,kTRUE);
 179       //hped1->SetDirectory(0);
 180       hped1->Add((const TH2S *)fLinearFitterHistoArray.UncheckedAt(idet));
 181     }
 182   }
 183
 184   TObject::Copy(c);
 185
 186 }
 187 //_____________________________________________________________________________
 188 Long64_t AliTRDCalibraVdriftLinearFit::Merge(const TCollection* list)
 189 {
 190   // Merge list of objects (needed by PROOF)
 191
 192   if (!list)
 193     return 0;
 194
 195   if (list->IsEmpty())
 196     return 1;
 197
 198   TIterator* iter = list->MakeIterator();
 199   TObject* obj = 0;
 200
 201   // collection of generated histograms
 202   Int_t count=0;
 203   while((obj = iter->Next()) != 0)
 204     {
 205       AliTRDCalibraVdriftLinearFit* entry = dynamic_cast<AliTRDCalibraVdriftLinearFit*>(obj);
 206       if (entry == 0) continue;
 207
 208       // Copy only the histos
 209       for (Int_t idet = 0; idet < 540; idet++){
 210         if(entry->GetLinearFitterHisto(idet)){
 211           TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
 212           Double_t entriesa = hped1->GetEntries();
 213           Double_t entriesb = ((TH2S *)entry->GetLinearFitterHisto(idet))->GetEntries();
 214           if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetLinearFitterHisto(idet));
 215         }
 216       }
 217
 218       count++;
 219     }
 220
 221   return count;
 222 }
 223 //_____________________________________________________________________
 224 void AliTRDCalibraVdriftLinearFit::Add(const AliTRDCalibraVdriftLinearFit *ped)
 225 {
 226   //
 227   // Add histo
 228   //
 229
 230   fVersion++;
 231
 232   for (Int_t idet = 0; idet < 540; idet++){
 233     const TH2S         *hped        = (TH2S*)ped->GetLinearFitterHistoNoForce(idet);
 234     //printf("idet %d\n",idet);
 235     if ( hped != 0x0 ){
 236       //printf("add\n");
 237       TH2S *hped1 = (TH2S *)GetLinearFitterHisto(idet,kTRUE);
 238       Double_t entriesa = hped1->GetEntries();
 239       Double_t entriesb = hped->GetEntries();
 240       if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
 241     }
 242   }
 243 }
 244 //______________________________________________________________________________________
 245 TH2S* AliTRDCalibraVdriftLinearFit::AddAll()
 246 {
 247     //
 248     // return pointer to TH2S of all added histos
 249   //
 250
 251   TH2S *histo = 0x0;
 252   for(Int_t k=0; k < 540; k++) {
 253     TH2S * u = GetLinearFitterHistoForce(k);
 254     if(k == 0) {
 255       histo = new TH2S(*u);
 256       histo->SetName("sum");
 257     }
 258     else histo->Add(u);
 259   }
 260
 261   return histo;
 262
 263 }
 264 //______________________________________________________________________________________
 265 TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHisto(Int_t detector, Bool_t force)
 266 {
 267     //
 268     // return pointer to TH2F histo
 269     // if force is true create a new histo if it doesn't exist allready
 270     //
 271     if ( !force || fLinearFitterHistoArray.UncheckedAt(detector) )
 272         return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
 273
 274     return GetLinearFitterHistoForce(detector);
 275
 276 }
 277 //______________________________________________________________________________________
 278 TH2S* AliTRDCalibraVdriftLinearFit::GetLinearFitterHistoForce(Int_t detector)
 279 {
 280   //
 281   // return pointer to TH2F histo
 282   // if NULL create a new histo if it doesn't exist allready
 283   //
 284   if (fLinearFitterHistoArray.UncheckedAt(detector))
 285     return (TH2S*)fLinearFitterHistoArray.UncheckedAt(detector);
 286
 287   // if we are forced and TLinearFitter doesn't yes exist create it
 288
 289   // new TH2F
 290   TString name("LFDV");
 291   name += detector;
 292   name += "version";
 293   name +=  fVersion;
 294
 295   TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name
 296                         ,(Int_t)fNbBindx,-fRangedx,fRangedx,(Int_t)fNbBindy
 297                         ,-fRangedy,fRangedy);
 298   lfdv->SetXTitle("tan(phi_{track})");
 299   lfdv->SetYTitle("dy/dt");
 300   lfdv->SetZTitle("Number of clusters");
 301   lfdv->SetStats(0);
 302   lfdv->SetDirectory(0);
 303
 304   fLinearFitterHistoArray.AddAt(lfdv,detector);
 305   return lfdv;
 306 }
 307 //______________________________________________________________________________________
 308 Bool_t AliTRDCalibraVdriftLinearFit::GetParam(Int_t detector, TVectorD *param)
 309 {
 310     //
 311     // return param for this detector
 312     //
 313   if ( fLinearFitterPArray.UncheckedAt(detector) ){
 314     const TVectorD     *vectorP     = (TVectorD*)fLinearFitterPArray.UncheckedAt(detector);
 315     if(!param) param = new TVectorD(2);
 316     for(Int_t k = 0; k < 2; k++){
 317       (*param)[k] = (*vectorP)[k];
 318     }
 319     return kTRUE;
 320   }
 321   else return kFALSE;
 322
 323 }
 324 //______________________________________________________________________________________
 325 Bool_t AliTRDCalibraVdriftLinearFit::GetError(Int_t detector, TVectorD *error)
 326 {
 327     //
 328     // return error for this detector
 329     //
 330   if ( fLinearFitterEArray.UncheckedAt(detector) ){
 331     const TVectorD     *vectorE     = (TVectorD*)fLinearFitterEArray.UncheckedAt(detector);
 332     if(!error) error = new TVectorD(3);
 333     for(Int_t k = 0; k < 3; k++){
 334       (*error)[k] = (*vectorE)[k];
 335     }
 336     return kTRUE;
 337   }
 338   else return kFALSE;
 339
 340 }
 341 //______________________________________________________________________________________
 342 void AliTRDCalibraVdriftLinearFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
 343 {
 344     //
 345     // Fill the 2D histos for debugging
 346     //
 347
 348   TH2S *h = ((TH2S *) GetLinearFitterHisto(detector,kTRUE));
 349   Double_t nbentries = h->GetEntries();
 350   if(nbentries < 5*32767) h->Fill(tnp,pars1);
 351
 352 }
 353 //____________Functions fit Online CH2d________________________________________
 354 void AliTRDCalibraVdriftLinearFit::FillPEArray()
 355 {
 356   //
 357   // Fill fLinearFitterPArray and fLinearFitterEArray from inside
 358   //
 359
 360
 361   Int_t *arrayI = new Int_t[540];
 362   for(Int_t k = 0; k< 540; k++){
 363     arrayI[k] = 0;
 364   }
 365
 366   // Loop over histos
 367   for(Int_t cb = 0; cb < 540; cb++){
 368     const TH2S *linearfitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
 369     //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) linearfitterhisto);
 370
 371     if ( linearfitterhisto != 0 ){
 372
 373       // Fill a linearfitter
 374       TAxis *xaxis = linearfitterhisto->GetXaxis();
 375       TAxis *yaxis = linearfitterhisto->GetYaxis();
 376       TLinearFitter linearfitter = TLinearFitter(2,"pol1");
 377       //printf("test\n");
 378       Double_t integral = linearfitterhisto->Integral();
 379       //printf("Integral is %f\n",integral);
 380       Bool_t securitybreaking = kFALSE;
 381       if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
 382       for(Int_t ibinx = 0; ibinx < linearfitterhisto->GetNbinsX(); ibinx++){
 383         for(Int_t ibiny = 0; ibiny < linearfitterhisto->GetNbinsY(); ibiny++){
 384           if(linearfitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
 385             Double_t x = xaxis->GetBinCenter(ibinx+1);
 386             Double_t y = yaxis->GetBinCenter(ibiny+1);
 387
 388             for(Int_t k = 0; k < (Int_t)linearfitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
 389               if(!securitybreaking){
 390                 linearfitter.AddPoint(&x,y);
 391                 arrayI[cb]++;
 392               }
 393               else {
 394                 if(arrayI[cb]< 1198){
 395                   linearfitter.AddPoint(&x,y);
 396                   arrayI[cb]++;
 397                 }
 398               }
 399             }
 400
 401           }
 402         }
 403       }
 404
 405       //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
 406
 407       // Eval the linear fitter
 408       if(arrayI[cb]>10){
 409         TVectorD  *par  = new TVectorD(2);
 410         TVectorD   pare = TVectorD(2);
 411         TVectorD  *parE = new TVectorD(3);
 412         //printf("Fit\n");
 413         if(((fRobustFit) && (linearfitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (linearfitter.Eval()==0))) {
 414           //if((linearfitter.Eval()==0)) {
 415           //printf("Take the param\n");
 416           linearfitter.GetParameters(*par);
 417           //printf("Done\n");
 418           //linearfitter.GetErrors(pare);
 419           //Float_t  ppointError =  TMath::Sqrt(TMath::Abs(linearfitter.GetChisquare())/arrayI[cb]);
 420           //(*parE)[0] = pare[0]*ppointError;
 421           //(*parE)[1] = pare[1]*ppointError;
 422
 423           (*parE)[0] = 0.0;
 424           (*parE)[1] = 0.0;
 425           (*parE)[2] = (Double_t) arrayI[cb];
 426           fLinearFitterPArray.AddAt(par,cb);
 427           fLinearFitterEArray.AddAt(parE,cb);
 428
 429           //par->Print();
 430           //parE->Print();
 431         }
 432         //printf("Finish\n");
 433       }
 434
 435       //delete linearfitterhisto;
 436
 437     }// if something
 438
 439   }
 440
 441   delete [] arrayI;
 442
 443 }
 444
 445 //____________Functions fit Online CH2d________________________________________
 446 void AliTRDCalibraVdriftLinearFit::FillPEArray2()
 447 {
 448   //
 449   // Fill fFitterPArray and fFitterEArray from inside
 450   //
 451
 452   // Loop over histos
 453   TF1 *f1 = new TF1("f1","[0]+[1]*x",-1,1);
 454
 455   for(Int_t cb = 0; cb < 540; cb++){
 456     //const TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
 457     TH2S *fitterhisto = (TH2S*)fLinearFitterHistoArray.UncheckedAt(cb);
 458     //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) fitterhisto);
 459
 460     if ( fitterhisto != 0 ){
 461
 462       Int_t nEntries=0;
 463       TGraphErrors *gg=DrawMS(fitterhisto,nEntries);
 464       if (!gg) continue;
 465       // Number of points of the TGraphErrors
 466       //printf("Number of points %d for detector %d\n",gg->GetN(),cb);
 467       if(gg->GetN() <  fMinNpointsFit) {
 468         if(gg) delete gg;
 469         continue;
 470       }
 471       //printf("det %d, number of points %d, nEntries %d\n",cb,gg->GetN(),nEntries);
 472       gg->Fit(f1,"Q0");
 473
 474       TVectorD  *par  = new TVectorD(2);
 475       TVectorD  *parE = new TVectorD(3);
 476       (*parE)[0] = 0.0;
 477       (*parE)[1] = 0.0;
 478       (*parE)[2] = (Double_t) nEntries;
 479       (*par)[0] = f1->GetParameter(0);
 480       (*par)[1] = f1->GetParameter(1);
 481       fLinearFitterPArray.AddAt(par,cb);
 482       fLinearFitterEArray.AddAt(parE,cb);
 483       //printf("Value %f for detector %d\n",(*par)[1],cb);
 484
 485       if(fDebugLevel==0) {
 486         if(gg) delete gg;
 487       }
 488       else {
 489         if(cb==fSeeDetector) {
 490           gStyle->SetPalette(1);
 491           gStyle->SetOptStat(1111);
 492           gStyle->SetPadBorderMode(0);
 493           gStyle->SetCanvasColor(10);
 494           gStyle->SetPadLeftMargin(0.13);
 495           gStyle->SetPadRightMargin(0.01);
 496           TCanvas *stat = new TCanvas("stat","",50,50,600,800);
 497           stat->cd(1);
 498           fitterhisto->Draw("colz");
 499           gg->Draw("P");
 500           f1->Draw("same");
 501           break;
 502         }
 503       }
 504     }
 505   }
 506   if(fDebugLevel==0) delete f1;
 507 }
 508
 509 //_________Helper function__________________________________________________
 510 TGraphErrors* AliTRDCalibraVdriftLinearFit::DrawMS(const TH2 *const h2, Int_t &nEntries)
 511 {
 512   //
 513   // Debug function
 514   //
 515
 516   TF1 fg("fg", "gaus", -10., 30.);
 517   TGraphErrors *gp = new TGraphErrors();
 518
 519   TAxis *ax(h2->GetXaxis());
 520   TAxis *ay(h2->GetYaxis());
 521   TH1D *h1(NULL);
 522   for(Int_t ipt(0), jpt(1), ig(0); ipt<ax->GetNbins(); ipt++, jpt++){
 523     h1 = h2->ProjectionY("py", jpt, jpt);
 524     fg.SetParameter(1, h1->GetMean());
 525     //Float_t x(ax->GetBinCenter(jpt));
 526     Int_t n=Int_t(h1->Integral(1, h1->GetNbinsX()));
 527     nEntries+=n;
 528     if(n < 15){
 529       //Warning("drawMS()", Form("reject x[%d]=%f on n=%d", jpt, x, n));
 530       continue;
 531     }
 532     h1->Fit(&fg, "WWQ0");
 533     if(fg.GetNDF()<2){
 534       //Warning("drawMS()", Form("reject x[%d]=%f on NDF=%d", jpt, x, fg.GetNDF()));
 535       continue;
 536     }
 537     if(((fg.GetParameter(1)+fg.GetParameter(2)/2)>ay->GetXmax()) || ((fg.GetParameter(1)-fg.GetParameter(2)/2)<ay->GetXmin()) || (TMath::Abs(fg.GetParameter(0))< 0.00001)) continue;
 538     gp->SetPoint(ig, ax->GetBinCenter(jpt), fg.GetParameter(1));
 539     gp->SetPointError(ig, 0, TMath::Sqrt(pow(fg.GetParError(1),2) + (1/pow(fg.GetParameter(0),2))));
 540     ig++;
 541   }
 542   delete h1;
 543   return gp;
 544 }