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