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