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