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