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