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 <TTreeStream.h>
  38
  39
  40 //header file
  41 #include "AliTRDCalibraExbAltFit.h"
  42
  43 ClassImp(AliTRDCalibraExbAltFit) /*FOLD00*/
  44
  45 //_____________________________________________________________________
  46 AliTRDCalibraExbAltFit::AliTRDCalibraExbAltFit() : /*FOLD00*/
  47   TObject(),
  48   fVersion(0),
  49   fFitterHistoArray(540),
  50   fFitterPArray(540),
  51   fFitterEArray(540),
  52   fRobustFit(kFALSE)
  53 {
  54   //
  55   // default constructor
  56   //
  57 }
  58 //_____________________________________________________________________
  59 AliTRDCalibraExbAltFit::AliTRDCalibraExbAltFit(const AliTRDCalibraExbAltFit &ped) : /*FOLD00*/
  60   TObject(ped),
  61   fVersion(ped.fVersion),
  62   fFitterHistoArray(540),
  63   fFitterPArray(540),
  64   fFitterEArray(540),
  65   fRobustFit(kFALSE)
  66 {
  67     //
  68     // copy constructor
  69     //
  70   for (Int_t idet = 0; idet < 540; idet++){
  71
  72     const TVectorD     *vectorE     = (TVectorD*)ped.fFitterEArray.UncheckedAt(idet);
  73     const TVectorD     *vectorP     = (TVectorD*)ped.fFitterPArray.UncheckedAt(idet);
  74     const TH2S         *hped        = (TH2S*)ped.fFitterHistoArray.UncheckedAt(idet);
  75
  76     if ( vectorE != 0x0 ) fFitterEArray.AddAt(new TVectorD(*vectorE), idet);
  77     if ( vectorP != 0x0 ) fFitterPArray.AddAt(new TVectorD(*vectorP), idet);
  78     if ( hped != 0x0 ){
  79       TH2S *hNew = (TH2S *)hped->Clone();
  80       //hNew->SetDirectory(0);
  81       fFitterHistoArray.AddAt(hNew,idet);
  82     }
  83   }
  84 }
  85 //_____________________________________________________________________
  86 AliTRDCalibraExbAltFit::AliTRDCalibraExbAltFit(const TObjArray &obja) : /*FOLD00*/
  87   TObject(),
  88   fVersion(0),
  89   fFitterHistoArray(540),
  90   fFitterPArray(540),
  91   fFitterEArray(540),
  92   fRobustFit(kFALSE)
  93 {
  94   //
  95   // constructor from a TObjArray
  96   //
  97   for (Int_t idet = 0; idet < 540; idet++){
  98     const TH2S         *hped        = (TH2S*)obja.UncheckedAt(idet);
  99     if ( hped != 0x0 ){
 100       TH2S *hNew = (TH2S *)hped->Clone();
 101       //hNew->SetDirectory(0);
 102       fFitterHistoArray.AddAt(hNew,idet);
 103     }
 104   }
 105 }
 106 //_____________________________________________________________________
 107 AliTRDCalibraExbAltFit& AliTRDCalibraExbAltFit::operator = (const  AliTRDCalibraExbAltFit &source)
 108 {
 109   //
 110   // assignment operator
 111   //
 112   if (&source == this) return *this;
 113   new (this) AliTRDCalibraExbAltFit(source);
 114
 115   return *this;
 116 }
 117 //_____________________________________________________________________
 118 AliTRDCalibraExbAltFit::~AliTRDCalibraExbAltFit() /*FOLD00*/
 119 {
 120   //
 121   // destructor
 122   //
 123   fFitterHistoArray.SetOwner();
 124   fFitterPArray.SetOwner();
 125   fFitterEArray.SetOwner();
 126
 127   fFitterHistoArray.Delete();
 128   fFitterPArray.Delete();
 129   fFitterEArray.Delete();
 130
 131 }
 132 //_____________________________________________________________________________
 133 void AliTRDCalibraExbAltFit::Copy(TObject &c) const
 134 {
 135   //
 136   // Copy function
 137   //
 138
 139   AliTRDCalibraExbAltFit& target = (AliTRDCalibraExbAltFit &) c;
 140
 141   // Copy only the histos
 142   for (Int_t idet = 0; idet < 540; idet++){
 143     if(fFitterHistoArray.UncheckedAt(idet)){
 144       TH2S *hped1 = (TH2S *)target.GetFitterHisto(idet,kTRUE);
 145       //hped1->SetDirectory(0);
 146       hped1->Add((const TH2S *)fFitterHistoArray.UncheckedAt(idet));
 147     }
 148   }
 149
 150   TObject::Copy(c);
 151
 152 }
 153 //_____________________________________________________________________________
 154 Long64_t AliTRDCalibraExbAltFit::Merge(const TCollection* list)
 155 {
 156   // Merge list of objects (needed by PROOF)
 157
 158   if (!list)
 159     return 0;
 160
 161   if (list->IsEmpty())
 162     return 1;
 163
 164   TIterator* iter = list->MakeIterator();
 165   TObject* obj = 0;
 166
 167   // collection of generated histograms
 168   Int_t count=0;
 169   while((obj = iter->Next()) != 0)
 170     {
 171       AliTRDCalibraExbAltFit* entry = dynamic_cast<AliTRDCalibraExbAltFit*>(obj);
 172       if (entry == 0) continue;
 173
 174       // Copy only the histos
 175       for (Int_t idet = 0; idet < 540; idet++){
 176         if(entry->GetFitterHisto(idet)){
 177           TH2S *hped1 = (TH2S *)GetFitterHisto(idet,kTRUE);
 178           Double_t entriesa = hped1->GetEntries();
 179           Double_t entriesb = ((TH2S *)entry->GetFitterHisto(idet))->GetEntries();
 180           if((entriesa + entriesb) < 5*32767) hped1->Add(entry->GetFitterHisto(idet));
 181         }
 182       }
 183
 184       count++;
 185     }
 186
 187   return count;
 188 }
 189 //_____________________________________________________________________
 190 void AliTRDCalibraExbAltFit::Add(const AliTRDCalibraExbAltFit *ped)
 191 {
 192   //
 193   // Add histo
 194   //
 195
 196   fVersion++;
 197
 198   for (Int_t idet = 0; idet < 540; idet++){
 199     const TH2S         *hped        = (TH2S*)ped->GetFitterHistoNoForce(idet);
 200     //printf("idet %d\n",idet);
 201     if ( hped != 0x0 ){
 202       //printf("add\n");
 203       TH2S *hped1 = (TH2S *)GetFitterHisto(idet,kTRUE);
 204       Double_t entriesa = hped1->GetEntries();
 205       Double_t entriesb = hped->GetEntries();
 206       if((entriesa + entriesb) < 5*32767) hped1->Add(hped);
 207     }
 208   }
 209 }
 210 //______________________________________________________________________________________
 211 TH2S* AliTRDCalibraExbAltFit::GetFitterHisto(Int_t detector, Bool_t force)
 212 {
 213     //
 214     // return pointer to TH2F histo
 215     // if force is true create a new histo if it doesn't exist allready
 216     //
 217     if ( !force || fFitterHistoArray.UncheckedAt(detector) )
 218         return (TH2S*)fFitterHistoArray.UncheckedAt(detector);
 219
 220     return GetFitterHistoForce(detector);
 221
 222 }
 223 //______________________________________________________________________________________
 224 TH2S* AliTRDCalibraExbAltFit::GetFitterHistoForce(Int_t detector)
 225 {
 226   //
 227   // return pointer to TH2F histo
 228   // if NULL create a new histo if it doesn't exist allready
 229   //
 230   if (fFitterHistoArray.UncheckedAt(detector))
 231     return (TH2S*)fFitterHistoArray.UncheckedAt(detector);
 232
 233   // if we are forced and TLinearFitter doesn't yes exist create it
 234
 235   // new TH2F
 236   TString name("LFEXB");
 237   name += detector;
 238   name += "version";
 239   name +=  fVersion;
 240
 241   TH2S *lfdv = new TH2S((const Char_t *)name,(const Char_t *) name,
 242                         30, -TMath::DegToRad()*45, TMath::DegToRad()*45,
 243                         30, 0.3, 1.4);
 244   lfdv->SetXTitle("tan(phi_{track})");
 245   lfdv->SetYTitle("rms");
 246   lfdv->SetZTitle("Number of tracklets");
 247   lfdv->SetStats(0);
 248   lfdv->SetDirectory(0);
 249
 250   fFitterHistoArray.AddAt(lfdv,detector);
 251   return lfdv;
 252 }
 253 //______________________________________________________________________________________
 254 Bool_t AliTRDCalibraExbAltFit::GetParam(Int_t detector, TVectorD *param)
 255 {
 256     //
 257     // return param for this detector
 258     //
 259   if ( fFitterPArray.UncheckedAt(detector) ){
 260     const TVectorD     *vectorP     = (TVectorD*)fFitterPArray.UncheckedAt(detector);
 261     if(!param) param = new TVectorD(vectorP->GetNoElements());
 262     for(Int_t k = 0; k < vectorP->GetNoElements(); k++){
 263       (*param)[k] = (*vectorP)[k];
 264     }
 265     return kTRUE;
 266   }
 267   else return kFALSE;
 268
 269 }
 270 //______________________________________________________________________________________
 271 Bool_t AliTRDCalibraExbAltFit::GetError(Int_t detector, TVectorD *error)
 272 {
 273     //
 274     // return error for this detector
 275     //
 276   if ( fFitterEArray.UncheckedAt(detector) ){
 277     const TVectorD     *vectorE     = (TVectorD*)fFitterEArray.UncheckedAt(detector);
 278     if(!error) error = new TVectorD(vectorE->GetNoElements());
 279     for(Int_t k = 0; k < vectorE->GetNoElements(); k++){
 280       (*error)[k] = (*vectorE)[k];
 281     }
 282     return kTRUE;
 283   }
 284   else return kFALSE;
 285
 286 }
 287 //______________________________________________________________________________________
 288 void AliTRDCalibraExbAltFit::Update(Int_t detector, Float_t tnp, Float_t pars1)
 289 {
 290     //
 291     // Fill the 2D histos for debugging
 292     //
 293
 294   TH2S *h = ((TH2S *) GetFitterHisto(detector,kTRUE));
 295   Double_t nbentries = h->GetEntries();
 296   if(nbentries < 5*32767) h->Fill(tnp,pars1);
 297
 298 }
 299 //____________Functions fit Online CH2d________________________________________
 300 void AliTRDCalibraExbAltFit::FillPEArray()
 301 {
 302   //
 303   // Fill fFitterPArray and fFitterEArray from inside
 304   //
 305
 306
 307   Int_t *arrayI = new Int_t[540];
 308   for(Int_t k = 0; k< 540; k++){
 309     arrayI[k] = 0;
 310   }
 311
 312   // Loop over histos
 313   for(Int_t cb = 0; cb < 540; cb++){
 314     const TH2S *fitterhisto = (TH2S*)fFitterHistoArray.UncheckedAt(cb);
 315     //printf("Processing the detector cb %d we find %d\n",cb, (Bool_t) fitterhisto);
 316
 317     if ( fitterhisto != 0 ){
 318
 319       // Fill a fitter
 320       TAxis *xaxis = fitterhisto->GetXaxis();
 321       TAxis *yaxis = fitterhisto->GetYaxis();
 322       TLinearFitter fitter = TLinearFitter(3,"pol2");
 323       //printf("test\n");
 324       Double_t integral = fitterhisto->Integral();
 325       //printf("Integral is %f\n",integral);
 326       Bool_t securitybreaking = kFALSE;
 327       if(TMath::Abs(integral-1199) < 0.00001) securitybreaking = kTRUE;
 328       for(Int_t ibinx = 0; ibinx < fitterhisto->GetNbinsX(); ibinx++){
 329         for(Int_t ibiny = 0; ibiny < fitterhisto->GetNbinsY(); ibiny++){
 330           if(fitterhisto->GetBinContent(ibinx+1,ibiny+1)>0){
 331             Double_t x = xaxis->GetBinCenter(ibinx+1);
 332             Double_t y = yaxis->GetBinCenter(ibiny+1);
 333
 334             for(Int_t k = 0; k < (Int_t)fitterhisto->GetBinContent(ibinx+1,ibiny+1); k++){
 335               if(!securitybreaking){
 336                 fitter.AddPoint(&x,y);
 337                 arrayI[cb]++;
 338               }
 339               else {
 340                 if(arrayI[cb]< 1198){
 341                   fitter.AddPoint(&x,y);
 342                   arrayI[cb]++;
 343                 }
 344               }
 345             }
 346
 347           }
 348         }
 349       }
 350
 351       //printf("Find %d entries for the detector %d\n",arrayI[cb],cb);
 352
 353       // Eval the  fitter
 354       if(arrayI[cb]>15){
 355         TVectorD  *par  = new TVectorD(3);
 356         TVectorD  pare = TVectorD(2);
 357         TVectorD  *parE = new TVectorD(3);
 358         //printf("Fit\n");
 359         //if((fitter.EvalRobust(0.8)==0)) {
 360         //if((fitter.EvalRobust()==0)) {
 361         if(((fRobustFit) && (fitter.EvalRobust(0.8)==0)) || ((!fRobustFit) && (fitter.Eval()==0))) {
 362           //if((fitter.Eval()==0)) {
 363           //printf("Take the param\n");
 364           fitter.GetParameters(*par);
 365           //printf("Done\n");
 366           //fitter.GetErrors(*pare);
 367           //Float_t  ppointError =  TMath::Sqrt(TMath::Abs(fitter.GetChisquare())/arrayI[cb]);
 368           //(*parE)[0] = pare[0]*ppointError;
 369           //(*parE)[1] = pare[1]*ppointError;
 370           (*parE)[0] = 0.0;
 371           (*parE)[1] = 0.0;
 372           (*parE)[2] = (Double_t) arrayI[cb];
 373           fFitterPArray.AddAt(par,cb);
 374           fFitterEArray.AddAt(parE,cb);
 375
 376           //par->Print();
 377           //parE->Print();
 378         }
 379         //printf("Finish\n");
 380       }
 381
 382       //delete fitterhisto;
 383
 384     }// if something
 385
 386   }
 387
 388   delete [] arrayI;
 389
 390 }
 391